TP 4 - Probabilités, simulation de variables aléatoires.

import numpy as np
import matplotlib.pyplot as plt

Pour simuler des nombres aléatoires on peut utiliser le module scipy.stats, mais aussi numpy.random pour certaines lois. On trouvera dans ce dernier tout un tas de méthodes pour les différentes lois (densité (pdf), fonction de répartition (cdf),...)

Les documentations pour ses deux modules sont disponibles en ligne

• (scipy.stats) https://docs.scipy.org/doc/scipy/reference/stats.html
• (numpy.random) https://docs.scipy.org/doc/numpy/reference/routines.random.html

I) Avec scipy.stats

import scipy.stats as st

Avec ce module, on peut définir une loi st.dist comme un objet python, où dist correspond à une certaine loi.

Par exemple, pour une loi Normale $N(m,\sigma)$ on écrira

m, sig = 0., 1.
X = st.norm(m, sig)

type(X)

scipy.stats._distn_infrastructure.rv_continuous_frozen

Ici $X\sim N(m,\sigma)$.

1) Simulation d'un échantillon de $X$

Pour génrer un échantillon de $X$ (un ensemble de v.a iid de même loi que $X$) on utilisera la méthode .rvs(n), où n est la taille de l'échantillon souhaité.

n = 10
Xn =  X.rvs(n)

Xn

array([ 0.15760552,  0.05140919, -0.61400417, -0.04234742, -0.03610134,
       -2.00245731, -1.12421968,  1.0838306 , -0.19920081, -0.88713999])

Ces nombres constituent les réalisations d'un échantillon de $X$.

Xn =  X.rvs((5,3))
Xn

array([[-0.91312213, -0.61684394, -0.75905904],
       [-0.855824  ,  0.21172854, -0.78341573],
       [ 1.38459715, -0.54957006,  0.15892713],
       [ 0.69717903,  0.24113369,  0.77728617],
       [-0.40141927, -0.02478887, -0.89587404]])

On a aussi accès à la densité de $X$, mais aussi à sa fonction de répartition ou sa fonction quantile.

2) Densité de $X$

Pour obtenir la densité de la loi de $X$ on utilise la méthode .pdf(x) (pour probability density function) où x est le point auquelle on évalue la densité de $X$.

x = np.linspace(-5, 5, 100)
y = X.pdf(x)

plt.plot(x, y)
plt.show()

3) Fonction de répartition de $X$

Pour obtenir la fonction de répartition de la loi de $X$ on utilise la méthode .cdf(x) (pour cumulative distribution function) où x est le point auquel on évalue la fonction de répartition de $X$.

x = np.linspace(-5, 5, 100)
y = X.cdf(x)

plt.plot(x, y)
plt.show()

4) Fonction quantile de $X$

La fonction quantile, définie sur l'intervalle $[0,1]$, est en quelque sorte la fonction réciproque de la fonction de répartition. Autrement dit, Si $F$ est la fonction de répartition d'une v.a $X$, la fonction quantile, qu'on notera $G$ vérifie

$$ F(G(u)) = u \qquad \text{et}\qquad G(F(x))=x \qquad\text{pour}\quad x\in \mathbb{R}\quad\text{et}\quad u\in [0,1]. $$

Pour obtenir la fonction quantile de la loi de $X$ on utilise la méthode .ppf(y) (pour purcent point function) où y est le point auquel on évalue la fonction quantile de $X$.

u = np.linspace(0, 1, 100)
y = X.ppf(u)

plt.plot(u, y)
plt.show()

On peut illustrer graphiquement que la fonction quantile est la fonction réciproque de la function de répartition.

u = np.linspace(0, 1, 100)
v = X.cdf(X.ppf(u))

plt.plot(u, v)
plt.show()

x = np.linspace(-5, 5, 100)
y = X.ppf(X.cdf(x))

plt.plot(x, y)
plt.show()

5) Histogramme

La définition d'un histogramme est disponible à l'adresse suivante

https://fr.wikipedia.org/wiki/Histogramme

En gros: on définie un certain nombre de classes/d'intervalles (on découpe l'ensemble des valeurs que peut prendre un variable aléatoire); pour notre loi normale on prendra

x = np.linspace(-5, 5, 100)

on découpe l'intervalle $[-5,5]$ en 100 parties égales; et à partir de l'échatillon de $X$

n = 1000
Xn = X.rvs(n)

on compte le nombre d'observation qui tombe dans chaque intervalle

plt.hist(Xn, bins=x)
plt.show()

On peut normaliser cette histogramme pour qu'il ressemble à la densité de la loi de $X$

n = 5000
Xn = X.rvs(n)

plt.hist(Xn, bins=x, density=True, label='histogramme')
plt.plot(x, X.pdf(x), lw=3, label='densité')
plt.legend()
plt.show()

6) D'autres propriétés

On peut extraire d'autres informations, comme la moyen (.mean()), la variance (.var()), ou l'écart type (.std() pour standard deviation).

X.mean()

0.0

X.var()

1.0

X.std()

1.0

II) Avec numpy.random

Ce module permet aussi de simuler des nombres aléatoires. Par exemple, la loi uniforme sur $]0,1[$

import numpy.random as npr

n = 10
npr.rand(n)

array([0.808454  , 0.80982503, 0.29787396, 0.50498299, 0.15302792,
       0.69448877, 0.29631666, 0.21776241, 0.91195552, 0.06281807])

npr.rand(5,3)

array([[0.94744846, 0.66979692, 0.52266157],
       [0.52191044, 0.29328341, 0.34686021],
       [0.51081151, 0.40679886, 0.82406177],
       [0.16357571, 0.13403762, 0.71481972],
       [0.00801446, 0.72985478, 0.86807771]])

la loi $N(0,1)$

n = 10
npr.randn(n)

array([ 0.42408263,  1.74279721, -0.59219621, -0.69633437,  1.22445549,
       -3.30614026,  0.26835891, -1.82431504, -0.47148839, -0.5080102 ])

la loi uniforme sur les entiers de $p$ à $q$

p = 2 
q = 5
n = 10

npr.randint(p, q, n)

array([3, 4, 2, 4, 4, 4, 2, 3, 2, 2])

où la loi uniforme sur un ensemble donné

S = ['a', 'b', 'c', 'd']

n = 10

npr.choice(S, 10)

array(['c', 'c', 'c', 'd', 'a', 'a', 'b', 'c', 'a', 'b'], dtype='<U1')

III) Quelques propriétés

On vient de simuler des nombres aléatoires de loi normale, mais comment sont ils produit? En général, la brique de base est la simulation de nombres aléatoires de loi uniforme sur $]0,1[$, via un générateur de nombres pseudo-aléatoires. Ce générateur se base sur une suite numérique non aléatoire, mais qui a un comportement tellement complexe, qu'on a l'impression que les nombres qu'elle donne sont tirés aléatoirement.

1) Simulation de la loi normale

• Génerer deux échantillons $(U_1,\dots, U_n)$ et $(V_1,\dots,V_n)$ de loi uniform sur $]0,1[$ de taille $n=5000$.
• Calculer $$ X_j = \sqrt{-2\log(U_j)}\, \cos(2\pi V_j) \qquad j=1,\dots, n $$
• Comparer l'histogramme normalisé de l'échantillon $(X_1,\dots, X_n)$ à la densité de la loi $N(0,1)$.

n = 5000
Un, Vn = npr.rand(n), npr.rand(n)

Xn = np.sqrt(-2*np.log(Un)) * np.cos(2*np.pi*Vn)

plt.hist(Xn, bins=x, density=True, label='histogramme')
plt.plot(x, X.pdf(x), lw=3, label='densité')
plt.legend()
plt.show()

2)

• Génerer deux échantillons $U=(U_1,\dots, U_n)$ et $V=(V_1,\dots,V_n)$ de loi $N(0,1)$ de taille $n=5000$. Comparer l'histogramme de l'échantillon $U+V$ à celui d'une loi normale dont on précisera les paramètres. Répeter l'opération pour 3 et 4 échantillons de la loi $N(0,1)$.

n = 5000
U, V = npr.randn(n), npr.randn(n)
W = U + V

x = np.linspace(-5, 5, 100)
m, sig = 0, np.sqrt(2)

plt.hist(W, bins=x, density=True, label='histogramme')
plt.plot(x, st.norm(m,sig).pdf(x), lw=3,  label='density N(%s, %f)' % (m, sig))

plt.legend()
plt.show()

On voit que l'histogramme normé de l'échantillon $W$ correspond à la densité de la loi $N(0,2)$. Ceci nous dit que la loi de $W$ est proche de la loi $N(0,2)$.

• Génerer un échantillon $(X_1,\dots, X_n)$ de loi $N(0,1)$ de taille $n=20000$, et évaluer $Y=(F(X_1),\dots, F(X_n))$ où $F$ est la fonction de répartition de la loi $N(0,1)$. Comparer lhistogramme normalisée à la densité de la loi uniforme sur $]0,1[$.

n = 20000
X = npr.randn(n)
Y = st.norm(0,1).cdf(X)

plt.hist(Y, bins=80, density=True, label='histogramme')
y = np.linspace(0, 1, 100)
plt.plot(y, st.uniform(0,1).pdf(y), lw=3, label='densité $\mathcal{U}(0,1)$')
plt.legend()
plt.show()

On voit que l'histogramme normé de l'échantillon $F(X)$ correspond à la densité de la loi $\mathcal{U}(0,1)$. Ceci nous dit que la loi de $F(X)$, où $X\sim N(0,1)$, est (proche de) la loi $\mathcal{U}(0,1)$.

• Génerer un échantillon $U=(U_1,\dots, U_n)$ de loi uniform sur $]0,1[$ de taille $n=5000$, et évaluer $V=-\log(U)$. Comparer l'hitogramme de l'échantillon $V$ à la desnité de la loi exponentielle de paramètre $1$.

n = 5000
U = npr.rand(n)
V = -np.log(U)

v = np.linspace(-1, 6, 100)
plt.hist(V, bins=v, density=True, label='histogramme')
plt.plot(v, st.expon(0,1).pdf(v), lw=3, label='densité $\mathcal{E}(1)$')
plt.legend()
plt.show()

On voit que l'histogramme normé de l'échantillon $V$ correspond à la densité de la loi $\mathcal{E}(1)$. Ceci nous dit que la loi de $V$ est (proche de) la loi $\mathcal{E}(1)$.

3)

Simuler un échantillon $(X_1,\dots, X_n)$ de $X$.

Illustrer l'évolution en fonction de $n$ de la moyenne

$$ M_n = \frac{1}{n} \sum_{j=1}^n X_j,$$

et tracer sur le même graphique la ligne vertical d'ordonnée $y=\mathbb{E}[X]$.

Recommencer plusieur. fois Que peut on observer? A-t-on un phénomène similaire si on change de loi pour l'échantillon de $X$. Avez vous une idée de la raison de ce phénomène?

n = 1000
Xn = X.rvs(n)

Mn = np.cumsum(Xn) / np.arange(1,n+1)

plt.plot(np.arange(1,n+1), Mn, label="Mn")
plt.hlines(X.mean(), 0, n+1, color='orange', label="$\mathbb{E}[X]$")
plt.legend()
plt.show()

On peut observer que quand $n$ devient grand, $M_n$ est plus stable tout en se rapprochant de $\mathbb{E}[X]$.

n = 1000
N = 10

X = st.norm(0,1)
#X = st.expon(1)
#X = st.gamma(1)
#X = st.cauchy(1)

Xn = X.rvs((N,n))

Mn = np.cumsum(Xn, axis = 1) / np.arange(1,n+1)

cmap = plt.get_cmap("tab10")

for j in range(N-1):
    plt.plot(np.arange(1,n+1), Mn[j], color=cmap(0))
plt.plot(np.arange(1,n+1), Mn[N-1], color=cmap(0), label="Mn")
    
plt.hlines(X.mean(), 0, n+1, color='orange', lw=3, label="$\mathbb{E}[X]$")
plt.legend()
plt.show()

On observe toujour le même phénomène en changeant la loi, sauf pour la loi de Cauchy. La raison est que cette loi n'admet pas d'espérance.

Si on calcul l'espérance de $M_n$ on a

$$ \mathbb{E}[M_n] = \mathbb{E}\Big[\frac{1}{n}\sum_{j=1}^n X_j \Big] = \frac{1}{n}\sum_{j=1}^n\mathbb{E}[X_j] = \frac{1}{n}\sum_{j=1}^n\mathbb{E}[X_1] = \frac{n}{n}\mathbb{E}[X_1] = \mathbb{E}[X_1], $$

et pour la variance

$$ Var[M_n] = Var\Big[\frac{1}{n}\sum_{j=1}^n X_j \Big] = \frac{1}{n^2}Var\Big[\sum_{j=1}^n X_j\Big] = \frac{1}{n^2}\sum_{j=1}^n Var[X_j] = \frac{1}{n^2}\sum_{j=1}^n Var[X_1] = \frac{n}{n^2}Var[X_1] = \frac{1}{n} Var[X_1]. $$

Ainsi, plus $n$ est grand, plus la variance de $M_n$ devient petite $Var[M_n]$, c'est à dire qu'en moyenne $M_n$ s'éloigne de moins en moins de sa moyenne $\mathbb{E}[M_n] = \mathbb{E}[X]$.

4)

• Simuler $N=3000$ échantillons $(X_1,\dots, X_n)$ de $X$ de taille $n=5000$ (On fera une boucle for).
• Pour chaque échantillon évaluer

$$ U = \sqrt{\frac{n}{Var(X)}}(M_n - \mathbb{E}[X]). $$

• Tracer l'histogramme normalisée des différentes valeurs de $U$, et le comparer avec la densité de la loi $N(0,1)$.

• Qu'observe-t-on? Recommencer plusieurs fois, et essayer avec d'autre loi pour l'échantillon de $X$.

n = 5000
N = 3000

#X = st.norm(0,1)
X = st.expon(1)
#X = st.gamma(1)

Xn = X.rvs((N,n))

U = np.sqrt(n/X.var())*(np.sum(Xn, axis = 1) / n - X.mean())

x = np.linspace(-5, 5, 100)
plt.hist(U, bins=x, density=True, label='histogramme')
plt.plot(x, st.norm(0,1).pdf(x), lw=3, label="pdf N(0,1)")
plt.legend()
plt.show()

On voit que l'histogramme normé des $U$ correspond à la densité de la loi $N(0,1)$. Ceci nous dit que la loi de $U$ est proche de la loi $N(0,1)$ et ceci peu importe la loi de $X$ formant l'échantillon.