Estimation du maximum de vraisemblance avec la loi normale

Maximum de vraisemblance

\begin{equation}
L(\theta_1,\theta_2,\theta_3,...,\theta_m) = \prod_{i=1}^n f(x_i;\theta_1,\theta_2,\theta_3,...,\theta_m)
\end{equation}

Loi normale

\begin{equation}
f(x_i,\theta_1,\theta_2) = \frac{1}{\sqrt{2\pi\theta_2}}exp\big[-\frac{(x_i-\theta_1)^2}{2\theta_2}\big]
\end{equation}

\begin{equation}
L(\theta_2,\theta_2)=\prod_{i=1}^n f(x_i,\theta_1,\theta_2) = \theta_2^{-n/2} (2\pi)^{-n/2} exp \big[ -\frac{1}{2\theta_2}\sum_{i=1}^n(x_i-\theta_1)^2 \big]
\end{equation}

\begin{equation}
log L(\theta_2,\theta_2) = - \frac{n}{2}log\theta_2 - \frac{n}{2}log(2\pi)-\frac{\sum(x_i-\theta_1)^2}{2\theta_2}
\end{equation}

\begin{equation}
\require{cancel}
\frac{\partial log L(\theta_1,\theta_2) }{\partial \theta_1} = \frac{\cancel{-2}\sum(x_i-\theta_1)\cancel{(-1)}}{\cancel{2}\theta_2}=0
\end{equation}

\begin{equation}
\sum x_i - n\theta_1 = 0
\end{equation}

\begin{equation}
\boxed{ \hat{\theta}_1 = \hat{\mu} = \frac{\sum x_i}{n} = \bar{x} }
\end{equation}

\begin{equation}
\require{cancel}
\frac{\partial log L(\theta_1,\theta_2) }{\partial \theta_2} = -\frac{n}{2\theta_2}+\frac{\sum(x_i-\theta_1)^2}{2\theta_2^2}=0
\end{equation}

\begin{equation}
-n\theta_2+\sum(x_i-\theta_1)^2 = 0
\end{equation}

\begin{equation}
\boxed{ \hat{\theta}_2 = \frac{\sum(x_i-\bar{x})^2}{n} }
\end{equation}

Code python

from random import randint

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats

#----------------------------------------------------------------------------------------#
# data

nb_data = 100

mu, sigma = 10, 2

x = np.linspace(0,20,1000)

plt.plot(x,stats.norm.pdf(x,mu,sigma))

x = np.random.normal(mu, sigma, nb_data)
y = np.zeros((nb_data))

plt.scatter(x,y)
plt.title('Random Numbers from a Normal Distribution \n (10,2)')

#plt.savefig('normal_distribution_random_numbers.png', bbox_inches='tight')
plt.show()

plt.close()

#----------------------------------------------------------------------------------------#
# theta 1 and theta 2 estimate

theta1_estimate = 0
theta2_estimate = 0

for i in range(nb_data):
    theta1_estimate = theta1_estimate + x[i]

theta1_estimate = theta1_estimate / nb_data

for i in range(nb_data):
    theta2_estimate = theta2_estimate + ( x[i] - theta1_estimate )**2

theta2_estimate = theta2_estimate / nb_data

#----------------------------------------------------------------------------------------#
# Maximum Likelihood (known variance)

theta2 = sigma ** 2

theta1_list = []
log_likelihood_list = []

for theta1 in np.arange(5,15,0.1):

    part1 = - (nb_data / 2.0) * np.log(theta2)
    part2 = - (nb_data / 2.0) * np.log(2.0*np.pi)

    part3 = 0
    for j in range(nb_data):
        part3 = part3 + ( x[j] - theta1 )**2
    part3 = - part3 / (2.0 * theta2 )

    theta1_list.append(theta1)
    log_likelihood_list.append(part1 + part2 + part3)

f = plt.figure()
ax = f.add_subplot(111)

plt.axvline(x=mu)

plt.plot(theta1_list,log_likelihood_list)
plt.title('Log-Likelihood')

plt.text(0.65,0.5,r'$\theta_1=$'+str(round(theta1_estimate,2)),horizontalalignment='center',
     verticalalignment='center', transform = ax.transAxes)

#plt.savefig('maximum_likelihood_gaussian_distribution_1.png', bbox_inches='tight')
plt.show()

Recherches associées