In [1]:
import numpy as np
In [2]:
d = 7 # dimension des observations
N = 100 # nombre d'observations

La loi a posteriori des paramètres $\beta\in\mathbb R^{d+1}$ sachant les données $X=(X^{(1)},\ldots,X^{(n)})\in\mathbb (R^d)^n$ et $y=(y_1,\ldots,y_n)\in\{0,1\}^n$ d'un modèle de régression logisitique Logit avec un prior gaussien $\mathcal N(0,I_{d+1})$ est donnée par $$ p(\beta|X,y)\propto\prod_{i=1}^n\sigma(X_i, \beta)^{y_i}(1-\sigma(X_i, \beta))^{1-y_i}\; e^{-\|\beta\|^2/2} $$ où $\propto$ veut dire à une constante multiplicative près et avec $$ \sigma(x,\beta):=\frac1{1+e^{-(\beta_0+\beta_1 x_1+ \ldots +\beta_d x_d)}},\qquad x=(x_1,\ldots,x_d). $$ On veut simuler de la loi a posteriori $p(\beta|X,y)$ et, par exemple, obtenir une approximation de $$ \hat \theta := \int_{\mathbb R^d} \theta p(\beta|X,y) d\theta, $$ qui est l'estimateur qui minimise l'erreur quadratique $$\eta\mapsto \int_{\mathbb R^d} \ \|\eta-\theta\|^2 p(\beta|X,y) d\theta $$

In [3]:
def sigma(x,b):
    y = np.concatenate((np.ones(1),x)) # on rajoute 1 au vecteur x
    return 1/(1+np.exp(-np.dot(y,b)))


def UNposterior(b,X,Y): # "densité" a posteriori sans renormalisation
    # X = matrice N lignes d colonnes, Y = N-vecteur de 0-1, b = d-vecteur
    N = len(X[:,0]) # nombre d'observations
    P = 1
    for i in range(N):   # utiliser np.prod plutôt ? 
        if Y[i]==1:
            P *= sigma(X[i],b)
        else:
            P *= 1-sigma(X[i],b) 
    return P*np.exp(-np.dot(b,b)/2)
In [4]:
X = np.random.rand(N,d)
Y = np.random.randint(2,size=N)
b = np.random.multivariate_normal(np.zeros(d+1),np.identity(d+1))
In [5]:
# MCMC

# On essaie avec une variance d'ordre 1, mais trop peu d'acceptations. 
# On modifie alors la variance de manière à viser un taux d'acceptation de 25%
In [6]:
def MCMCSample(Nstep):
    
    b_old = np.random.multivariate_normal(np.zeros(d+1),np.identity(d+1))   
    count = 0

    for i in range(Nstep):
        b_new = b_old + np.random.multivariate_normal(np.zeros(d+1),np.identity(d+1)/3)  # check variance
        if np.random.rand() < UNposterior(b_new,X,Y) / UNposterior(b_old,X,Y):
            b_old = b_new
            count +=1
    return(b_new,count)
In [7]:
MCMCSample(200)
Out[7]:
(array([ 1.11815872, -0.6443458 , -1.24872617,  0.42942326, -0.76397127,
         0.52680061,  2.21922648, -1.55142666]), 19)
In [8]:
def MHestimator(Nstep):
    
    b_old = np.random.multivariate_normal(np.zeros(d+1),np.identity(d+1))   
    S = np.zeros(d+1)
    
    for i in range(Nstep):
        b_new = b_old + np.random.multivariate_normal(np.zeros(d+1),np.identity(d+1)/3)
        if np.random.rand() < UNposterior(b_new,X,Y) / UNposterior(b_old,X,Y):
            b_old = b_new
        S += b_old
    
    return(S/Nstep)
In [9]:
MHestimator(10000)
Out[9]:
array([-0.01202576,  0.11415232, -0.62951689,  0.85977793, -0.75911585,
        0.0861205 ,  1.16950562, -0.29390117])
In [10]:
# Prediction de la proba de succès de Y pour un nouvel individu X

Xnew = np.random.rand(d)
bEst = MHestimator(20)
ProbY = sigma(Xnew,bEst)
ProbY
Out[10]:
0.5289922707577207
In [ ]: