Calcular la banda de confianza del ajuste por mínimos cuadrados

Tengo una pregunta con la que peleo por días con ahora.

¿Cómo calculo la banda de confianza (95%) de un ajuste?

Ajustar las curvas a los datos es el trabajo diario de todos los físicos, por lo que creo que esto debería implementarse en algún lugar, pero no puedo encontrar una implementación para esto, tampoco sé cómo hacerlo matemáticamente.

Lo único que encontré es el de los seaborn que hace un buen trabajo para el cuadrado mínimo lineal .

 import numpy as np from matplotlib import pyplot as plt import seaborn as sns import pandas as pd x = np.linspace(0,10) y = 3*np.random.randn(50) + x data = {'x':x, 'y':y} frame = pd.DataFrame(data, columns=['x', 'y']) sns.lmplot('x', 'y', frame, ci=95) plt.savefig("confidence_band.pdf") 

lineal-menos cuadrado

Pero esto es solo un cuadrado mínimo lineal. Cuando quiero ajustar, por ejemplo, una curva de saturación como ecuación de saturación Estoy jodido

Claro, puedo calcular la distribución t a partir del error scipy.optimize.curve_fit de un método de mínimos cuadrados como scipy.optimize.curve_fit pero eso no es lo que estoy buscando.

¡¡Gracias por cualquier ayuda!!

Puedes lograrlo fácilmente usando el módulo StatsModels .

También vea este ejemplo y esta respuesta .

Aquí hay una respuesta para su pregunta:

 import numpy as np from matplotlib import pyplot as plt import statsmodels.api as sm from statsmodels.stats.outliers_influence import summary_table x = np.linspace(0,10) y = 3*np.random.randn(50) + x X = sm.add_constant(x) res = sm.OLS(y, X).fit() st, data, ss2 = summary_table(res, alpha=0.05) fittedvalues = data[:,2] predict_mean_se = data[:,3] predict_mean_ci_low, predict_mean_ci_upp = data[:,4:6].T predict_ci_low, predict_ci_upp = data[:,6:8].T fig, ax = plt.subplots(figsize=(8,6)) ax.plot(x, y, 'o', label="data") ax.plot(X, fittedvalues, 'r-', label='OLS') ax.plot(X, predict_ci_low, 'b--') ax.plot(X, predict_ci_upp, 'b--') ax.plot(X, predict_mean_ci_low, 'g--') ax.plot(X, predict_mean_ci_upp, 'g--') ax.legend(loc='best'); plt.show() 

Ejemplo

El trust_band confidence_band() kmpfit calcula la banda de confianza para los mínimos cuadrados no lineales . Aquí para su curva de saturación:

 from pylab import * from kapteyn import kmpfit def model(p, x): a, b = p return a*(1-np.exp(b*x)) x = np.linspace(0, 10, 100) y = .1*np.random.randn(x.size) + model([1, -.4], x) fit = kmpfit.simplefit(model, [.1, -.1], x, y) a, b = fit.params dfdp = [1-np.exp(b*x), -a*x*np.exp(b*x)] yhat, upper, lower = fit.confidence_band(x, dfdp, 0.95, model) scatter(x, y, marker='.', color='#0000ba') for i, l in enumerate((upper, lower, yhat)): plot(x, l, c='g' if i == 2 else 'r', lw=2) savefig('kmpfit confidence bands.png', bbox_inches='tight') 

Los dfdp son las derivadas parciales ∂f / ∂p del modelo f = a * (1-e ^ (b * x)) con respecto a cada parámetro p (es decir, a y b), vea mi respuesta a una pregunta similar para enlaces de fondo. Y aquí la salida:

bandas de confianza kmpfit