指数曲线拟合的置信区间

Question

我正在尝试为一些 x,y 数据（可以在这里找到）获取一个指数拟合的置信区间。以下是我用来找到最佳指数拟合的数据的最小工作示例（MWE）：

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

这个代码运行后会生成：

我该如何获得这个拟合的95%（或者其他值）的置信区间，最好是使用纯 python、numpy 或 scipy（这些都是我已经安装的包）？

Answer 1

我一直喜欢用简单的方法来获取置信区间。如果你有 n 个数据点，那么可以使用 random 这个包，从你的数据中随机选择 n 个点，并且要允许重复选择（也就是说，程序可以多次选择同一个点，这一点非常重要）。完成这个步骤后，绘制这些重新抽样的点，并找到最佳拟合线。然后重复这个过程 10,000 次，每次都会得到一条新的拟合线。最后，你的 95% 置信区间就是包围住 95% 最佳拟合线的那两条线。

在 Python 中编程实现这个方法相对简单，但从统计学的角度来看，这个方法的效果可能有点不清楚。如果能提供更多关于你为什么要这样做的信息，可能会得到更合适的答案来帮助你完成任务。

Answer 2

curve_fit() 这个函数会返回一个叫做协方差矩阵的东西，简称为 pcov。这个矩阵里面包含了估算的不确定性（也就是1个标准差）。不过，这个计算是基于误差是正态分布的假设，有时候这个假设可能不太靠谱。

你也可以考虑使用 lmfit 这个包（它是纯Python写的，基于scipy），它为scipy.optimize中的拟合方法提供了一个封装（包括leastsq()，也就是curve_fit()使用的那个）。这个包可以做很多事情，其中之一就是可以明确计算置信区间。

Answer 3

注意: 获取拟合曲线置信区间的实际答案可以参考Ulrich在这里的内容。

---

经过一些研究（可以查看这里、这里和1.96），我找到了自己的解决方案。

这个方案可以接受任意的X%置信区间，并绘制出上下的曲线。

以下是最小可重现示例（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()

Answer 4

你可以使用 uncertainties 这个模块来进行不确定性计算。uncertainties 可以帮助你跟踪不确定性和相关性。你可以直接从 curve_fit 的输出中创建相关的 uncertainties.ufloat。

如果你想对一些不是内置的操作，比如 exp，进行这些计算，你需要使用 uncertainties.unumpy 中的函数。

另外，尽量避免使用 from pylab import * 这种导入方式。这会覆盖掉 Python 的一些内置函数，比如 sum。

下面是一个完整的例子：

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)

结果是：

Answer 5

Gabriel的回答是错误的。这里用红色标出了他数据的95%置信带，这是通过GraphPad Prism计算出来的：

背景知识：通常我们把“拟合曲线的置信区间”称为置信带。对于95%的置信带，我们可以有95%的把握它包含了真实的曲线。（这和上面灰色部分的预测带是不同的。预测带是关于未来数据点的。想了解更多，可以查看GraphPad曲线拟合指南的这一页。）

在Python中，kmpfit可以计算非线性最小二乘法的置信带。这里是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()

dfdp是模型f = a*e^(b*x)对每个参数p（也就是a、b和c）的偏导数∂f/∂p。想了解更多背景知识，可以查看kmpfit教程或者GraphPad曲线拟合指南的这一页。（与我的示例代码不同，kmpfit教程没有使用库中的confidence_band()，而是用它自己稍微不同的实现。）

最后，Python绘制的图和Prism的图是匹配的：