[PYTHON] Essayons à nouveau La plupart des estimations de probabilité et ajustement du modèle (distribution de probabilité) ② Distribution de probabilité continue

introduction

Contenu précédent: Estimation la plus probable de la distribution de probabilité discrète et de l'ajustement du modèle Cette fois, nous appliquerons l'estimation la plus probable de la distribution de probabilité continue et du modèle (distribution de probabilité).

Distribution de probabilité continue: distribution normale

La fonction de probabilité qui représente la distribution normale est la suivante. $ P (y) = \ frac {1} {\ sqrt {2πσ ^ 2}} e ^ {(- \ frac {(y_i- \ mu) ^ 2} {2 \ sigma ^ 2})} (\ mu , \ sigma = paramètre) $

Les données de test pour estimer et ajuster le modèle sont créées à l'avance à partir de la distribution de probabilité cible afin que les réponses puissent être appariées.

Python

"""Nombres aléatoires de la distribution normale"""
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(seed=10)
norm_values = np.random.normal(0, 1, size=1000) #Moyenne 0,Extraire 1000 pièces d'une distribution normale avec un écart type de 1
#dessin
plt.hist(norm_values, bins=50, range=[-5, 5], normed=True)
bin_edges

Puisqu'il s'agit d'une variable continue, elle ne peut pas être comptée comme un type discret. Par conséquent, il est compté en le divisant en 0,2.

Maintenant, appliquons une distribution normale à ces données et estimons ses paramètres les plus probables.

Python

"""Estimation des paramètres par différenciation partielle de la fonction de vraisemblance logarithmique"""
import sympy
#Définir des variables (v=σ**2)
sympy.var('μ v y')
#Responsabilité p(Paramètres|x)Définir
fe=(1/sympy.sqrt(2*sympy.pi*v))*sympy.exp(-(y-μ)**2/(2*v))
#Logistique
logf=sympy.log(fe)
#Différencier partiellement f et développer l'équation
pdff1 = sympy.expand(sympy.diff(logf, μ)) #Différenciation partielle pour μ
pdff2 = sympy.expand(sympy.diff(logf, v)) #Différenciation partielle pour v

La formule pdff1 qui différencie partiellement la vraisemblance logarithmique par rapport à μ est $ \ frac {y} {v} - \ frac {\ mu} {v} $, et la formule pdff2 qui est partiellement différenciée par rapport à v est $ \ frac {-1} {2v. } + \ frac {y ^ 2} {2v ^ 2} - \ frac {y \ mu} {v ^ 2} + \ frac {\ mu ^ 2} {2v ^ 2} $ C'est vrai. $ \ sum pdff1 = \ sum pdff2 = 0 $ car vous pouvez trouver μ, v

Python

def L_sympy(fmu,fs,var,values):       
    likelihood_mu = 0 #Valeur initiale de la vraisemblance
    likelihood_s = 0 #Valeur initiale de la vraisemblance
    for i in np.arange(len(values)):
        # likelihood
        likelihood_mu += fmu.subs(var,values[i]) #Remplacez x par une valeur
        likelihood_s += fs.subs(var,values[i]) #Remplacez x par une valeur
    param = sympy.solve([likelihood_mu,likelihood_s]) #Résous l'équation
    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}]

C'est presque le même que le réglage de la moyenne μ = 0 et σ = 1.

Quand le modèle statistique se complique

Pour plus de détails, voir le précédent ①. Utilisez la méthode de résolution d'une équation non linéaire.

Python

"""Fonction de densité de probabilité"""
def probability_function(x,param):
    from scipy.special import factorial
    # param[0]s,param[1]mu:Paramètres
    # y:Données y
    # return:Densité de probabilité P des données y(y|Paramètres) 
    return (1/np.sqrt(2*np.pi*param[1]**2))*np.exp(-0.5*(y-param[0])**2/param[1]**2)

"""Fonction de vraisemblance logarithmique (continue)"""
def L_func_c(param,y):       
    likelihood = 0 #Valeur initiale de la vraisemblance
    for i in np.arange(len(y)):
        # model output 
        p = probability_function(y[i], param)
        #vraisemblance vraisemblance
        likelihood += -np.log(p) #Responsabilité
    return likelihood

"""Estimation des paramètres"""
"""
Méthode Semi-Newton (BFGS, L)-Méthode BFGS: la mémoire peut être sauvegardée dans les cas compliqués)
"""
from scipy import optimize
x0 = [0,0.1] #Valeur initiale du paramètre
bound = [(-100, 100),(0, None)]#,(0,None)] #Portée de la recherche de paramètres optimale(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})
#Paramètre d'estimation le plus probable
print('Paramètre μ,σ：',params_MLE.x)
#Nombre de paramètres
k=3
#AIC (le modèle avec la plus petite valeur est un bon ajustement)
print('AIC：',params_MLE.fun*(2)+2*k)

Paramètre μ,σ： [-0.01455655  0.93795781]
AIC： 2715.7763344850605

Je veux connaître les paramètres rapidement ...

Pour plus de détails, voir le ① précédent . Utilisez la méthode de résolution d'une équation non linéaire.

Python

"""scipy.stats.Il existe également une estimation des paramètres avec ajustement"""
from scipy.stats import norm
fit_parameter = norm.fit(norm_values)
fit_parameter

#Paramètre μ,σ
(-0.014556635615470447, 0.9379577571961389)

Appliquer aux données et essayer d'illustrer

Appliquons le modèle (distribution normale) aux données avec les paramètres estimés μ = -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)

Ça fait du bien.