Example of how to calculate a log-likelihood using a normal distribution in python:
See the note: How to estimate the mean with a truncated dataset using python ? to understand the interest of calculating a log-likelihood using a normal distribution in python.
1 -- Generate random numbers from a normal distribution
Let's for example create a sample of 100000 random numbers from a normal distribution of mean $\mu_0 = 3$ and standard deviation $\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 -- Plot the data
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 -- Calculate the log-likelihood
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 -- Find the mean
Mean estimation using numpy:
print('mean ---> ', np.mean(data))
print('std deviation ---> ', np.std(data))
returns for example
mean ---> 3.0009174745755143
std deviation ---> 0.49853007155264806
Mean estimated from the maximum of the log-likelihood:
y_min = y.index(max(y))
print('mean (from max log likelohood) ---> ', x[y_min])
returns for example
mean (from max log likelohood) ---> 2.9929929929929937