计算最小二乘拟合的置信带。

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

但这只是最小二乘法。当我想要拟合像这样的饱和曲线时，我就束手无策了。

当然，我可以从最小二乘法的标准误差计算t分布，例如使用scipy.optimize.curve_fit，但那不是我要找的。

感谢任何帮助！！

- Suuuehgi

---

https://dev59.com/j14d5IYBdhLWcg3wG_aH#63560689 - Marco Cerliani

Answer 1

你可以轻松地使用 StatsModels 模块来实现这一点。

此外，请参见此示例和这个答案。

以下是你问题的答案：

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

- Seyeong Jeong

---

很遗憾，目前这个功能只在statsmodels中适用于线性函数，并且将在下一个版本中适用于广义线性模型，但对于一般的非线性函数尚不可用。 - Josef

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（即a和b）的偏导数∂f/∂p，详见我回答中的相关链接。以下是输出结果：

- Ulrich Stern