Calcolo della banda di confidenza di adattamento minimo quadrato

Question

Ho una domanda con cui mi batto per giorni con ora.

Come si calcola la banda di confidenza (95%) di un adattamento?

adattamento delle curve ai dati è il lavoro di tutti i giorni di ogni fisico - quindi penso che questo dovrebbe essere attuata da qualche parte - ma non riesco a trovare un'implementazione per questo nemmeno io so come farlo matematicamente .

L'unica cosa che ho trovato è seaborn che fa un buon lavoro per lineare minimum-square.

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") 

Ma questo è solo lineare dei minimi quadrati. Quando voglio adattarmi, ad es. una curva di saturazione come , sono fregato.

Certo, posso calcolare la t-distribuzione dall'errore std di un metodo di minimo numero come scipy.optimize.curve_fit ma non è quello che sto cercando.

Grazie per qualsiasi aiuto !!

Answer 1

È possibile ottenere questo facilmente utilizzando il modulo StatsModels.

Vedere anche this example e this answer.

Ecco una risposta per la tua domanda:

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() 

Answer 2

s' confidence_band()kmpfit calcola la band di fiducia per non lineari minimi quadrati. Qui per visualizzare la curva di saturazione:

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') 

Il dfdp sono le derivate parziali ∂f/∂p del modello f = a * (1-e^(b * x)) rispetto a ciascun parametro p (cioè , a e b), vedere il mio answer a una domanda simile per i link in background. Ed ecco l'output: