Exemple de comment calculer une fonction de vraisemblance (log-likelihood) en python pour des données distribuées selon une loi normale:
Pour comprendre l'intérêt de savoir calculer le "log-likelihood" pour estimer une moyenne, voir Comment estimer une moyenne avec des données tronquées (censurées) en python (exemple avec une loi normale) ?
.
1 -- Générer des nombres aléatoires depuis une loi normale
Générons par exemple 100000 depuis une loi normale avec $\mu_0 = 3$ et $\sigma = 0.5$
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats
mu = 3.0
sigma = 0.5
data = np.random.randn(100000) * sigma + mu
2 -- Visualiser les données
hx, hy, _ = plt.hist(data, bins=50, normed=1,color="lightblue")
plt.ylim(0.0,max(hx)+0.05)
plt.title(r'Normal distribution $\mu_0 = 3$ and $\sigma_0 = 0.5$')
plt.grid()
plt.savefig("likelihood_normal_distribution_01.png", bbox_inches='tight')
#plt.show()
plt.close()
3 -- Calculer la fonction de vraisemblance
scipy.stats.norm.pdf(6,2.0,1.0)
print( np.log(scipy.stats.norm.pdf(data,2.0,1.0)).sum() )
x = np.linspace(-10, 10, 1000, endpoint=True)
y = []
for i in x:
y.append(np.log(scipy.stats.norm.pdf(data,i,0.5)).sum())
plt.plot(x,y)
plt.title(r'Log-Likelihood')
plt.xlabel(r'$\mu$')
plt.grid()
plt.savefig("likelihood_normal_distribution_02.png", bbox_inches='tight')
#plt.show()
3 -- Trouver la moyenne
Moyenne estimée avec numpy:
print('mean ---> ', np.mean(data))
print('std deviation ---> ', np.std(data))
donne par exemple
mean ---> 3.0009174745755143
std deviation ---> 0.49853007155264806
Moyenne estimée a partir du maximum du log-likelihood:
y_min = y.index(max(y))
print('mean (from max log likelohood) ---> ', x[y_min])
donne par exemple:
mean (from max log likelohood) ---> 2.9929929929929937