Examen de Python

Le rendu de ce devoir est sous forme d'un notebook avec les bouts de condes et les commentaires dans la partie texte (markdown).

À titre indicatif, la liste des fonctions à utiliser : np.log, scipy.special.gamma, scipy.special.gammaln, scipy.stats.norm.cdf, scipy.optimize.minimize, scipy.stats.t, np.var, np.std, scipy.stats.norm.pdf np.random.RandomState, ``

Estimation par maximum de vraisemblance du degré de liberté d'une loi de Student

Pour cet exercice, nous avons besoin d'appeler diverses fonctions. Il est indispensable consulter l'aide en ligne des différentes fonctions pour pouvoir les utiliser.

In [ ]:
import numpy as np
from scipy.special import gammaln
from scipy.optimize import minimize

rs = np.random.RandomState(19991231) ## pour fixer le générateur aléatoire
  • Simuler $1000$ pseudo-réalisations i.i.d d'une variable aléatoire $X$ de loi de Student de degré de liberté $\nu = 10$. Penser à faire appel à rs.
  • Sachant que $\text{Var}\big(X\big) = \frac{\nu}{\nu - 2}$, proposer une transformation permettant d'obtenir des réalisations d'une variable aléatoire standardisée (de variance $1$). Vérifier la transformation en affichant la variance empirique.
  • Proposer un estimateur de $\nu$ par la méthode des moments de $\nu$. Effectuer le calcul sur l'échantillon généré précédemment.

  • La densité d'une loi de Student est donnée par

    $$f(x;\nu,\mu,\sigma^{2})=\frac{\Gamma\left(\frac{\nu+1}{2}\right)}{\Gamma\left(\frac{\nu}{2}\right)}\,\frac{1}{\sqrt{\pi(\nu-2)}}\,\frac{1}{\sigma}\,\frac{1}{\left(1+\frac{\left(x-\mu\right)^{2}}{\sigma^{2}(\nu-2)}\right)^{\frac{\nu+1}{2}}}$$

    où $\mu = 0$, $\sigma = 1$ après standardisation et $\Gamma(\cdot)$ est la fonction gamma d'Euler. Écrire une fonction student_loglik qui prend en arguments un nombre positif (correspondant au degré de liberté) et un vecteur d'observations qui retourne la log-vraisemblance. Tester la fonction sur l'échantillon généré précédemment avec le vrai degré de liberté $\nu = 10$.

  • Proposer une légère modification de la fonction précédente qui permet d'utiliser la fonction minimize du module scipy.optimize pour calculer un estimateur de $\nu$ par maximum de vraisemblance. Utiliser l'estimateur par la méthode des moments comme point initial en limitant les bornes de minimisation à l'intervalle [5, 15]. Implémenter la procédure et afficher le résulat de l'estimation de $\nu$ par maximum de vraisemblance.

Modèle de régression Probit

Pour cet exercice, nous avons besoin de compléter les appels précédents par l'appel au sous-module 'stats' de la librairie 'scipy'

In [ ]:
from scipy import stats

Le modèle probit est un modèle de régression binaire (similaire à la régression logistique) qui vise à expliquer une variables aléatoire de Bernoulli $Y\sim B\big(p(x)\big)$ par l'intermédiaire de sa probabilité de succès. Dans le cas d'une seule variable explicative $x$, la fonction de lien associée à ce modèle est donnée par $$p(x)=\Phi\big(\beta_0 + \beta_1 x\big)$$ où $\Phi$ est la fonction de répartition d'une variable aléatoire normale centrée réduite $X \sim \mathcal N (0, 1)$ et $\big(\beta_0, \beta_1\big)$ sont les coefficients de régression.

  • Soit $U$ une variable aléatoire uniforme suivant une loi uniforme $U(0,1)$ et $p\in ]0,1[$. Vérifier que la loi de la variable aléatoire $Y = \mathbf{1}_{\{U \le p\}}$ où
$$Y = \begin{cases} 1 & \quad \text{si} \quad U \le p\\ 0 & \quad \text{sinon } \end{cases}$$

est une variables aléatoire suivant une loi de Bernoulli de paramètre $p$.

  • Simuler $1000$ couples $(x_i, y_i)$ suivant un modèle probit où $X_i \sim \mathcal N(0,1)$, $\beta_0 = 0$ et $\beta_1 = 1$. Pour simuler des réalisations d'une variable aléatoire normale centrée réduite ou suivant une loi uniforme, on peut faire appel à nouveau à l'objet rs.

  • Écrire une fonction probit_loglik qui prend en arguments le vecteur des deux coefficients, le vecteur des réalisations x et le vecteur des réalisations y qui retourne la log-vraisemblance des observations. Tester la fonction sur l'échantillon généré précédemment avec les vraies valeurs des coefficients de régresion $\beta_0 =0, \beta_1 = 1$. Penser à utiliser les fonctions listés précédemment.

  • Modifier la fonction précédente pour utiliser la fonction minimize du module scipy.optimize pour calculer un estimateur du couple $(\beta_0, \beta_1)$ par maximum de vraisemblance. Nous n'avons pas besoin de spécifier l'intervalle de recherche des paramètres. Tester plusieurs points de départ y compris les vraies valeurs du couple $(\beta_0, \beta_1)$.

  • Sachant que les dérivés de la log-vraisemblance d'une seule observation

$$ \frac{\partial \left\{y_i \ln\left(\Phi\left(\beta_0+\beta_1 x_i \right)\right) + \left(1-y_i\right) \ln\left(1-\Phi\left(\beta_0+\beta_1 x_i \right)\right)\right\}}{\partial \beta_j} $$

qui est

$$ y_i \frac{\phi\left(\beta_0+\beta_1 x_i \right)}{\Phi\left(\beta_0+\beta_1 x_i \right)} - \left(1-y_i\right)\frac{\phi\left(\beta_0+\beta_1 x_i \right)}{1-\Phi\left(\beta_0+\beta_1 x_i \right)} $$

pour $\beta_0$ et

$$ y_i x_i \frac{\phi\left(\beta_0+\beta_1 x_i \right)}{\Phi\left(\beta_0+\beta_1 x_i \right)} - \left(1-y_i\right)x_i\frac{\phi\left(\beta_0+\beta_1 x_i \right)}{1-\Phi\left(\beta_0+\beta_1 x_i \right)} $$

for $\beta_1$ où $\Phi(\cdot)$ est la fonction de répartition d'une variable aléatoire suivant une loi normale centrée réduite et $\phi(\cdot)$ sa densité de probabilité. Proposer une approximation de la matrice d'information de Fisher d'une observation. Rappeler la défintion de l'information de Fisher d'une observation. Proposer une approximation basée sur la loi des grands nombres.

Modules en python

Les modules permettent de combiner plusieurs fonctions dans un seul fichier Python et d'accéder à l'aide de import module puis de la syntaxe module.function. Supposons qu'un fichier nommé core.py (fichier à mettre dans le même répertoire du notebook) contienne le code suivant:

In [ ]:
r"""Demonstration module.
This is the module docstring.
"""


def square(x):
    r"""Returns the square of a scalar input
    """
    return x*x


def cube(x):
    r"""Returns the cube of a scalar input
    """ 
    return x*x*x

Les fonctions square et cube sont accessibles par d'autres fichiers dans le même répertoire comme suit

In [ ]:
import core
y = -3
print(core.square(y))
print(core.cube(y))

  • Écrire un module probit avec les fonctions internes probit_loglik (implémentée précédemment) et coef_estimate qui retourne les valeurs des coefficients de régression probit estimés par maximum de vraisemblance (penser à fixer le point initial de l'algorithme de minimisation).

  • (Bonus) Ajouter une fonction confint au module probit qui affiche un intervalle de confiance à $95\%$ de chaque paramètre de régression.