TP 6 - Probabilités, loi des grands nombres.

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as st

I) Estimation d'une densité de loi de probabilité.

On considère un échantilllon $(X_1,\dots,X_n)$ de la v.a $X$, et la fonction

$$ \hat f_{n,h} (x) =\frac{1}{n h} \sum_{j=1}^n K\Big(\frac{x-X_j}{h}\Big) $$

où $h>0$ est appelé le paramètre de lissage et $K$ un noyau. On considérera les noyaux

$$ K_1(x) = \mathbf{1}_{(|x|\leq 1/2)}=\begin{cases} 1&\text{if } |x|\leq 1/2 \\ 0 & \text{sinon} \end{cases}\qquad \text{and} \qquad K_2(x) = \frac{e^{-x^2/2}}{\sqrt{2\pi}}. $$

Remarque: En appliquant la loi des grands nombres on a

$$ \hat f_{n,h} (x) \simeq \frac{1}{h}\mathbb{E}\Big[K\Big(\frac{x-X}{h}\Big)\Big] = \int_{-\infty}^{+\infty} \frac{1}{h} K\Big(\frac{x-t}{h}\Big)f_{X}(t)dt, $$

où $f_{X}$ est la densité de $X$. Or pour $h$ petit on a $K(x/h)/h \simeq \delta(x)$ (la masse de Dirac en 0), et donc dans ce cas

$$ \hat f_{n,h} (x) \simeq \int_{-\infty}^{+\infty} \delta(x-t)f_{X}(t)dt = f_{X}(x). $$

1)

• Implémenter les noyaux $K_1$ et $K_2$, ainsi que la fonction $\hat f_{n,h}$ qui prend en entrée un échantillon.
• Génerer un échantillon de taille $n=5000$ de la loi $N(0,1)$.
• Illustrer l'histogramme normalisée de cette échantillon, et le comparer avec le graphe de la fonction $\hat f_{n,h}$. On jouera avec le paramètre $h$ pour voir qu'elle role il joue.
• Essayer avec d'autres lois.

K1 = lambda x: 1*(np.abs(x)<0.5)

K2 = lambda x: np.exp(-0.5 * x**2)/np.sqrt(2*np.pi)

def hat_f(x, h, X, K=K2):
    x_tmp = x.reshape((-1,1))
    X_tmp = X.reshape((1,-1))
    
    return np.mean(K((x_tmp - X_tmp)/h), axis=1) / h

X = st.norm(0,1)

n=5000
Xn = X.rvs(n)

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

y = hat_f(x, h, Xn, K2) 

plt.hist(Xn, bins=x, density=True, label='histogramme')
plt.plot(x, y, lw=3, label='estimation')

plt.legend()
plt.show()

Comme menstionné plus haut, $h$ est un paramètre de lissage. Plus $h$ est grand, plus l'estimation de la densité de $X$ est lisse.

Prenons maintenant un échantillon la concatenation d'un échantillon de loi $N(0,1)$ et de la loi $N(3,0.5^2)$.

Y = st.norm(3, 0.5)
Yn = Y.rvs(n)
Zn = np.concatenate((Xn, Yn))

x = np.linspace(-4, 6, 100)
h = 0.2
y = hat_f(x, h, Zn)

plt.hist(Zn, bins=x, density=True, label='histogramme')
plt.plot(x, y, lw=3, label='estimation')

plt.legend()
plt.show()

2)

Chercher à obtenir les mêmes résultats que dans la question précédente avec la fonction KernelDensity du module sklearn.neighbors de la librairie scikitlearn.

from sklearn.neighbors import KernelDensity

kde = KernelDensity(kernel='gaussian', bandwidth=0.2).fit(Zn[:,None])

log_dens = kde.score_samples(x[:, None])

plt.hist(Zn, bins=x, density=True, label='histogramme')
plt.plot(x, np.exp(log_dens), lw=3, label='scikit-learn KDE')

plt.legend()
plt.show()

II) Moyenne locale

On va considérer le signal suivant

$$ f(t) = e^{-0.2 \pi t^2}\cos(2\pi a t)\qquad a = 1. $$

sur l'intervalle de temps $[-4,4]$.

1)

• Illustrer ce signal sur une grille de temps de 200 points.
• Illustrer ce même signal sur lequel auquel on rajoute du bruit en chaque point à l'aide de la loi $N(0,0.1^2)$.

a = 1.
f = lambda t: np.exp(- 0.2 * np.pi * t**2)*np.cos(2*np.pi*a*t)

t = np.linspace(-4, 4, 200)

noise = st.norm(0,0.05).rvs(len(t))

f_noise = f(t) + noise

plt.plot(t, f_noise, "o", label='signal bruité')
plt.plot(t, f(t), label="vrai signal")

plt.legend()
plt.show()

2)

• Implementer une fonction qui calcule en chaque point $t_n$ de la grille de temps la moyenne locale des valeurs du signal bruité, c'est à dire la moyenne des valeurs du signal aux temps $t_{n-h},\dots, t_n, t_{n+h}$. Ici $h$ représente le rayon de la fenêtre pour calculer la moyenne. On pourra étendre le signal bruité par $0$.
• Jouer avec la taille de la grille de temps et le nombre de point qui servent à calculer la moyenne.
• Comment expliquez vous ce que vous observez?

def moyenne_locale(x, h):
    
    x_smoothed = np.zeros(len(x))
    x_tmp = np.concatenate((np.zeros(h), x, np.zeros(h)))
    
    for j in range(h, len(x)+h):
        x_smoothed[j-h] = np.mean(x_tmp[j-h:j+h])
        
    return x_smoothed

f_smoothed = moyenne_locale(f_noise, 4)

plt.plot(t, f_noise, "o", label='signal bruité', markersize=1.2)
plt.plot(t, f(t), lw=3, label="vrai signal")
plt.plot(t, f_smoothed, lw=3, label='signal moyenne locale')

plt.legend()
plt.show()

• Plus la fenêtre pour calculer la moyenne est large, plus le signal est lisse (moins bruité).

En effet,

$$ Var\Big( \frac{ f_{noise}(t_{n-h}) +\dots + f_{noise}(t_{n}) +\dots+ f_{noise}(t_{n+h})}{2h} \Big) = \frac{1}{(2h)^2}\big(Var( f_{noise}(t_{n-h}) ) +\dots + Var(f_{noise}(t_{n})) +\dots+ Var(f_{noise}(t_{n+h}))\big).$$

Mais pour chaque $Var(f_{noise}(t_{j}))$, comme on a

$$ f_{noise}(t_{j}) = f(t_j)+\varepsilon_j $$

ces variances sont données par

\begin{align*} Var(f_{noise}(t_{j})) & = \mathbb{E}[(f_{noise}(t_j)-\mathbb{E}[f_{noise}(t_j)])^2] \\ & = \mathbb{E}[(f(t_j)+\varepsilon_j - \mathbb{E}[ f(t_j)+\varepsilon_j])^2] \\ & = \mathbb{E}[(f(t_j)+\varepsilon_j - f(t_j)-\mathbb{E}[ \varepsilon_j])^2] \\ & = \mathbb{E}[( \varepsilon_j - \mathbb{E}[ \varepsilon_j])^2] \\ & = Var(\varepsilon_j)\\ & = \sigma^2. \end{align*}

Ainsi,

$$ Var\Big( \frac{ f_{noise}(t_{n-h}) +\dots + f_{noise}(t_{n}) +\dots+ f_{noise}(t_{n+h})}{2h} \Big) =\frac{\sigma^2}{2h} \ll 1 $$

pour un rayon de fenêtre $h$ suffisament grand. Dans ce cas, la variance est petite et la moyenne locale reste proche de son espérance, ce qui donne un effet plus lisse. Par contre l'espérance de cette moyenne locale est donnée par

\begin{align*} \mathbb{E}\Big[ \frac{ f_{noise}(t_{n-h}) +\dots + f_{noise}(t_{n}) +\dots+ f_{noise}(t_{n+h})}{2h} \Big] & = \frac{ f(t_{n-h}) +\dots + f(t_{n}) +\dots+ f(t_{n+h})}{2h} \\ & +\underbrace{\mathbb{E}\Big[ \frac{ \varepsilon_{n-h} +\dots + \varepsilon_{n} +\dots+ \varepsilon_{n+h}}{2h} \Big]}_{=0} \\ &= \frac{ f(t_{n-h}) +\dots + f(t_{n}) +\dots+ f(t_{n+h})}{2h}. \end{align*} Cependant, cette dernière moyenne peu ne pas être proche de $f(t_n)$ pour une grille de temps grossière, ou un rayon $h$ trop grand.

• Plus la grille de temps est fine et plus le signal en lui même est mieux reconstruit.

En effet car si la grille est fine, les différentes valeurs $f(t_{n-h}),\dots, f(t_n), \dots, f(t_{n+h})$ sont prochez, et plus la moyenne locale autour du point $t_n$ reste proche de $f(t_n)$:

$$ \frac{f(t_{n-h})+\dots + f(t_n) + \dots + f(t_{n+h})}{2h} \simeq f(t_n).$$