Previous content: Maximum likelihood estimation of discrete probability distribution and model fitting This time, we will apply the maximum likelihood estimation of the continuous probability distribution and the model (probability distribution).
The probability function that represents the normal distribution is as follows.
The test data for estimating and fitting the model is created in advance from the probability distribution of the target so that the answers can be matched.
Python
"""Random numbers from the normal distribution"""
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(seed=10)
norm_values = np.random.normal(0, 1, size=1000) #Average 0,1000 extracts from normal distribution with standard deviation 1
#drawing
plt.hist(norm_values, bins=50, range=[-5, 5], normed=True)
bin_edges
Since it is a continuous variable, it cannot be counted like a discrete type. Therefore, it is counted by dividing it into 0.2.
Now, let's apply a normal distribution to such data and estimate its parameters with maximum likelihood.
Python
"""Parameter estimation by partial differentiation of the log-likelihood function"""
import sympy
#Define variables (v=σ**2)
sympy.var('μ v y')
#Likelihood p(Parameters|x)Define
fe=(1/sympy.sqrt(2*sympy.pi*v))*sympy.exp(-(y-μ)**2/(2*v))
#Logarithmic
logf=sympy.log(fe)
#Partially differentiate f and expand the equation
pdff1 = sympy.expand(sympy.diff(logf, μ)) #Partial differential with respect to μ
pdff2 = sympy.expand(sympy.diff(logf, v)) #Partial differential with respect to v
The formula pdff1 that partially differentiates the log-likelihood with respect to μ is $ \ frac {y} {v}-\ frac {\ mu} {v} $, and the formula pdff2 that partially differentiates the log-likelihood with respect to v is $ \ frac {-1} {2v } + \ frac {y ^ 2} {2v ^ 2}-\ frac {y \ mu} {v ^ 2} + \ frac {\ mu ^ 2} {2v ^ 2} $ That's right. Since $ \ sum pdff1 = \ sum pdff2 = 0 $ μ, v can be found.
Python
def L_sympy(fmu,fs,var,values):
likelihood_mu = 0 #Initial value of likelihood
likelihood_s = 0 #Initial value of likelihood
for i in np.arange(len(values)):
# likelihood
likelihood_mu += fmu.subs(var,values[i]) #Assign a value to x
likelihood_s += fs.subs(var,values[i]) #Assign a value to x
param = sympy.solve([likelihood_mu,likelihood_s]) #Solve the equation
return param
parameters = L_sympy(pdff1,pdff2,"y",norm_values)
parameters[0]["σ"]=sympy.sqrt(parameters[0][v])
parameters
[{v: 0.879764754284410, μ: -0.0145566356154705, 'σ': 0.937957757196138}]
It's almost the same as the average μ = 0 and σ = 1.
For details, see the previous ①. Use the method of solving nonlinear equations.
Python
"""Probability density function"""
def probability_function(x,param):
from scipy.special import factorial
# param[0]s,param[1]mu:Parameters
# y:Data y
# return:Probability density P of data y(y|Parameters)
return (1/np.sqrt(2*np.pi*param[1]**2))*np.exp(-0.5*(y-param[0])**2/param[1]**2)
"""Log-likelihood function (continuous)"""
def L_func_c(param,y):
likelihood = 0 #Initial value of likelihood
for i in np.arange(len(y)):
# model output
p = probability_function(y[i], param)
#likelihood likelihood
likelihood += -np.log(p) #Likelihood
return likelihood
"""Parameter estimation"""
"""
Quasi-Newton method (BFGS, L)-BFGS method: Memory can be saved in complicated cases)
"""
from scipy import optimize
x0 = [0,0.1] #Initial values of parameters
bound = [(-100, 100),(0, None)]#,(0,None)] #Scope of optimal parameter search(min,max)
params_MLE = optimize.minimize(L_func_c,x0,args=(norm_values),method='l-bfgs-b',
jac=None, bounds=bound, tol=None, callback=None,
options={'disp': None, 'maxls': 20, 'iprint': -1,
'gtol': 1e-05, 'eps': 1e-08, 'maxiter': 15000,
'ftol': 2.220446049250313e-09, 'maxcor': 10,
'maxfun': 15000})
#Maximum likelihood estimation parameters
print('Parameter μ,σ:',params_MLE.x)
#Number of parameters
k=3
#AIC (the model with the smallest value is a good fit)
print('AIC:',params_MLE.fun*(2)+2*k)
Parameter μ,σ: [-0.01455655 0.93795781]
AIC: 2715.7763344850605
For details, see the previous ① . Use the method of solving nonlinear equations.
Python
"""scipy.stats.There is also a parameter estimation with fit"""
from scipy.stats import norm
fit_parameter = norm.fit(norm_values)
fit_parameter
#Parameter μ,σ
(-0.014556635615470447, 0.9379577571961389)
Let's apply the model (normal distribution) to the data with the estimated parameters μ = -0.01455655, σ = 0.93795781.
Python
acc_mle = probability_function(np.sort(norm_values), params_MLE.x)
plt.hist(norm_values, bins=50, range=[-5, 5], normed=True)
plt.plot(np.sort(norm_values), acc_mle)
It feels good.
Recommended Posts