Comment calculer la probabilité d'ajustement de courbe dans Scipy ?

Question

J'ai un ajustement de modèle non linéaire qui ressemble à ceci :

La ligne continue sombre correspond à l'ajustement du modèle et la partie grise correspond aux données brutes.

Version courte de la question :comment puis-je obtenir la probabilité d'ajustement de ce modèle, afin de pouvoir effectuer un test de rapport de vraisemblance logarithmique ?Supposons que le résidu soit normalement distribué.

Je suis relativement nouveau dans le domaine des statistiques et mes réflexions actuelles sont les suivantes :

1. Obtenez le résidu de l'ajustement de la courbe et calculez la variance du résidu ;

2. Utilisez cette équationEt branchez la variance du résidu dans sigma-carré, x_i comme expérience et mu comme ajustement du modèle ;

3. Calculez le rapport de log-vraisemblance.

Quelqu'un pourrait-il m'aider avec ces deux questions en version complète ?

1. Ma méthode est-elle correcte ?(Je pense que oui, mais ce serait vraiment génial de s'en assurer !)

2. Existe-t-il des fonctions prêtes à l'emploi dans python/scipy/statsmodels pour faire cela pour moi ?

Answer 1

Votre fonction de vraisemblance

qui est simplement la somme du log de la fonction de densité de probabilité de la distribution gaussienne.

est la probabilité de ajuster un mu et un sigma pour votre résidu, pas la probabilité de votre modèle compte tenu de vos données.En un mot, votre approche est faux.

Puisque vous faites des moindres carrés non linéaires, en suivant ce que @usethedeathstar a déjà mentionné, vous devriez aller directement F-test..Considérons l'exemple suivant, modifié à partir de http://www.walkingrandomly.com/?p=5254, et nous effectuons F-test en utilisant R.Et nous discuterons de la façon de le traduire en python à la fin.

# construct the data vectors using c()
> xdata = c(-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9)
> ydata = c(0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001)
# some starting values
> p1 = 1
> p2 = 0.2
> p3 = 0.01

# do the fit
> fit1 = nls(ydata ~ p1*cos(p2*xdata) + p2*sin(p1*xdata), start=list(p1=p1,p2=p2))
> fit2 = nls(ydata ~ p1*cos(p2*xdata) + p2*sin(p1*xdata)+p3*xdata, start=list(p1=p1,p2=p2,p3=p3))

# summarise
> summary(fit1)

Formula: ydata ~ p1 * cos(p2 * xdata) + p2 * sin(p1 * xdata)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
p1 1.881851   0.027430   68.61 2.27e-12 ***
p2 0.700230   0.009153   76.51 9.50e-13 ***
---
Signif. codes:  0 ?**?0.001 ?*?0.01 ??0.05 ??0.1 ??1

Residual standard error: 0.08202 on 8 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 2.189e-06

> summary(fit2)

Formula: ydata ~ p1 * cos(p2 * xdata) + p2 * sin(p1 * xdata) + p3 * xdata

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
p1  1.90108    0.03520  54.002 1.96e-10 ***
p2  0.70657    0.01167  60.528 8.82e-11 ***
p3  0.02029    0.02166   0.937     0.38    
---
Signif. codes:  0 ?**?0.001 ?*?0.01 ??0.05 ??0.1 ??1

Residual standard error: 0.08243 on 7 degrees of freedom

Number of iterations to convergence: 9 
Achieved convergence tolerance: 2.476e-06

> anova(fit2, fit1)
Analysis of Variance Table

Model 1: ydata ~ p1 * cos(p2 * xdata) + p2 * sin(p1 * xdata) + p3 * xdata
Model 2: ydata ~ p1 * cos(p2 * xdata) + p2 * sin(p1 * xdata)
  Res.Df Res.Sum Sq Df     Sum Sq F value Pr(>F)
1      7   0.047565                             
2      8   0.053813 -1 -0.0062473  0.9194 0.3696

ici nous avons deux modèles, fit1 a 2 paramètres, donc le résidu a 8 degrés de liberté ; fit2 a un paramètre supplémentaire et le résidu a 7 degrés de liberté.Le modèle 2 est-il nettement meilleur ?Non, la valeur F est de 0,9194, sur (1,7) degrés de liberté et ce n’est pas significatif.

Pour obtenir le tableau ANOVA :Le DF résiduel est facile.Résidu Somme des carrés : 0.08202*0.08202*8=0.05381 et 0.08243*0.08243*7=0.04756293 (avis: 'Erreur type résiduelle :0,08243 sur 7 degrés de liberté', etc).Dans python, vous pouvez l'obtenir par (y_observed-y_fitted)**2, depuis scipy.optimize.curve_fit() ne renvoie pas les résidus.

Le F-ratio est 0.0062473/0.047565*7 et pour obtenir la valeur P : 1-scipy.stats.f.cdf(0.9194, 1, 7).

Rassemblez-les, nous avons python équivalent:

In [1]:

import scipy.optimize as so
import scipy.stats as ss
xdata = np.array([-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9])
ydata = np.array([0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001])
def model0(x,p1,p2):
    return p1*np.cos(p2*x) + p2*np.sin(p1*x)
def model1(x,p1,p2,p3):
    return p1*np.cos(p2*x) + p2*np.sin(p1*x)+p3*x
p1, p2, p3 = 1, 0.2, 0.01
fit0=so.curve_fit(model0, xdata, ydata, p0=(p1,p2))[0]
fit1=so.curve_fit(model1, xdata, ydata, p0=(p1,p2,p3))[0]
yfit0=model0(xdata, fit0[0], fit0[1])
yfit1=model1(xdata, fit1[0], fit1[1], fit1[2])
ssq0=((yfit0-ydata)**2).sum()
ssq1=((yfit1-ydata)**2).sum()
df=len(xdata)-3
f_ratio=(ssq0-ssq1)/(ssq1/df)
p=1-ss.f.cdf(f_ratio, 1, df)
In [2]:

print f_ratio, p
0.919387419515 0.369574503394

Comme @usethedeathstar l'a souligné :lorsque le résidu est normalement distribué, moindres carrés non linéaires EST le maximum de vraisemblance.Par conséquent, le test F et le test du rapport de vraisemblance sont équivalents.Parce que, Le rapport F est une transformation monotone du rapport de vraisemblance λ.

Ou de manière descriptive, voir : http://www.stata.com/support/faqs/statistics/chi-squared-and-f-distributions/

Answer 2

Votre formule a l'air correcte pour moi.Il devrait vous donner les mêmes résultats que scipy.stats.norm.logpdf(x, loc=mu, scale=sigma)

Puisque vous avez déjà vos estimations de MU et Sigma, je ne pense pas qu'il existe une fonction pour le test de rapport de vraisemblance où vous pouvez brancher vos résultats.

Si vous avez les estimations de deux modèles, où l'on est imbriquée dans l'autre, vous pouvez facilement le calculer vous-même.

http://fr.wikipedia.org/wiki/likélium-ratio_test

Voici la partie d'une méthode dans StatsModels qui calcule le test LR pour comparer deux modèles linéaires imbriqués https://github.com/statsmodels/statsmodels/blob/master/statsmodels/regression/lineear_model.py#l1531