如何使用 kapteyn.kmpfit 计算具有 2 个或更多独立变量的模型的置信带

Question

我尝试计算具有两个自变量的模型的置信带。

我使用 kapteyn 包中的 kmpfit 来优化参数和计算置信区间。 confidence_band 函数似乎不适用于具有多个独立参数的模型。有谁知道如何解决这个问题？

对于具有一个独立参数的模型，一切都按预期工作，没有任何问题：

from kapteyn import kmpfit
import numpy as np
import matplotlib.pyplot as plt 


#linear 1d model
def model_1d(p, x):
    a = p[0]
    b = p[1]
    return a*x + b


#data
x = np.array([0, 10, 20, 30])
y  = np.array([1.1, 1.8, 3.2, 3.0])

#optimizing parameters a and b with kmpfit
p0 = [0, 0] #start values for parameters a and b
f = kmpfit.simplefit(model_1d, p0, x, y)    


#calculation of confidence band

derivatives_a = x #1st derivative of a is x
derivatives_b = np.array([1, 1, 1, 1]) #1st derivative of b is always 1
dfdp = [derivatives_a, derivatives_b]

yhat, upper, lower = f.confidence_band(x, dfdp, 0.95, model_1d)


#plots
plt.plot(x,yhat) #plot fitted model
plt.plot(x,upper,'--') #plot upper confidence band
plt.plot(x,lower,'--') #plot lower confidence band
plt.plot(x,y,'.r') #plot data points

脚本运行良好。因此，您会看到具有上下置信带的拟合模型： plot with upper and lower confidence bands

下面代码稍微改成二维模型，即有两个自变量。此代码不再有效：

from kapteyn import kmpfit
import numpy as np

#linear with two independent variables x[0] amd x[1]
def model_2d(p, x):
    a = p[0]
    b = p[1]
    return a*x[0] + b*x[1]


#data
x = np.array([[0, 10, 20, 30], [0, 10, 20, 30]])
y  = [1.1, 1.8, 3.2, 3.0]

#optimizing parameters a and b with kmpfit
p0 = [0, 0] #start values for parameters a and b
f = kmpfit.simplefit(model_2d, p0, x, y)    


#calculation of confidence band

derivatives_a = x[0,:] #1st derivative of a is x[0]
derivatives_b = x[1,:] #1st derivative of b is x[1]
dfdp = [derivatives_a, derivatives_b]

yhat, upper, lower = f.confidence_band(x, dfdp, 0.95, model_2d)

对于 2D 模型，使用 simplefit 进行参数优化仍然有效。但是，使用 confidence_band 计算置信区间不再有效。我收到以下错误消息，表明 numpy 数组 x 的形状会产生问题（但是，相同的 x 适用于 simplefit）：

File "tmp.py", line 25, in <module>
    yhat, upper, lower = f.confidence_band(x, dfdp, 0.95, model_2d)
  File "src/kmpfit.pyx", line 845, in kmpfit.Fitter.confidence_band 
(src/kmpfit.c:8343)

如果您知道如何解决这个问题，或者如果有可用的替代 python 软件包，我将不胜感激。谢谢！

Answer 1

您可能会发现 lmfit (http://lmfit.github.io/lmfit-py/) 对此很有用。它的模型 class 具有强大但易于使用的曲线拟合方法，并且还支持对每个参数设置边界和约束。

一般来说，所有变量的不确定性和变量对之间的相关性是使用快速且通常可靠的曲率矩阵求逆方法来计算的。对于更详细的不确定性（包括任何指定 'sigma' 级别的置信度），它还允许直接（如果稍微慢一点）探索任何变量的置信度级别，并且可以帮助您制作卡方的二维图一对的正方形变量。

更具体地说，我认为您要求的是根据变量中的不确定性生成模型 y 中的不确定性。 Lmfit 的模型 class 有一个 eval_uncertainty() 方法来执行此操作。看 http://lmfit.github.io/lmfit-py/model.html#calculating-uncertainties-in-the-model-function 举个例子。

如何使用 kapteyn.kmpfit 计算具有 2 个或更多独立变量的模型的置信带

How to calculate confidence bands for models with 2 or more independent variables with kapteyn.kmpfit

python

scipy

curve-fitting