Daidalos May 10, 2020

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