Méthode de Newton-Raphson pour la maximisation de la vraisemblance pour la régression logistique¶
Python works with packages. A large number of packages can be dowloaded all together in the Anaconda environment. Start downloading here (version >3) https://www.anaconda.com/download/ Anaconda comes with
- Spyder, a GUI for python
- Jupyter Notebook, how to create great documents
Enjoy!
Bee careful :
- Python is sensible to spaces in the beginning of lines,
- Python is sensible to number type : $3$ and $3.$ are totally different: integer and float types.
import warnings
warnings.filterwarnings('ignore')
import numpy as np
import numpy.random as npr
import matplotlib.pyplot as plt
import pandas as pd
Les données¶
On a un dataset formé par deux variables qualitatives binaires $X$ et $Y$. Telles que les valeurs prises par :
- $X$ : fumeuse ou non fumeuse
- $Y$ : FPN ou poids normal
Soit donc les valeurs définissant les $n$ patients $\forall i=1..n,(x_i,y_i)$.
En fait on le génère :
X = [0]*115+[1]*74
Y = [1]*29+[0]*130+[1]*30
data = np.array([X,Y]).T
pd_data = pd.DataFrame(data=data,columns=["Tabagisme (X)","Poids (Y)"])
poids_tabagisme = pd.crosstab(index=pd_data["Tabagisme (X)"],
columns=pd_data["Poids (Y)"])
poids_tabagisme
La méthode de Newton-Raphson¶
C'est une méthode qui permet de trouver une valeur, proche de la valeur initiale, qui annule la fonction. Elle nécessite que la fonction soit dérivable sur tous les points de l'intervalle considéré et que cette dérivée ne s'annule jamais.Intuition en uni-dimensionnel¶
L'idée est de trouver à chaque itération un point qui annule l'approximation linéaire locale de la fonction au point courant.
Equation de mise à jour de la suite :
- $x_0$ un point initial
- Pour la recherce d'un :
- minimum : $x_{i+1}=x_i-\frac{f(x_i)}{f'(x_i)}$
- maximum : $x_{i+1}=x_i+\frac{f(x_i)}{f'(x_i)}$
Le script suivant permet de représenter la méthode de construction de la suite. Sur ce graphique, $F_{x_i}(x)$ est la fonction tangeante à la courbe en $x_i$
import numpy as np
import matplotlib.pyplot as plt
def f_(t):
return(.1*t**3 + 1.5*t + 1)
u = np.linspace(-1, 4)
t_0 = 3
fig, ax = plt.subplots()
ax.plot(u, f_(t_0)+(u-t_0)*(0.3*t_0**2+1.5), 'b:')
ax.plot(u, f_(u))
ax.plot(t_0, f_(t_0), marker='o')
ax.plot(t_0-f_(t_0)/(0.3*t_0**2+1.5),0, marker='o')
ax.axhline(color='black')
ax.axvline(color='black')
x=[t_0,t_0-f_(t_0)/(0.3*t_0**2+1.5),-0.2,-0.75]
y=[f_(t_0),0,0.6,-8]
labels=["$x_i$","$x_{i+1}$","$f(x)$","$F_{x_i}(x)$"]
for i, txt in enumerate(labels):
ax.annotate(txt, (x[i], y[i]),ha='right', va='bottom')
plt.show()
QUESTION Démontrer que l'équation de mise à jour de la suite trouve bien à chaque itération le point qui annule l'approximation linéaire locale de la fonction.
Si l'approximation linéaire locale de la courbe, notée $F_{x_i}:x\rightarrow F_{x_i}(x)$ qui vérifie donc $f(x_i)=F_{x_i}(x_i)$ et $f'(x_i)=F_{x_i}'(x_i)$, s'annule en $x_{i+1}$ alors on peut écrire l'équation de la pente en $x_i$ telle que : $$ f'(x_i)=\frac{F_{x_i}(x_i)-F_{x_i}(x_{i+1})}{x_i-x_{i+1}}=\frac{f(x_i)}{x_i-x_{i+1}}, $$ en réordonnant les termes il vient $$ x_{i+1}=x_i-\frac{f(x_i)}{f'(x_i)}, $$ qui a bien la forme recherchée.
Application au calcul de la vraisemblance¶
Soit un jeu de données (dataset) formé par deux variables X et Y telles que Y doit être expliqué par X grâce au modèle logistique suivant :
$$ \pi_i=\pi(x_i)=\frac{e^{\beta_0+\beta_1x_i}}{1+e^{\beta_0+\beta_1x_i}}, $$
on peut alors écrire
$$ (Y|X=x)\sim \mathcal{B}(\pi(x))\\ \mathbb{P}(Y=1|X=x)=\pi(x)\quad\text{et}\quad\mathbb{P}(Y=0|X=x)=1-\pi(x) $$
Le but du jeu est de qualibrer $\beta_0$ et $\beta_1$ afin que la vraissemblance du modèle soit la plus forte sachant les observations déjà récoltées. En supposant les échantillons indépendants, il vient que la probabilité d'observer la variable Y sachant la variable X s'écrit :
$$ \mathbb{P}(Y_1=y_1,\cdots,Y_n=y_n|X_1=x_1,\cdots,X_n=x_n) = \prod_{i=1}^n\mathbb{P}(Y_i=y_i|X_i=x_i), $$
que nous appelons vraisemblance et qui s'écrit ici, $\forall i\in\{1,\cdots,n\}\quad \pi_i^{y_i}(1-\pi_i)^{1-y_i}$ et donc en global
$$ \mathcal{L}(\beta_0,\beta_1|X,Y)=\prod_{i=1}^n\pi_i^{y_i}(1-\pi_i)^{1-y_i}. $$
On passe par la méthode du maximum de vraissemblance pour écrire la condition d'optimalité. On apelle alors $l(\beta_0,\beta_1|X,Y)$ la log-vraisemblance du problème.
On définit la matrice de Score associée à la fonction de $2$ variables $l(\beta_0,\beta_1|X,Y)$,
$$ \mathcal{U}(\beta_0,\beta_1|X,Y)= \begin{pmatrix} \frac{\partial l}{\partial \beta_0}(\beta_0,\beta_1|X,Y)\\ \frac{\partial l}{\partial \beta_1}(\beta_0,\beta_1|X,Y) \end{pmatrix}, $$ estimée au point de minimum $(\beta_0^\star,\beta_1^\star)$ on doit avoir $$ \mathcal{U}(\beta_0^\star,\beta_1^\star|X,Y)= \begin{pmatrix} 0\\ 0 \end{pmatrix} =0 $$
En pratique il est très compliqué de résoudre ce problème on passe donc par une méthode itérative d'approximation : La méthode de Newton-Raphson matricielle.
Si l'on doit résoudre le problème $\mathcal{U}(\beta_0^\star,\beta_1^\star|X,Y)=0$ et que l'on peut calculer la matrice Hessienne de la fonction $l$, qui s'écrit
$$ \mathcal{H}(\beta_0,\beta_1|X,Y)= \begin{pmatrix} \frac{\partial^2 l}{\partial \beta_0^2} & \frac{\partial^2 l}{\partial \beta_0\partial \beta_1}\\ \frac{\partial^2 l}{\partial \beta_0\partial \beta_1} & \frac{\partial^2 l}{\partial \beta_1^2} \end{pmatrix}, $$
alors la mise à jour par la méthode de Newton-Raphson s'écrit, en notant $\beta^{(i)}=\begin{pmatrix} \beta_0^{(i)}\\ \beta_1^{(i)} \end{pmatrix}$ , dans le cas à plusieurs variables :
$$\beta^{(i+1)}=\beta^{(i)}+\mathcal{H}(\beta_0^{(i)},\beta_1^{(i)}|X,Y)^{-1}\mathcal{U}(\beta_0^{(i)},\beta_1^{(i)}|X,Y).$$
ATTENTION : On recherche ici le maximum de la vraisemblance donc la descente est une montée, d'où le + et non pas -.
Il convient de trouver une initialisation raisonnable. Par exemple, on sait que $\mathbb{P}(Y=1|X=0)=\pi(0)=\frac{e^{\beta_0}}{1+e^{\beta_0}}$ et donc par inversion il vient $\beta_0=log(\frac{\pi_0}{1-\pi_0})$. par la règle de Bayes il vient :
$$ \mathbb{P}(Y=1|X=0)=\frac{\mathbb{P}(Y=1\cap X=0)}{\mathbb{P}(X=0)}, $$ où l'on sait que $\mathbb{P}(Y=1\cap X=0)=\frac{29}{30+29+44+86}$, $\mathbb{P}(X=0)=\frac{86+29}{30+29+44+86}$ et donc pour valeur initiale on a
$$ \beta_0^{(0)}=log(\frac{29}{86})\approx -1.087051 $$
Les différentes fonctions s'expriment tel que :
- La log-vraissemblance $$ l(\beta_0,\beta_1|X,Y)=\sum_{i=1}^ny_i(\beta_0+\beta_1x_i)-log(1+e^{\beta_0+\beta_1x_i}). $$
- Son Score $$ \mathcal{U}(\beta_0,\beta_1|X,Y)= \begin{pmatrix} \frac{\partial l}{\partial \beta_0}\\ \frac{\partial l}{\partial \beta_1} \end{pmatrix} $$
$$ \mathcal{U}(\beta_0,\beta_1|X,Y)= \begin{pmatrix} \sum_{i=1}^ny_i-\frac{e^{\beta_0+\beta_1x_i}}{1+e^{\beta_0+\beta_1x_i}}\\ \sum_{i=1}^nx_i(y_i-\frac{e^{\beta_0+\beta_1x_i}}{1+e^{\beta_0+\beta_1x_i}}) \end{pmatrix}. $$
- Et la Hessienne associée $$ \mathcal{H}(\beta_0,\beta_1|X,Y)= \begin{pmatrix} \frac{\partial^2 l}{\partial \beta_0^2} & \frac{\partial^2 l}{\partial \beta_0\partial \beta_1}\\ \frac{\partial^2 l}{\partial \beta_0\partial \beta_1} & \frac{\partial^2 l}{\partial \beta_1^2} \end{pmatrix} $$
$$ \mathcal{H}(\beta_0,\beta_1|X,Y)= \begin{pmatrix} \sum_{i=1}^n\frac{e^{\beta_0+\beta_1x_i}}{(1+e^{\beta_0+\beta_1x_i})^2}& \sum_{i=1}^nx_i\frac{e^{\beta_0+\beta_1x_i}}{(1+e^{\beta_0+\beta_1x_i})^2} \\ \sum_{i=1}^nx_i\frac{e^{\beta_0+\beta_1x_i}}{(1+e^{\beta_0+\beta_1x_i})^2}& \sum_{i=1}^nx_i^2\frac{e^{\beta_0+\beta_1x_i}}{(1+e^{\beta_0+\beta_1x_i})^2} \end{pmatrix}. $$
Application : programmation sous Python¶
On définit plusieurs fonctions permettant de calculer les valeurs de la fonction aux points courants.
La fonction de vraissemblance :
def log_like(data,th):
b0 = th[0]
b1 = th[1]
X_in = data[:,0:1]
Y_real = data[:,1:2]
Y_est = b0+b1*X_in
out = np.dot(Y_real.T,Y_est)-np.sum(np.log(1+np.exp(Y_est)))[0,0]
return(out)
def U_score(data,th):
b0 = th[0]
b1 = th[1]
X_in = data[:,0:1]
Y_real = data[:,1:2]
exp_Y_est = np.exp(b0+b1*X_in)
p = exp_Y_est/(1+exp_Y_est)
delta = Y_real-p
L_prime_0 = np.sum(delta)
L_prime_1 = np.dot(X_in.T,delta)[0,0]
return(np.array([[L_prime_0],[L_prime_1]]))
def Hess(data,th):
b0 = th[0]
b1 = th[1]
X_in = data[:,0:1]
exp_Y_est = np.exp(b0+b1*X_in)
p_2 = exp_Y_est/(1+exp_Y_est)**2
H_0_0 = np.sum(p_2)
H_0_1 = np.dot(X_in.T,p_2)
H_1_1 = np.dot(X_in.T**2,p_2)
return(np.array([[H_0_0,H_0_1],[H_0_1,H_1_1]]))
Une fonction qui fait réellement le calcul itératif
def newton_raphson(data,errorMin=1e-9,maxIter=100):
pd_data = pd.DataFrame(data=data,columns=["X","Y"])
crosstab = pd.crosstab(index=pd_data["X"],
columns=pd_data["Y"]).values
b_0_init = np.log(crosstab[0,1]/(crosstab[1,1]))
theta = np.array([[b_0_init],[0]])
err = 10*errorMin
iteration = 1
crit = True
while crit:
U = U_score(data,theta)
H = Hess(data,theta)
theta = theta + np.dot(np.linalg.inv(H),U)
err = np.linalg.norm(U)
crit_1 = err>errorMin
crit_2 = iteration<maxIter
crit = crit_1 & crit_2
iteration = iteration + 1
if crit_2==False:
print("Maximum number of iterations reached")
return(theta)
Que l'on peut appliquer à notre jeu de données
theta = newton_raphson(data)
print(theta)
Le signe de $\beta_1$ renseigne sur le lien entre X et Y : la probabilité de $Y=1$ augmente avec celle de $X=1$ lorsque $\beta_1>0$ et diminue dans l'autre cas.
Afin d'estimer la variance on utilise la matrice d'information de Fisher
Info_Fisher = np.linalg.inv(Hess(data,[theta[0,0],theta[1,0]]))
std_b0 = np.sqrt(Info_Fisher[0,0])
std_b1 = np.sqrt(Info_Fisher[1,1])
print(std_b0)
print(std_b1)