計算最小二乘擬合的置信度帶

Question

我得到了一個問題，我現在已經打了好幾天了。

如何計算擬合的（95％）置信度？

數據擬合曲線是每一個物理學家的每一天的工作 - 所以我認爲這應該是地方實現 - 但我無法找到一個實現這個我也不知道該怎麼做數學這。

我發現的唯一的東西是seaborn，線性好，最少。

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

但這僅僅是線性最小平方。當我想適合例如像這樣的飽和曲線，我擰了。

當然，我可以從像scipy.optimize.curve_fit這樣的最小平方法的std-error中計算t分佈，但這不是我正在尋找的。

感謝您的任何幫助！

Answer 1

您可以使用StatsModels模塊輕鬆實現此功能。

另見this example和this answer。

這裏是你的問題的答案：

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

kmpfit的confidence_band()計算置信帶爲非線性最小二乘。此處爲你的飽和曲線：

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

的dfdp是偏導數模型F = A *（1-E ^（b *的X））相對於每個參數p的∂F/∂p（即，a和b），請參閱我的answer針對背景鏈接的類似問題。而這裏的輸出：