Error estándar en la regresión no lineal

Question

He estado haciendo algunas simulaciones físicas de Monte Carlo con Python y no puedo determinar el error estándar para los coeficientes de un ajuste no lineal de mínimos cuadrados.

Inicialmente, estaba usando SciPy's scipy.stats.linregress para mi modelo, ya que pensé que sería un modelo lineal, pero noté que en realidad es una especie de función de alimentación. Luego utilicé NumPy's polyfit con los grados de libertad siendo 2 pero no puedo encontrar de todos modos para determinar el error estándar de los coeficientes.

Sé que gnuplot puede determinar los errores por mí, pero tengo que hacer ajustes para más de 30 casos diferentes. Me preguntaba si alguien sabe de todos modos si Python lee el error estándar de gnuplot o si hay alguna otra biblioteca que pueda usar.

Answer 1

se ve como gnuplot utiliza levenberg - marquardt y hay una pitón implementation available - puede obtener las estimaciones de error del atributo mpfit.covar (dicho sea de paso, debe preocuparse por lo que las estimaciones de error "significan" - son otros parámetros permitidos para ajustar para compensar, por ejemplo?)

Answer 2

¡Finalmente encontré la respuesta a esta pregunta tan larga! Espero que esto al menos le ahorre a alguien unas pocas horas de investigación sin esperanza para este tema. Scipy tiene una función especial llamada curve_fit en su sección de optimización. Utiliza el método de mínimos cuadrados para determinar los coeficientes y, lo mejor de todo, le da la matriz de covarianza. La matriz de covarianza contiene la varianza de cada coeficiente. Más exactamente, la diagonal de la matriz es la varianza y por enraizamiento cuadrado de los valores, ¡se puede determinar el error estándar de cada coeficiente! Scipy no tiene mucha documentación para esto, así que aquí está un código de ejemplo para una mejor comprensión:

import numpy as np 
from scipy.optimize import curve_fit 
import matplotlib.pyplot as plot 


def func(x,a,b,c): 
    return a*x**2 + b*x + C#Refer [1] 

x = np.linspace(0,4,50) 
y = func(x,2.6,2,3) + 4*np.random.normal(size=len(x)) #Refer [2] 


coeff, var_matrix = curve_fit(func,x,y) 
variance = np.diagonal(var_matrix) #Refer [3] 

SE = np.sqrt(variance) #Refer [4] 

#======Making a dictionary to print results======== 
results = {'a':[coeff[0],SE[0]],'b':[coeff[1],SE[1]],'c':[coeff[2],SE[2]]} 

print "Coeff\tValue\t\tError" 
for v,c in results.iteritems(): 
    print v,"\t",c[0],"\t",c[1] 
#========End Results Printing================= 

y2 = func(x,coeff[0],coeff[1],coeff[2]) #Saves the y values for the fitted model 

plot.plot(x,y) 
plot.plot(x,y2) 

plot.show() 

1. Lo que esta función devuelve es fundamental, ya que define lo que se utiliza para encajar para el modelo
2. Utilizando el función para crear algunos datos arbitrarios + algo de ruido
3. países de diagonal a una matriz 1D la matriz de covarianza que es sólo una matriz normal de
4. Square enraizamiento la varianza para obtener el error estándar (sE)