Code source de PC_Mat2_2015_27

#! /usr/bin/env python3
# -*- coding: utf-8; mode: python -*-
r""" Correction d'un exercice de maths avec Python, issu des annales pour les oraux du concours Centrale-Supélec.

- *Sujet :* PC 1, http://www.concours-centrale-supelec.fr/CentraleSupelec/MultiY/C2015/PC-Mat2-2015-27.pdf

- *Dates :* séances de révisions, vendredi 10 et samedi 11 juin (2016).
- *Auteur :* Lilian Besson, pour le cours d'informatique pour tous en prépa MP au Lycée Lakanal (http://perso.crans.org/besson/infoMP/),
- *Licence :* MIT Licence (http://lbesson.mit-license.org).

-----------------------------------------------------------------------------

"""

from __future__ import print_function, division  # Python 2 compatibility

import matplotlib.pyplot as plt


# %% Question 0.
[docs]def q0(): r""" Préliminaires : - Comme pour tous ces exercices, on commence par **importer** les modules requis, ici uniquement `matplotlib <http://matplotlib.org/users/beginner.html>`_ (pour la question `Q5 <#PC_Mat2_2015_27.q5>`_) : >>> import matplotlib.pyplot as plt """ pass
# %% Question 1.
[docs]def q1(): r""" Question 1. Question purement mathématique. - Les variables aléatoires :math:`X_n` suivent toutes une loi de Bernoulli de paramètre :math:`p = \frac{1}{2}` (:math:`X_n \thicksim \mathcal{B}(\frac{1}{2})`), et sont indépendantes, - Et la variable aléatoire :math:`S_n = \sum_{i=1}^{n} X_i` suit donc une loi binomiale de paramètres :math:`(n, p)`, par définition (:math:`S_n \thicksim \mathcal{Bin}(n, p)`. - Donc on peut utiliser le cours directement : .. math:: \mathbb{E}(S_n) = n \mathbb{E}(X_1) = n p = n \frac{1}{2} = \frac{n}{2}, \\ \mathbb{V}(S_n) = n \mathbb{V}(X_1) = n p (1-p) = n \frac{1}{4} = \frac{n}{4}. """ pass
# %% Question 2.
[docs]def q2(): r""" Question 2. - On va noter :math:`f(n)` la fonction définie par :math:`\lfloor \frac{n}{2} \rfloor` si :math:`n` est impair (si :math:`n = 2 k + 1`, :math:`f(n) = k`), et :math:`\lfloor \frac{n-1}{2} \rfloor = \frac{n}{2} - 1` si :math:`n` est pair (si :math:`n = 2 k`, :math:`f(n) = k - 1`), - :math:`f : \mathbb{N}^{*} \to \mathbb{N}` est bien définie, et vérifie :math:`\forall n > 0, f(n) = \min {i \in \mathbb{N} : i < \frac{n}{2}}`. - En fait, :math:`f(n) = \lfloor \frac{n-1}{2} \rfloor` comme me l'avait fait remarqué un élève. - Ainsi, comme :math:`S_n` est à support discret dans :math:`[0, n] \cap \mathbb{N}`, on peut écrire :math:`\mathbb{P}(S_n < \frac{n}{2})` comme une somme finie : .. math:: \mathbb{P}(S_n < \frac{n}{2}) &= \sum_{i = 0}^{i < \frac{n}{2}} \mathbb{P}(S_n = i) \\ &= \sum_{i = 0}^{f(n)} \mathbb{P}(S_n = i) \\ &= \sum_{i = 0}^{f(n)} {n \choose k} p^{k} (1-p)^{n - k} \\ &= \sum_{i = 0}^{f(n)} {n \choose k} (\frac{1}{2})^{k} (\frac{1}{2})^{n - k} \\ &= (\frac{1}{2})^{n} \sum_{i = 0}^{f(n)} {n \choose k}. - Cette dernière expression semble satisfaisante : c'est une expression sommatoire explicite, et l'énoncé demandait de ne pas chercher à la simplifier. """ pass
# %% Question 3.
[docs]def q3(): r""" Question 3. Ici, pour écrire une fonction pour calculer :math:`{n \choose k}`, `plusieurs approches <https://en.wikipedia.org/wiki/Binomial_coefficient#Binomial_coefficient_in_programming_languages>`_ sont possibles : 1. on écrit une fonction :py:func:`factorial` (de :math:`\mathbb{N} \to \mathbb{N}`) qui calcule :math:`n!`, et on utilise :math:`{n \choose k} = \frac{n!}{k! (n-k)!}`. Ça marche, mais ce sera assez lent, et très redondant en calculs (cf. :py:func:`binom_factorial`). 2. ou on écrit une fonction :py:func:`binom` directement, en se basant sur la propriété élémentaire :math:`{n \choose k} = (\frac{n-k+1}{k}) {n \choose k-1}`. C'est la deuxième approché, itérative et facile à écrire, que je recommande. - Note: on pense à l'optimisation gratuite qui consiste à remplacer :math:`k` par :math:`\min(k, n-k)`. Ça donne quelque chose comme ça : :: def binom(n, k): if k < 0 or k > n: # Cas de base (n k) = 0 return 0 if k == 0 or k == n: # Cas de base (n k) = 1 return 1 k = min(k, n - k) # Utilise la symétrie ! produit = 1 for i in range(k): produit *= (n - i) for i in range(k): produit //= (i + 1) return int(produit) # Force à être entier ! .. seealso:: En fait, on peut aussi utiliser :py:func:`scipy.special.binom` du module :py:mod:`scipy.special`, qui fonctionne très bien (elle renvoit un flottant attention) : >>> from scipy.special import binom as binom2 >>> binom2(10, 2) 45.0 """ pass
[docs]def binom(n, k): r""" Calcule efficacement en :math:`\mathcal{O}(\min(k, n-k))` multiplication et division entières le nombre ``binom(n, k)`` :math:`= {n \choose k}` (et en :math:`\mathcal{O}(1)` en mémoire !). - Exemples : >>> binom(6, -3) 0 >>> binom(6, 0) 1 >>> binom(6, 1) 6 >>> binom(6, 2) 15 >>> binom(6, 3) 20 >>> binom(6, 4) 15 >>> binom(6, 5) 6 >>> binom(6, 6) 1 """ if k < 0 or k > n: # Cas de base (n k) = 0 return 0 if k == 0 or k == n: # Cas de base (n k) = 1 return 1 k = min(k, n - k) # Utilise la symétrie ! produit = 1 for i in range(k): produit *= (n - i) for i in range(k): produit //= (i + 1) return int(produit) # Force à être entier !
from math import factorial
[docs]def binom_factorial(n, k): r""" Coefficient binomial calculé en utilisant des factorielles (via :py:func:`math.factorial`) : .. math: {n \choose k} = \frac{n!}{k! (n-k)!}. """ return factorial(n) // (factorial(k) * factorial(n - k))
# %% Question 4.
[docs]def q4(): r""" Question 4. - La question 2 donne la formule, cf. `Q2 <#PC_Mat2_2015_27.q2>`_, - Et la question précédente donne une fonction pour calcul :math:`{n \choose k}`, :py:func:`binom`, - Il s'agit donc juste de combiner les deux : :: def sn2(n): n2 = ((n - 1) // 2) # f(n) somme = sum(binom(n, k) for k in range(1, n2 + 1)) return somme / 2**n .. hint:: Pour éviter des problèmes de conversions entre entiers (``int``, en précision arbitraire) et flottants (``float``, en précision à 16 chiffres), il vaut mieux calculer ``somme / 2**n`` (division entre deux gros entiers) que ``(1/2)**n * somme`` (petit flottant fois gros entier... = 0 !). """ pass
[docs]def sn2(n): r""" Calcule la formule obtenue en `Q2 <#PC_Mat2_2015_27.q2>`_, en sommant les valeurs de :math:`{n \choose k}` pour :math:`k = 1, \dots, f(n)` (en utilisant :py:func:`binom`). - Exemples : >>> sn2(1) 0.0 >>> sn2(2) 0.0 >>> sn2(3) 0.375 >>> sn2(4) 0.25 >>> sn2(5) 0.46875 >>> sn2(50) # doctest: +ELLIPSIS 0.44386... >>> sn2(500) # doctest: +ELLIPSIS 0.48216... >>> sn2(1000) # doctest: +ELLIPSIS 0.48738... >>> sn2(2000) # Ca commence à prendre du temps # doctest: +ELLIPSIS 0.49108... """ # n2 = (n // 2) if n % 2 == 1 else (n // 2) - 1 # f(n) n2 = ((n - 1) // 2) # f(n) # somme = sum(binom(n, k) for k in range(1, n2 + 1)) somme = sum(binom_factorial(n, k) for k in range(1, n2 + 1)) # return (1 / 2)**n * somme # FIXED Marche pas, cf. explication plus haut return somme / 2**n
# %% Question 5.
[docs]def q5(): r""" Question 5. - Avec `Q4 <#PC_Mat2_2015_27.q4>`_, on dispose de la fonction :py:func:`sn2`, et donc on peut calculer :math:`u_{20k} =` ``sn2(20k)``, très facilement. - On crée une liste de valeurs de :math:`u_{20k}` pour :math:`k = 1,\dots,100`, - Enfin, on utilise simplement la fonction :py:func:`matplotlib.pyplot.plot` (``plt.plot(X, Y)``) pour afficher ces points : :: # On créé deux listes valeurs_k = [ 20 * k for k in range(1, 100 + 1) ] valeurs_u = [ sn2(k) for k in valeurs_k ] plt.figure() # Nouvelle figure plt.plot(valeurs_k, valeurs_u) # On affiche les points (k, sn2(20*k)) plt.title("Valeurs de u(20k) pour k = 1 .. 100") plt.show() # On affiche la figure - Cela devrait donner une figure comme ça : .. image:: PC_Mat2_2015_27__figure1.png - On peut donc raisonnablement conjecturer que la suite converge vers :math:`\frac{1}{2}`, ce qu'on va monter en `Q7 <#PC_Mat2_2015_27.q7>`_ et `Q8 <#PC_Mat2_2015_27.q8>`_ (de deux façons différentes), mais la convergence semble très lente. """ # On créé deux listes valeurs_k = [20 * k for k in range(1, 100 + 1)] valeurs_u = [sn2(k) for k in valeurs_k] plt.figure() # Nouvelle figure # plt.hold(True) # On affiche les points (k, sn2(20*k)) plt.plot(valeurs_k, valeurs_u, 'b+-', label=r'$u_{20k}$') # plt.plot(valeurs_k, [0.50] * len(valeurs_k), 'r:', label=r'Limite $\frac{1}{2}$') plt.title(r"Valeurs de $u_{20k}$ pour $k = 1,...,100$") # plt.legend() plt.savefig("PC_Mat2_2015_27__figure1.png") plt.show() # On affiche la figure
# %% Question 6.
[docs]def q6(): r""" Question 6. - Par symétrie, on obtient directement que :math:`\mathbb{P}(S > \frac{n}{2}) = \mathbb{P}(S < \frac{n}{2}) = u_n`. - Soit :math:`n \geq 1`. De façon immédiate, on peut décomposer l'intervalle entier :math:`[0, n] \cap \mathbb{N}` en trois morceaux disjoints : .. math:: [0, n] \cap \mathbb{N} &= { i \in \mathbb{N} : 0 \leq i \leq n} \\ &= \{ 0 \leq i < \frac{n}{2} \} \cup \{ i = \frac{n}{2} \} \cup \{ \frac{n}{2} < i \leq n \}. - Donc en considérant les probabilités pour la v.a. :math:`S_n`, on a directement (par la formule de l'union, cas disjoint) : .. math:: 1 &= \mathbb{P}(S \in [0, n] \cap \mathbb{N}) \\ &= \mathbb{P}(S \in \{ 0 \leq i < \frac{n}{2} \} \cup \{ i = \frac{n}{2} \} \cup \{ \frac{n}{2} < i \leq n \}) \\ &= \mathbb{P}(S \in \{ 0 \leq i < \frac{n}{2} \}) + \mathbb{P}(S \in \{ i = \frac{n}{2} \}) + \mathbb{P}(S \in \{ \frac{n}{2} < i \leq n \}) \\ &= \mathbb{P}(S < \frac{n}{2}) + \mathbb{P}(S_n = \frac{n}{2}) + \mathbb{P}(S > \frac{n}{2}) \\ &= u_n + \mathbb{P}(S_n = \frac{n}{2}) + \mathbb{P}(S > \frac{n}{2}). - Et de même, on procède pareil pour montrer, trivialement, que : .. math:: v_n &= \mathbb{P}(S \geq \frac{n}{2}) \\ &= \mathbb{P}(S_n = \frac{n}{2}) + \mathbb{P}(S > \frac{n}{2}) \\ &= \mathbb{P}(S_n = \frac{n}{2}) + \mathbb{P}(S < \frac{n}{2}) \\ &= \mathbb{P}(S_n = \frac{n}{2}) + u_n. - **Conclusion :** on a bien montré que :math:`\forall n \geq 1, v_n = u_n + \mathbb{P}(S_n = \frac{n}{2})`. """ pass
# %% Question 7.
[docs]def q7(): r""" Question 7. - D'après les calculs faits en `Q6 <#PC_Mat2_2015_27.q6>`_, :math:`u_n \to \frac{1}{2}` *ssi* :math:`\mathbb{P}(S_n = \frac{n}{2}) \to 0` pour :math:`n \to +\infty`. - Mais, pour :math:`n` impair, on a déjà que :math:`\frac{n}{2} \not\in \mathbb{N}` donc :math:`\mathbb{P}(S_n = \frac{n}{2}) = 0`. - Et pour :math:`n = 2k` pair, on peut écrire exactement la valeur de :math:`\mathbb{P}(S_n = \frac{n}{2}) = \mathbb{P}(S_n = k) = (\frac{1}{2})^n {n \choose k}`. - On étudie la limite de :math:`\frac{1}{4^k}{2k \choose k}`, et sans même avoir besoin d'un équivalent (ie. formule de Stirling) on devrait (sa)voir que ça tend vers :math:`0`. Mais soyons rigoureux : .. note:: *Preuve :* (avec Stirling). On a : .. math:: {2k \choose k} &= \frac{(2k)!}{(k!)^2} \\ &\thicksim \frac{\sqrt{2 \pi 2k}(\frac{2k}{\mathrm{e}})^{2k}}{(\sqrt{2 \pi k}(\frac{k}{\mathrm{e}})^{k})^2} \\ &\thicksim \frac{2\sqrt{\pi k} 2^{2k} (\frac{k}{\mathrm{e}})^{2k}}{2 \pi k(\frac{k}{\mathrm{e}})^{2k}} \\ &\thicksim \frac{\sqrt{\pi k} 4^k}{\pi k} \\ &\thicksim \frac{4^k}{\sqrt{\pi k}}. Et donc :math:`\frac{1}{4^k}{2k \choose k} \thicksim \frac{1}{\sqrt{\pi k}}` pour :math:`k \to +\infty`, ce qui prouve que :math:`\mathbb{P}(S_{2k} = \frac{2k}{2}) \to 0` aussi. - **Conclusion :** que :math:`n` soit pair ou impair, on a bien montré que :math:`\mathbb{P}(S_n = \frac{n}{2}) \to 0` pour :math:`n \to +\infty`. - Et donc, cela prouve le résultat observé empiriquement en `Q5 <#PC_Mat2_2015_27.q5>`_ : .. math:: \lim_{n \to +\infty} u_n = \frac{1}{2}. .. note:: On a une estimation de la vitesse de convergence de :math:`\mathbb{P}(S_n = \frac{n}{2})` vers :math:`0` (en :math:`\mathcal{O}(\frac{1}{\sqrt{n}}`), donc cela donne directement un ordre de grandeur de convergence de :math:`u_n` vers sa limite ! Cool non ? En fait, c'est une conséquence immédiate du `théorème central limite <https://fr.wikipedia.org/wiki/Théorème_central_limite>`_, qui est LE résultat central en statistiques (cf. `ce point là <https://fr.wikipedia.org/wiki/Th%C3%A9or%C3%A8me_central_limite#Convergence_vers_la_limite>`_). """ pass
# %% Question 8.
[docs]def q8(): r""" Question 8. .. hint:: On remarque que l'énoncé nous aide en suggérant ici de se restreindre aux :math:`n` pairs, ce qui est exactement ce qu'on a fait en `Q7 <#PC_Mat2_2015_27.q7>`_. - De même qu'au dessus, on obtient que : .. math:: w_n := \mathbb{P}(S_{2n} = n) = \frac{1}{4^n}{2n \choose n} - Donc on peut calculer explicitement :math:`\ln(\frac{w_{n+1}}{w_n})` : .. math:: \ln(\frac{w_{n+1}}{w_n}) &= \ln( \frac{1}{4} {2(n+1) \choose (n+1)}/{2n \choose n} ) \\ &= - \ln(4) + \ln( \frac{(2(n+1))!}{((n+1)!)^2} / \frac{(2n)!}{(n!)^2} ) \\ &= - \ln(4) + \ln( \frac{(2n+2)(2n+1)}{(n+1)^2} ) \\ &= - 2\ln(2) + \ln( 2\frac{2n+1}{n+1} ) \\ &= - 2\ln(2) + \ln(2) + \ln(\frac{2n+1}{n+1}) \\ &= - \ln(2) + \ln(\frac{(2n+1)}{n+1}) \\ &= \ln(\frac{(2n+1)}{2(n+1)}) \\ &= -\ln(\frac{2n+2}{2n+1}) = -\ln(\frac{(2n+1) + 1}{2n+1}) \\ &= -\ln(1 + \frac{1}{2n+1}). - Ensuite, en notant :math:`x = \frac{1}{2n+1}`, :math:`x \to 0` si :math:`n \to \infty`, et comme :math:`\ln(1 + x) \thicksim x`, on obtient : .. math:: \ln(\frac{w_{n+1}}{w_n}) &= -\ln(1 + \frac{1}{2n+1}) \\ &\thicksim - \frac{1}{2n+1} \to 0. - Dès lors, on peut sommer la série de terme général :math:`\ln(\frac{w_{n+1}}{w_n})`, pour obtenir, par comparaison avec la série divergente de terme général de signe constant (négatif) :math:`- \frac{1}{2n+1}`, une série divergente vers :math:`-\infty`. - Mais on reconnaît une somme téléscopique : :math:`\sum_{n=0}^{N} \ln(\frac{w_{n+1}}{w_n}) = \ln(w_{N+1}) - \ln(w_{0})`, donc ça montre que :math:`\ln(w_{N+1}) \to -\infty` pour :math:`N \to +\infty`. - Ainsi on obtient directement que :math:`w_n \to 0` pour :math:`n \to +\infty`, c'est ce qui était demandé. - Enfin on obtient comme voulu que :math:`u_n \to \frac{1}{2}`. .. note:: Ici aussi, cet équivalent :math:`\ln(\frac{w_{n+1}}{w_n}) \thicksim - \frac{1}{2n+1}` donne une indication sur la vitesse de convergence de :math:`w_n \to 0` (en :math:`\mathcal{O}(\sqrt{n})`), qui est la même que celle obtenue en `Q7 <#PC_Mat2_2015_27.q7>`_. Cool. """ pass
# %% Fin de la correction. if __name__ == '__main__': from doctest import testmod print("Test automatique de toutes les doctests ecrites dans la documentation (docstring) de chaque fonction :") testmod(verbose=True) # testmod() print("\nPlus de details sur ces doctests peut etre trouve dans la documentation de Python:\nhttps://docs.python.org/3/library/doctest.html (en anglais)") # Fin de PC_Mat2_2015_27.py