python - 指数曲线拟合的置信区间

Question

我正在尝试获取某些 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%(或其他值)置信区间>(我已经安装了哪些软件包)？

Answer 1

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

为了能够对非内置操作(例如 exp)进行这些计算，您需要使用 uncertainties.unumpy 中的函数。

您还应该避免使用 from pylab import * 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)

结果: