Utiliser l'algorithme espérance-maximisation avec python pour calculer une moyenne sur des données manquantes

Un simple exemple de code écrit en python pour tester le calcul de la moyenne et la déviation standard avec des données manquantes, en utilisant un estimateur du maximum de vraisemblance et l'algorithme espérance-maximisation (en anglais expectation-maximization algorithm, souvent abrégé EM), proposé par Dempster et al. (1977).

Pour rappel, en statistique il est important de savoir si on a affaire à des données manquantes ou pas. Car dans le cas de données manquantes les estimateurs comme:

\begin{equation}
\hat{\mu}=\frac{1}{n}\sum_1^n x_i
\end{equation}

pour calculer la moyenne par exemple est un estimateur fortement biaisé.

Comment faire dans le cas de données manquantes ? une approche possible est de visualiser les données à l'aide d'un histogramme par exemple, puis de regarder si les données suivent une distribution connue (comme une distribution normale) et d'appliquer l'algorithme espérance-maximisation (EM qui consiste en deux étapes: (1): une étape d'évaluation de l'espérance (E), où l'on calcule l'espérance de la vraisemblance en tenant compte des dernières variables observées et (2): une étape de maximisation (M), où l'on estime le maximum de vraisemblance des paramètres en maximisant la vraisemblance trouvée à l'étape E)

Exemple (ci dessus), avec des données manquantes (censurées légèrement à droite: histogramme rouge du milieu) en s'inspirant de la video "Expectation Maximization: how it works". Dans cet exemple on peut supposer (on visualisant l'histogramme) que les données suivent une distribution normale. Pour appliquer l'algorithme espérance-maximisation on commence avec une moyenne et une standard déviation quelconques (par exemple $\mu=10$ et $\sigma=10$), puis on calcule pour chaque points ('E-step'):

\begin{equation}
P(x_i)=\frac{1}{\sqrt{2\pi\sigma^2}}exp(-\frac{(x_i-\mu)^2}{2\sigma^2})
\end{equation}

et on en déduit ('M-step') une nouvelle estimation de $\mu$ et $\sigma$:

\begin{equation}
\hat{\mu}=\frac{\sum P(x_i).x_i}{\sum P(x_i)}
\end{equation}

\begin{equation}
\hat{\sigma}=\frac{\sum P(x_i).(x_i-\hat{\mu})^2}{\sum P(x_i)}
\end{equation}

On recommence alors les étapes E et M pour maximiser l'espérance. Code source écrit avec python 3 pour tester l'algorithme espérance-maximisation:

import numpy as np
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import math

import scipy.stats

#----------------------------------------------------------------------------------------#
# define the truth

mu = 100 
sigma = 25

#----------------------------------------------------------------------------------------#
# create a random sample

nb_point = 10000
x = mu + sigma * np.random.randn(nb_point)

mu_estimate = 0
for i in range(nb_point):
    mu_estimate = mu_estimate + x[i]
mu_estimate = mu_estimate / nb_point

sigma_estimate = 0
for i in range(nb_point):
    sigma_estimate = sigma_estimate + (x[i]-mu_estimate)**2
sigma_estimate = math.sqrt( sigma_estimate / (nb_point - 1) )

print(mu_estimate, sigma_estimate)

num_bins = 50
n, bins, patches = plt.hist(x, num_bins, normed=1, color="royalblue", alpha=0.5)

plt.xlim(mu-4*sigma,mu+4*sigma)

plt.title('Uncensored Data ' + r'$\mu$ = '+'{:2.2f}'.format(mu_estimate)+ \
                                         r'; $\sigma$ = '+'{:2.2f}'.format(sigma_estimate))

plt.savefig("censored_data_1.png")
#plt.show()

#----------------------------------------------------------------------------------------#
# create (right) censored data

x_censored = x[x < mu+1*sigma]

nb_point_censored = x_censored.shape[0]

mu_estimate = 0
for i in range(nb_point_censored):
    mu_estimate = mu_estimate + x_censored[i]
mu_estimate = mu_estimate / nb_point_censored

sigma_estimate = 0
for i in range(nb_point_censored):
    sigma_estimate = sigma_estimate + (x_censored[i]-mu_estimate)**2
sigma_estimate = math.sqrt( sigma_estimate / (nb_point_censored - 1) )

print(mu_estimate, sigma_estimate)

num_bins = 50
n, bins, patches = plt.hist(x_censored, num_bins, normed=1, color="lightcoral", alpha=0.5)

plt.xlim(mu-4*sigma,mu+4*sigma)

labels= ['Uncensored','Right Censored']

plt.legend(labels)

plt.title('Censored Data ' + r'$\mu$ = '+'{:2.2f}'.format(mu_estimate)+ \
                             r'; $\sigma$ = '+'{:2.2f}'.format(sigma_estimate))

plt.savefig("censored_data_2.png")
#plt.show()

#----------------------------------------------------------------------------------------#
# EM

mu_estimate = 10
sigma_estimate = 10

def gaussian_function(x):
    return ( 1.0 / math.sqrt(2.0*math.pi*sigma_estimate**2) ) * math.exp( - (x-mu_estimate)**2 / (2.0*sigma_estimate**2))

print(scipy.stats.norm.pdf(12,mu_estimate,sigma_estimate))
print(gaussian_function(12))

for j in range(10):
    num_expectation = 0.0
    sigma_expectation = 0.0
    normalizing_constant = 0.0
    #----- new estimate of mu -----#
    for i in range(nb_point_censored):
        num_expectation = num_expectation + x_censored[i] * scipy.stats.norm.pdf(x_censored[i],mu_estimate,sigma_estimate)
        sigma_expectation = sigma_expectation + (x_censored[i]-mu_estimate)**2 * scipy.stats.norm.pdf(x_censored[i],mu_estimate,sigma_estimate)
        normalizing_constant = normalizing_constant + scipy.stats.norm.pdf(x_censored[i],mu_estimate,sigma_estimate)
    mu_estimate = num_expectation / normalizing_constant
    #----- new estimate of sigma -----#
    for i in range(nb_point_censored):
        sigma_expectation = sigma_expectation + (x_censored[i]-mu_estimate)**2 * scipy.stats.norm.pdf(x_censored[i],mu_estimate,sigma_estimate) 
    sigma_estimate = math.sqrt( sigma_expectation / normalizing_constant )

print(mu_estimate, sigma_estimate, sigma_estimate**2)

estimated_gaussian_x = list(np.linspace(mu-4*sigma, mu+4*sigma, 100))
estimated_gaussian_y = [scipy.stats.norm.pdf(i,mu_estimate,sigma_estimate) for i in estimated_gaussian_x]

plt.plot(estimated_gaussian_x, estimated_gaussian_y, color='black')

plt.xlim(mu-4*sigma,mu+4*sigma)

plt.title('Censored Data (MLE + EM): ' + r'$\mu$ = '+'{:2.2f}'.format(mu_estimate)+ \
                                         r'; $\sigma$ = '+'{:2.2f}'.format(sigma_estimate))

plt.savefig("censored_data_3.png")

Autre exemple avec des données fortement censurées à droite:

Références