Intervalo de confianza para ajuste de curva exponencial

Estoy tratando de obtener un intervalo de confianza en un ajuste exponencial de algunos datos x,y (disponibles aquí ). Aquí está el MWE que tengo que encontrar el mejor ajuste exponencial a los datos:

 from pylab import * from scipy.optimize import curve_fit # Read data. x, y = np.loadtxt('exponential_data.dat', unpack=True) def func(x, a, b, c): '''Exponential 3-param function.''' return a * np.exp(b * x) + c # Find best fit. popt, pcov = curve_fit(func, x, y) print popt # Plot data and best fit curve. scatter(x, y) x = linspace(11, 23, 100) plot(x, func(x, *popt), c='r') show() 

que produce:

introduzca la descripción de la imagen aquí

¿Cómo puedo obtener el intervalo de confianza del 95% (o algún otro valor) en este ajuste, preferiblemente usando python puro, numpy o scipy (cuáles son los paquetes que ya tengo instalados)?

La respuesta de Gabriel es incorrecta. Aquí, en rojo, la banda de confianza del 95% para sus datos, calculada por GraphPad Prism: Prisma de confianza y bandas de predicción.

Antecedentes: el “intervalo de confianza de una curva ajustada” generalmente se denomina banda de confianza . Para una banda de confianza del 95%, uno puede estar seguro al 95% de que contiene la curva verdadera. (Esto es diferente de las bandas de predicción , que se muestran arriba en gris. Las bandas de predicción se refieren a puntos de datos futuros. Para obtener más detalles, consulte, por ejemplo, esta página de la Guía de ajuste de curvas GraphPad).

En Python, kmpfit puede calcular la banda de confianza para los mínimos cuadrados no lineales. Aquí para el ejemplo de Gabriel:

 from pylab import * from kapteyn import kmpfit x, y = np.loadtxt('_exp_fit.txt', unpack=True) def model(p, x): a, b, c = p return a*np.exp(b*x)+c f = kmpfit.simplefit(model, [.1, .1, .1], x, y) print f.params # confidence band a, b, c = f.params dfdp = [np.exp(b*x), a*x*np.exp(b*x), 1] yhat, upper, lower = f.confidence_band(x, dfdp, 0.95, model) scatter(x, y, marker='.', s=10, color='#0000ba') ix = np.argsort(x) for i, l in enumerate((upper, lower, yhat)): plot(x[ix], l[ix], c='g' if i == 2 else 'r', lw=2) show() 

Los dfdp son las derivadas parciales ∂f / ∂p del modelo f = a * e ^ (b * x) + c con respecto a cada parámetro p (es decir, a, b, y c). Para obtener más información, consulte el Tutorial de kmpfit o esta página de la Guía de ajuste de curvas GraphPad. (A diferencia de mi código de muestra, el Tutorial de kmpfit no utiliza el depósito de confidence_band() de la biblioteca, sino su propia implementación, ligeramente diferente).

Finalmente, la ttwig de Python coincide con la de Prism:

bandas de confianza kmpfit

Puede usar el módulo de incertidumbres para hacer los cálculos de incertidumbre. uncertainties un seguimiento de las incertidumbres y la correlación. Puede crear uncertainties.ufloat correlacionadas.ufloat directamente desde la salida de curve_fit .

Para poder realizar esos cálculos en operaciones no integradas, como exp , debe utilizar las funciones desde uncertainties.unumpy .

También debe evitar su from pylab import * . Esto incluso sobrescribe las incorporaciones de Python, como la sum .

Un ejemplo completo:

 import numpy as np from scipy.optimize import curve_fit import uncertainties as unc import matplotlib.pyplot as plt import uncertainties.unumpy as unp def func(x, a, b, c): '''Exponential 3-param function.''' return a * np.exp(b * x) + c x, y = np.genfromtxt('data.txt', unpack=True) popt, pcov = curve_fit(func, x, y) a, b, c = unc.correlated_values(popt, pcov) # Plot data and best fit curve. plt.scatter(x, y, s=3, linewidth=0, alpha=0.3) px = np.linspace(11, 23, 100) # use unumpy.exp py = a * unp.exp(b * px) + c nom = unp.nominal_values(py) std = unp.std_devs(py) # plot the nominal value plt.plot(px, nom, c='r') # And the 2sigma uncertaintie lines plt.plot(px, nom - 2 * std, c='c') plt.plot(px, nom + 2 * std, c='c') plt.savefig('fit.png', dpi=300) 

Y el resultado: resultado

Aviso : la respuesta real para obtener el intervalo de confianza de la curva ajustada la proporciona Ulrich aquí .


Después de algunas investigaciones (vea aquí , aquí y 1.96 ) se me ocurrió mi propia solución.

Acepta un intervalo de confianza de X% arbitrario y traza curvas superiores e inferiores.

introduzca la descripción de la imagen aquí

Aquí está el MWE:

 from pylab import * from scipy.optimize import curve_fit from scipy import stats def func(x, a, b, c): '''Exponential 3-param function.''' return a * np.exp(b * x) + c # Read data. x, y = np.loadtxt('exponential_data.dat', unpack=True) # Define confidence interval. ci = 0.95 # Convert to percentile point of the normal distribution. # See: https://en.wikipedia.org/wiki/Standard_score pp = (1. + ci) / 2. # Convert to number of standard deviations. nstd = stats.norm.ppf(pp) print nstd # Find best fit. popt, pcov = curve_fit(func, x, y) # Standard deviation errors on the parameters. perr = np.sqrt(np.diag(pcov)) # Add nstd standard deviations to parameters to obtain the upper confidence # interval. popt_up = popt + nstd * perr popt_dw = popt - nstd * perr # Plot data and best fit curve. scatter(x, y) x = linspace(11, 23, 100) plot(x, func(x, *popt), c='g', lw=2.) plot(x, func(x, *popt_up), c='r', lw=2.) plot(x, func(x, *popt_dw), c='r', lw=2.) text(12, 0.5, '{}% confidence interval'.format(ci * 100.)) show() 

curve_fit() devuelve la matriz de covarianza, pcov, que contiene las incertidumbres estimadas (1 sigma). Esto supone que los errores se distribuyen normalmente, lo que a veces es cuestionable.

También puede considerar el uso del paquete lmfit (python puro, construido sobre scipy), que proporciona un envoltorio alrededor de las rutinas de ajuste scipy.optimize (incluyendo leastsq (), que es lo que usa curve_fit ()) y puede, entre otras cosas, calcular intervalos de confianza explícitamente.

Siempre me han gustado los bootstrapping simples para obtener intervalos de confianza. Si tiene n puntos de datos, use el paquete random para seleccionar n puntos de sus datos CON RESAMPLIFICACIÓN (es decir, permita que su progtwig obtenga el mismo punto varias veces si eso es lo que quiere hacer, muy importante). Una vez que hayas hecho eso, traza los puntos remuestreados y obtén el mejor ajuste. Haz esto 10,000 veces, obteniendo una nueva línea de ajuste cada vez. Luego, su intervalo de confianza del 95% es el par de líneas que encierran el 95% de las mejores líneas de ajuste que realizó.

Es un método bastante fácil de progtwigr en Python, pero no está claro cómo funcionaría esto desde un punto de vista estadístico. Más información sobre por qué quiere hacer esto probablemente conduciría a respuestas más apropiadas para su tarea.