L’énoncé du devoir peut être téléchargé en cliquant sur le bouton PDF ci-dessus.
Le jeu de données du devoir peut être téléchargé en cliquant sur le bouton Dataset ci-dessus.
import numpy as np
import pandas as pd
import scipy.stats as stats
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(style="darkgrid")
from statsmodels.nonparametric.api import KernelReg
import statsmodels.nonparametric.api as nparam
import statsmodels.stats.stattools as stattools
import statsmodels.api as sm
df = pd.read_csv("DataReg.csv")
Exploration des propriétés de $g(x)$
Construction d’un estimateur non-paramétrique de la densité des $X_i$
Nous avons chargé les données dans la variable df. Nous réalisons un premier graphique pour visulaliser les $X_i$ et les $Y_i$ :
plt.scatter(df["X"], df["Y"], alpha=0.1);
Nous commençons par tracer une densité grâce à la méthode du noyau gaussien. La variable de lissage est de 0.002. Ci-dessous, le paramètre bw détermine la bandwidth :
sns.distplot(df["X"], kde=True,
kde_kws={"color": "k", "lw": 1, "label": "KDE", "bw": 0.002});
Ce niveau laisse apparaître beaucoup trop de variance dans l’estimation de notre densité. En effet, nous pouvons observer des pics et des creux autour de l’histogramme des $X_i$. Nous pouvons essayer avec un niveau de 0.4 :
sns.distplot(df["X"], kde=True,
kde_kws={"color": "k", "lw": 1, "label": "KDE", "bw": 0.4});
Ici, le problème inverse apparaît. Il y a un fort biais et la densité estimée épouse mal l’histogramme. Si nous laissons la fonction estimer bw, le résultat semble satisfaisant :
sns.distplot(df["X"], kde=True,
kde_kws={"color": "k", "lw": 1, "label": "KDE"});
Estimation de $h$ par validation croisée ou par une autre méthode.
La fonction sns.distplot estime une valeur satisfaisante de h, grâce à une formule :
$h=(\frac{4}{5})^\frac{1}{5}\sigman^\frac{-1}{5}$
La formule ci-dessus, donne, toute chose égale par ailleurs :
- un $h$ élevé si $\sigma$ est élevé. En effet si les points sont dispersés, la fenêtre de lissage devra être plus large,
- un $h$ élevé si le nombre de points est faible. Plus le nombre de points est élevé et plus la fenêtre peut être fine. Si la fenêtre est fine avec peu de points, l’estimateur pourrait avoir trop de variance.
sigma = np.std(df["X"])
mu = np.mean(df["X"])
n = len(df["X"])
h = (4/5)**(1/5)*sigma*n**(-1/5)
print("La valeur de h obtenue grâce à la formule :", round(h,6))
La valeur de h obtenue grâce à la formule : 0.027037
sns.distplot(df["X"], kde=True,
kde_kws={"color": "k", "lw": 1, "label": "KDE", "bw": h});
Nous trouvons une densité très similaire à celle obtenue sans forcer $h$. Cette méthode donne de bons résultats en pratique et est plus légère à implémenter qu’une validation croisée. Pour l’estimation de la densité en elle même, il est possible de se référer à ce document que j’ai réalisé en parallèle du cours.
Implémenter un QQ-plot pour vérifier empiriquement l’hypothèse $g(x) = 1$
sm.qqplot(df["X"], stats.uniform, fit=True,
line='45', markersize=1)
plt.show()
Ci-dessus, nous avons réalisés un QQ-plot avec en abscisses les quantiles théoriques d’une loi uniforme. Il est raisonnable de dire que la loi des $X_i$ ne suit pas une loi uniforme. Nous pouvons vérifier en générant 1000 points à partir d’une loi uniforme et en refaisant le même graphique :
sm.qqplot(stats.uniform.rvs(size=1000), stats.uniform, fit=True,
line='45', markersize=1)
plt.show()
La densité n’étant pas uniforme, estimer un même paramètre de lissage pour l’ensemble de l’espace des $X_i$ (dans le cadre de la régression) n’est pas conseillé. Il serait mieux d’estimer une régression avec des polynômes locaux (LOWESS).
En effet le paramètre $h$ est sensible au nombre de points observés au voisinage du point à estimer. Si le nombre de points varie significativement d’une fenêtre à l’autre, alors la qualité de l’estimation peut également varier.
Reconstruction de $r(x)$
Est-il plausible de penser que la fonction $r$ est linéaire ?
Afin de répondre à cette question nous utiliserons la fonction KernelReg. La documentation, disponible à cette adresse, nous informe que l’implémentation de Nadaraya-Watson est disponible via la paramètre reg_type='lc'. Le comportement par défaut de la fonction est de faire un ajustement via polynômes locaux de degré supérieur à 1 accessible via reg_type='ll'.
Pour choisir le paramètre de lissage, nous demanderons à la fonction de réaliser une validation croisée sur les moindres carrés via bw='cv_ls'. Cet appel prend environ 35 secondes sur mon système :
%%time
model = KernelReg(endog=df["Y"].values, exog=df["X"].values,
reg_type='lc', var_type='c', bw='cv_ls',
defaults=nparam.EstimatorSettings(efficient=True))
/opt/conda/lib/python3.6/site-packages/statsmodels/nonparametric/kernel_regression.py:215: RuntimeWarning: invalid value encountered in true_divide
G = G_numer / G_denom
/opt/conda/lib/python3.6/site-packages/statsmodels/nonparametric/kernel_regression.py:227: RuntimeWarning: invalid value encountered in true_divide
B_x = d_mx / f_x - G * d_fx / f_x
/opt/conda/lib/python3.6/site-packages/statsmodels/nonparametric/kernel_regression.py:228: RuntimeWarning: invalid value encountered in true_divide
B_x = (G_numer * d_fx - G_denom * d_mx) / (G_denom**2)
CPU times: user 35.3 s, sys: 72 ms, total: 35.3 s
Wall time: 35.5 s
La méthode .fit() permet d’estimer les $\hat{Y}$ à partir des $X_i$ :
plt.scatter(df["X"], df["Y"], alpha=0.1)
plt.scatter(df["X"].values, model.fit(df["X"].values)[0],
color="r", alpha=.5, s=.05);
En observant ce graphique, nous voyons que $r$ n’est pas linéaire. Si on voulait modéliser le comportement des points à l’aide d’un modèle linéaire, nous pourrions diviser les $X_i$ en au moins deux sous-espaces. Le premier irait de 0 à environ 0.5 et le second irait de 0.5 à 0.7. L’espace allant de 0.7 à 1 semble plus difficile à modéliser car il y a moins de points dans cette zone. En ce sens, la régression non-paramétrique est une première étape permettant de repérer un changement de tendance.
Construire un estimateur non-paramétrique $\hat{r}_{n,h}(x)$ de $r(x)$
Tout comme le problème d’estimation de la densité, nous faisons face à un compromis biais-variance. Un $h$ trop petit va créer trop de variance dans l’estimation de notre régression. Ci-dessous, nous observons que les points estimés varient trop; la courbe rouge capture du bruit :
model_var = KernelReg(endog=df["Y"].values, exog=df["X"].values,
reg_type='lc', var_type='c', bw=np.array([0.002]),
defaults=nparam.EstimatorSettings(efficient=True))
plt.scatter(df["X"], df["Y"], alpha=0.1)
plt.scatter(df["X"].values, model_var.fit(df["X"].values)[0],
color="r", alpha=.5, s=.05);
A l’inverse, un paramètre de lissage trop grand aura tendance à biaiser le modèle et à empêcher la régression de capturer la tendance :
model_bias = KernelReg(endog=df["Y"].values, exog=df["X"].values,
reg_type='lc', var_type='c', bw=np.array([0.2]),
defaults=nparam.EstimatorSettings(efficient=True))
plt.scatter(df["X"], df["Y"], alpha=0.1)
plt.scatter(df["X"].values, model_bias.fit(df["X"].values)[0],
color="r", alpha=.5, s=.05);
print("h trop petit (variance):", model_var.bw[0])
print("h trop grand (biais):", model_bias.bw[0])
print("h optimal :", model.bw[0])
h trop petit (variance): 0.002
h trop grand (biais): 0.2
h optimal : 0.0179801912776
Nous conserverons donc la fenêtre choisie par l’implémentation de la validation croisée prévue par le package.
Discuter du choix automatique de la fenêtre
Un cas où le $h$ est le même pour tout $x$
Ce cas est cohérent si les $X_i$ sont distribués de manière uniforme dans l’espace. La variable de lissage optimale dépend du nombre d’observations sur une fenêtre donnée. Si le nombre de points varie, alors il est probable que l’estimation de la largeur de la fenêtre soit une moyenne entre les zones avec beacoup de points observés et les zones avec peu de points observés. Le risque est que la fenêtre estimée de manière automatique (grâce à une validation croisée par exemple) minimise l’erreur globale tout en n’étant pas optimale sur les sous-espaces (localement).
Nous pouvons simuler un exemple pour nous rendre compte que Nadaraya-Watson n’arrive pas à reconstituer une régression simple :
x_simul = np.linspace(start=0.1, stop=1, num=500)
y_simul = x_simul + np.random.randn(500)/50
mask = (x_simul<0.4) | (x_simul>0.6)
model_example = KernelReg(endog=y_simul[mask], exog=x_simul[mask],
reg_type='lc', var_type='c', bw='cv_ls',
defaults=nparam.EstimatorSettings(efficient=True))
plt.scatter(x_simul[mask], y_simul[mask], alpha=0.1)
plt.scatter(x_simul, model_example.fit(x_simul)[0],
color="r", alpha=.5, s=.05);
/opt/conda/lib/python3.6/site-packages/statsmodels/nonparametric/kernel_regression.py:215: RuntimeWarning: invalid value encountered in true_divide
G = G_numer / G_denom
/opt/conda/lib/python3.6/site-packages/statsmodels/nonparametric/kernel_regression.py:227: RuntimeWarning: invalid value encountered in true_divide
B_x = d_mx / f_x - G * d_fx / f_x
/opt/conda/lib/python3.6/site-packages/statsmodels/nonparametric/kernel_regression.py:228: RuntimeWarning: invalid value encountered in true_divide
B_x = (G_numer * d_fx - G_denom * d_mx) / (G_denom**2)
le cas ou le choix du $h$ dépend localement de la fonction à estimer
Le lissage par polynômes locaux est une extension naturelle du lissage par moyenne locale de Nadaraya-Watson. Une régression par polynômes locaux implique de modéliser la réponse avec un polynôme en utilisant une méthode des moindres carrés ordinaires pondérée localement. Les polynômes de degrés supérieur à un ont de meilleurs propriétés (moins de biais) que le polynôme d’ordre zéro utilisé dans le cadre de Nadaraya-Watson.
Pour nous en convaincre, nous pouvons refaire tourner le même exemple en utilisant l’option 'll' à la place de 'lc'. Nous voyons que le fit est meilleur dans l’espace sans points ainsi que sur les côtés (aux alentours de 0 et de 1) :
x_simul = np.linspace(start=0.1, stop=1, num=500)
y_simul = x_simul + np.random.randn(500)/50
mask = (x_simul<0.4) | (x_simul>0.6)
model_example = KernelReg(endog=y_simul[mask], exog=x_simul[mask],
reg_type='ll', var_type='c', bw='cv_ls',
defaults=nparam.EstimatorSettings(efficient=True))
plt.scatter(x_simul[mask], y_simul[mask], alpha=0.1)
plt.scatter(x_simul, model_example.fit(x_simul)[0],
color="r", alpha=.5, s=.05);
Etude de la loi des $\xi_i$
Implémenter ces estimateurs sur le jeu de données
$U_n =\frac{1}{2(n-1)}\sum_{i=1}^{n-1}(Y_{i+1} - Y_i)^2$
df["Y+1"] = df["Y"].shift()
df.head(5)
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
un = ((df["Y+1"] - df["Y"])**2).sum()/(2*9999)
print("La valeur de Un :", un)
La valeur de Un : 3.91335271843
$$V_n =\frac{1}{2(n-1)}\sum_{i=1}^{n-1}(\tilde{Y}_{i+1} - \tilde{Y}_i)^2$$
df_bis = df.loc[:,["X", "Y"]].copy()
df_bis.sort_values(by="X", inplace=True)
df_bis["Y+1"] = df_bis["Y"].shift()
vn = ((df_bis["Y+1"] - df_bis["Y"])**2).sum()/(2*9999)
print("La valeur de Vn :", vn)
La valeur de Vn : 1.00567446149
Que représentent ces deux quantités ?
Nous trouvons qu’empiriquement $U_n$ est la variance des $Y_i$ :
np.var(df["Y"])
3.918446241816131
Estimation de $x\mapsto{\mu(x)}$
Pour répondre à cette question, nous commençons par construire un estimateur pour les $X_i$ dont $i$ est compris entre 1 et 5000. Nous forcerons $h$ à 0.0179801912776 :
model_1_500 = KernelReg(endog=df["Y"][:5000].values, exog=df["X"][:5000].values,
reg_type='lc', var_type='c', bw=np.array([0.0179801912776]),
defaults=nparam.EstimatorSettings(efficient=True))
plt.scatter(df["X"], df["Y"], alpha=0.1)
plt.scatter(df["X"].values, model_1_500.fit(df["X"].values)[0],
color="r", alpha=.5, s=.05);
La seconde étape est d’estimer les résidus pour les $i$ compris entre 501 et 10000 :
df["Y_tilde"] = df["Y"] - model_1_500.fit(df["X"].values)[0]
df.loc[:5000, "Y_tilde"] = np.nan
df[4990:5010]
.dataframe tbody tr th {
vertical-align: top;
}
.dataframe thead th {
text-align: right;
}
En déduire un estimateur de $x\mapsto{\mu(x)}$ et l’implémenter graphiquement
Nous pouvons tracer la densité de $\tilde{Y}_i$ pour les $i$ allant de 501 à 10000 :
sns.distplot(df["Y_tilde"][5001:], kde=True,
kde_kws={"color": "k", "lw": 1, "label": "KDE"});
Grâce à notre estimateur de la densité, nous observons que la distribution est symétrique et que la forme générale ressemble à une distribution gaussienne.
Approximativement, les données semblent suivre une distribution gaussienne.
Nous pouvons également tracer un qqplot pour s’assurer que les quantiles empiriques correspondent aux quantiles théoriques :
sm.qqplot(df["Y_tilde"][5001:], stats.norm, fit=True,
line='45', markersize=1)
plt.show()
Protocole numérique pour vérifier la normalité
Nous pouvons implémenter un test de Jarque-Bera afin de vérifier si les résidus sont normaux. L’implémentation du test renvoie la valeur de la statistique de test, la p-valeur associée à la statistique de test, le skewness et le kurtosis :
stattools.jarque_bera(df["Y_tilde"][5001:])
(1.1921682131185092,
0.5509649372774027,
0.023340207515106737,
3.05953562787107)
D’après la documentation : “The Jarque-Bera test statistic tests the null that the data is normally distributed against an alternative that the data follow some other distribution.”
Nous ne pouvons donc pas rejeter l’hypothèse nulle au regard de la p-valeur de 55%. De plus le skewness et le kurtosis correspondent à une loi normale.
Il est donc probable que les résidus suivent une loi normale.
Pour vérifier, nous pouvons faire le test en injectant des données générées à l’aide d’une loi normale :
stattools.jarque_bera(stats.norm.rvs(size=100000))
(2.1677044529506437,
0.33828984327689721,
0.0015437569229027144,
2.977400966238352)
