¿Cómo calcular la probabilidad de ajuste de curvas en scipy?

Question

Tengo un ajuste de modelo no lineal que se ve así:

La línea continua oscura es el ajuste del modelo y la parte gris son los datos sin procesar.

Versión corta de la pregunta:¿Cómo obtengo la probabilidad de que este modelo se ajuste para poder realizar una prueba de relación de probabilidad logarítmica?Suponga que el residuo se distribuye normalmente.

Soy relativamente nuevo en las estadísticas y mis pensamientos actuales son:

1. Obtenga el residual del ajuste de la curva y calcule la varianza del residual;

2. Usa esta ecuaciónY conecte la varianza del residuo en sigma cuadrado, x_i como experimento y mu como ajuste del modelo;

3. Calcule la relación logarítmica de verosimilitud.

¿Alguien podría ayudarme con estas dos preguntas de la versión completa?

1. ¿Es correcto mi método?(Creo que sí, ¡pero sería genial estar seguro!)

2. ¿Hay funciones listas para usar en python/scipy/statsmodels para hacer esto por mí?

Answer 1

Tu función de probabilidad

que es simplemente la suma del logaritmo de la función de densidad de probabilidad de la distribución gaussiana.

es la probabilidad de Colocando un mu y un sigma para sus residuos., no la probabilidad de tu modelo dados tus datos.En una palabra, su enfoque es equivocado.

Dado que estás haciendo mínimos cuadrados no lineales, siguiendo lo que @usethedeathstar ya mencionó, deberías ir directamente a F-test..Considere el siguiente ejemplo, modificado de http://www.walkingrandomly.com/?p=5254, y llevamos a cabo F-test usando R.Y discutiremos cómo traducirlo al python al final.

# 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

aquí tenemos dos modelos, fit1 tiene 2 parámetros, por lo tanto el residuo tiene 8 grados de libertad; fit2 tiene un parámetro adicional y el residuo tiene 7 grados de libertad.¿Es el modelo 2 significativamente mejor?No, el valor F es 0,9194, en (1,7) grados de libertad y no es significativo.

Para obtener la tabla ANOVA:El DF residual es fácil.Residuo Suma de cuadrados: 0.08202*0.08202*8=0.05381 y 0.08243*0.08243*7=0.04756293 (aviso: 'Error estándar residual:0,08243 en 7 grados de libertad', etc).En python, puedes conseguirlo por (y_observed-y_fitted)**2, desde scipy.optimize.curve_fit() no devuelve los residuos.

El F-ratio es 0.0062473/0.047565*7 y para obtener el valor P: 1-scipy.stats.f.cdf(0.9194, 1, 7).

Ponlos juntos tenemos python equivalente:

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

Como señaló @usethedeathstar:cuando el residuo se distribuye normalmente, mínimos cuadrados no lineales ES la máxima probabilidad.Por lo tanto, la prueba F y la prueba de razón de verosimilitud son equivalentes.Porque, La relación F es una transformación monótona de la relación de verosimilitud λ.

O de forma descriptiva, ver: http://www.stata.com/support/faqs/statistics/chi-squared-and-f-distributions/

Answer 2

Tu fórmula se ve correcta para mí.Debe darle los mismos resultados que los scipy.stats.norm.logpdf(x, loc=mu, scale=sigma)

Dado que ya tiene sus estimaciones de MU y SIGMA, no creo que haya una función para la prueba de relación de probabilidad en la que pueda conectar sus resultados.

Si tiene las estimaciones de dos modelos, donde uno está anidado en el otro, entonces puede calcularlo fácilmente.

http://en.wikipedia.org/wiki/likeliud-ratio_test

Aquí está la parte de un método en StatSModels que calcula la prueba de LR para comparar dos modelos lineales anidados https://github.com/statsmodels/statsmodels/blob/master/statsmodels/regression/ linear_model.py#l1531