指数曲线拟合的置信区间

Posted 2023-03-12

【中文标题】指数曲线拟合的置信区间

【英文标题】：Confidence interval for exponential curve fit

【发布时间】：2014-08-29 07:08:14

【问题描述】：

我正在尝试获得与某些x,y 数据（可用here）的指数拟合的置信区间。这是我必须找到最适合数据的指数拟合的 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()

产生：

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

【问题讨论】：

***.com/a/63560689/10375049

【参考方案1】：

curve_fit() 返回协方差矩阵 - pcov - 保存估计的不确定性 (1 sigma)。这假设错误是正态分布的，这有时是有问题的。

您也可以考虑使用lmfit 包（纯python，建立在scipy 之上），它提供了一个围绕scipy.optimize 拟合例程的包装器（包括 minimumsq()，这是curve_fit() 使用的）并且可以，除其他外，明确计算置信区间。

【参考方案2】：

您可以使用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)

结果：

【讨论】：

我不知道uncertainties 包，我会试一试，它看起来很有趣。非常感谢！

天哪，我什至不知道我多年来一直在寻找它。

【参考方案3】：

Gabriel 的answer 不正确。图中红色表示他的数据的 95% 置信带，由 Graphpad Prism 计算得出：

背景：“拟合曲线的置信区间”通常称为置信带。对于 95% 的置信带，可以有 95% 的置信度认为它包含真实曲线。 （这不同于 预测带，上面以灰色显示。预测带是关于未来数据点的。有关更多详细信息，请参阅 GraphPad 曲线拟合指南的 page。）

在 Python 中，kmpfit 可以计算非线性最小二乘的置信区间。这里是加布里埃尔的例子：

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) + c 关于每个参数 p（即 a、b 和 c）的偏导数 ∂f/∂p。有关背景，请参阅 GraphPad 曲线拟合指南的 kmpfit Tutorial 或 page。 （与我的示例代码不同，kmpfit 教程不使用库中的 confidence_band()，而是使用它自己的、略有不同的实现。）

最后，Python 图与 Prism 相匹配：

【讨论】：

很好的回答乌尔里希，非常感谢！事实上，我相信我的旧答案实际上获得了预测带，而不是拟合曲线的置信区间。您似乎对这些统计数据了如指掌，您能确认一下吗？

我刚刚在 Prism 图中添加了预测波段。因此，您的旧答案不会计算预测范围。 GraphPad 曲线拟合指南的 page 说明了它们是如何在 Prism 中计算的。

【参考方案4】：

注意：获得拟合曲线的置信区间的实际答案由 Ulrich here 给出。

---

经过一些研究（参见here、here 和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()

【讨论】：

当你已经拥有来自exponential_data.dat 的x,y 为什么还要说x = linspace(11,23,100)。最好将其称为X1 或其他名称，以免人们感到困惑。我可以理解这是为了信心线。

我还找到了另一种解决方案。我们的curve_fit 中的协方差矩阵pcov 具有1sigma 误差，也可以使用。检查这个WEBSITE

@ThePredator 因为如果我调用完整的x 而不是linspace(11,23,100)，该函数将尝试绘制拟合曲线以触及所有x 值。自己尝试一下，注释掉x = linspace(11, 23, 100)，看看会发生什么:)

@ThePredator 协方差矩阵pcov 正是我的答案使用的：perr = np.sqrt(np.diag(pcov))。这些是 1 sigma 错误，这就是您获取它们的方式（参见here）

我不认为这个解决方案是正确的。我在这里看到两个主要问题：（1）选择一个参数置信区间的边际可以达到 95%，而第二个参数的置信区间也可以达到 1-0.05**2 --> 99.75%。所以你在这里的置信区间要大得多。 (2) 你假设你的参数是独立的，只有当你的协方差很小时，什么才是合法的近似。

【参考方案5】：

我一直喜欢通过简单的引导来获得置信区间。如果您有n 数据点，则使用random 包从您的数据中选择n 点并重新采样（即，如果您的程序想要这样做，则允许您的程序多次获得相同的点 - 非常重要） .完成后，绘制重采样点并获得最佳拟合。这样做 10,000 次，每次都获得一条新的合身线。那么您的 95% 置信区间是包含您制作的 95% 的最佳拟合线的那对线。

这是一种在 Python 中编程的非常简单的方法，但从统计的角度来看，它是如何实现的有点不清楚。有关您为什么要这样做的更多信息可能会为您的任务带来更合适的答案。