logo标签列表

指数曲线拟合的置信区间

pythonnumpyscipyconfidence-interval
17

我正在尝试对一些数据（可以在此处获取）进行指数拟合的置信区间。以下是我使用的最小工作示例，以找到最佳指数拟合：

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包？

https://dev59.com/j14d5IYBdhLWcg3wG_aH#63560689 - Marco Cerliani
5个回答
11
您可以使用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)

以下为结果： result

我之前不知道uncertainties这个包，看起来非常有趣，我会试一试。非常感谢！ - Gabriel
哦，天啊，我甚至不知道这些年来我一直在寻找它。 - Moritz
7
加布里埃尔的答案（链接）是错误的。下面的红色部分是GraphPad Prism计算出的他数据的95%置信区间： Prism confidence and prediction bands 背景：拟合曲线的“置信区间”通常称为置信区带。对于95%置信区带，我们可以有95%的把握认为它包含了真实曲线。（这和上面显示的灰色预测区带不同。预测区带是关于未来数据点的。更多细节请参见GraphPad Curve Fitting Guide的这个页面。）
在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。有关背景，请参见kmpfit教程或GraphPad Curve Fitting Guide的此页面。(与我的示例代码不同，kmpfit教程不使用库中的confidence_band()，而是使用稍微不同的自己的实现。)
最后，Python图表与Prism图表匹配:

kmpfit confidence bands

太棒了，Ulrich，非常感谢你的回答！事实上，我认为我的旧答案实际上获得了预测带，而不是拟合曲线的置信区间。你似乎对这些统计数据很熟悉，能确认一下吗？ - Gabriel
3
我刚刚在Prism图中添加了预测带。因此，您之前的答案没有计算预测带。GraphPad曲线拟合指南的这个页面说明了Prism如何计算预测带。 - Ulrich Stern
6

注意: 获取拟合曲线置信区间的实际答案由Ulrich在这里给出。

经过一些研究（请参见这里，这里和1.96），我想到了自己的解决方案。

它可以接受任意X%的置信区间并绘制上下曲线。

enter image description here

以下是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或其他名称，以免让人们感到困惑。我可以理解这是为置信线而设定的。 - Srivatsan
我还发现另一种解决方案。curve_fit 中的协方差矩阵 pcov 也可以用作1sigma误差。请查看该网站。 - Srivatsan
@ThePredator 因为如果我调用完整的 x 而不是 linspace(11,23,100)，函数将尝试绘制拟合曲线以触及所有这些 x 值。 你可以自己试试，将 x = linspace(11, 23, 100) 注释掉，看看会发生什么 :) - Gabriel
@ThePredator 协方差矩阵pcov正是我答案所使用的： perr = np.sqrt(np.diag(pcov))。那些是1 sigma误差，这就是获取它们的方法（请参见此处）。 - Gabriel
我认为这个解决方案不正确。我看到两个主要问题：（1）选择一个参数置信区间的边际将使您达到95％，同时也考虑第二个参数，您将得到1-0.05 ** 2-> 99.75％。因此，您的置信区间要大得多。 （2）您假设您的参数是独立的，这只是在协方差很小的情况下才是合法的近似值。 - Suuuehgi
4

curve_fit()会返回协方差矩阵 - pcov，其中包含了估计的不确定性（1 sigma）。这假设误差是正态分布的，但有时是值得怀疑的。

您还可以考虑使用lmfit包（纯Python，基于scipy构建），它提供了对scipy.optimize拟合例程的包装器（包括curve_fit()使用的leastsq()）等功能，可以明确地计算置信区间。

1
我一直是简单引导获取置信区间的粉丝。如果您有n个数据点，则使用random包从您的数据中选择n个点进行重采样（即允许程序多次获取相同的点，如果程序想这么做的话-非常重要）。完成后，绘制重采样点并获得最佳拟合线。重复此过程10,000次，每次获得一个新的拟合线。然后，您的95％置信区间是包围95％最佳拟合线的线对。
在Python中编程这种方法相当容易，但从统计学角度来看，这将如何运作还不太清楚。一些关于为什么要这样做的更多信息可能会导致更适合您任务的答案。
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接