scipy.optimze curve_fit() 给出了高标准误差，即使拟合看起来不错

Question

所以我试图将线性函数拟合到磁滞曲线上。首先，我使用破裂库来查找变化点以确定线性区域。然后我使用

curve_fit()

来拟合该区域的线性方程。

#Ausgleichsgeraden berrechnen
def fit(x,a,b):
    return a*x+b

#Messung 1
popt_x11, pcov_x11  = curve_fit(fit, Magnet_x11, Kerr_x11)
popt_x21, pcov_x21  = curve_fit(fit, Magnet_x21, Kerr_x21)
popt_y11, pcov_y11  = curve_fit(fit, Magnet_y11, Kerr_y11)
popt_y21, pcov_y21  = curve_fit(fit, Magnet_y21, Kerr_y21)

Magnet_x..

和

Kerr_x..

是曲线的线性区域。 但

np.sqrt(np.diag(pcov))

产生的标准误差至少比测量值高出 10 倍，这似乎不适合，因为绘制的线性曲线似乎非常适合数据。现在我不知道如何解决这个问题。

有人知道高偏差从何而来吗？

这是我尝试拟合线性区域的示例曲线：

这是我用来确定每个线性区域范围的破裂图：

例如，对于所示的曲线，我得到了

Messung 2: HC_rechts = 4693.0262 +- 36583.3256

的错误和

Messung 2: HC_links = -3364.4523 +- 40484.9339

。我认为这对于合身质量来说太高了。

我尝试更改拟合曲线的范围，并且还尝试使用

scipy.optimize

查找错误，但没有成功。

我读到高标准偏差来自

curve_fit

中未使用的参数，但这里不应该是这样。

This是我正在绘制并尝试拟合的示例数据。在 X 轴上，我绘制了第一列，其中没有任何不确定性。在 Y 轴上，我绘制了第三列除以第二列的两倍的图。

Answer 1

关于协方差

协方差是使用梯度下降返回的解处的数据和模型之间的误差的 Hessian 的倒数来估计的。

当协方差较高时，意味着解处的 Hessian 矩阵较低。因此，误差函数的曲率很低，梯度下降很慢。

在以下情况下，不确定性可能比参数高几个数量级：

1. 找到的解决方案不是最佳的，但您落在了错误景观的平坦区域；
2. 测量误差非常高，模型被噪声淹没。

关于你的价值观

我无法使用您提供的数据集重现您的参数和不确定性。

如果没有您使用的确切代码，就不可能说出原因，但您可以说服自己，至少参数是不可能的，因为您的数据的斜率 4e3 无法跨越 4e-3，所以您很可能处于情况 1 中。上面。

MCVE

我在这里留下了一个 MCVE，展示了如何回归你的参数：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import optimize
from sklearn.metrics import r2_score

def model(x, a, b):
    return a * x + b

加载数据：

data = pd.read_csv("Hysteresis.txt", sep="\t", header=None, names=["x", "y", "z"])
data["u"] = data["z"] / (2. * data["y"])

标签数据：

data["t"] = 0
data.loc[488:530, "t"] = 1
data.loc[530:1390, "t"] = 2
data.loc[1390:1430, "t"] = 3

分别回归每个标签：

fig, axe = plt.subplots()
for i in data["t"].unique():
    
    q = data["t"] == i
    x = data.loc[q, "x"]
    y = data.loc[q, "u"]
    
    popt, pcov = optimize.curve_fit(model, x, y)
    yhat = model(x, *popt)
    score = r2_score(y, yhat)
    
    print(i, score, popt, np.sqrt(np.diag(pcov)))
    
    xlin = np.linspace(x.min(), x.max(), 100)
    ylin = model(xlin, *popt)
    
    axe.scatter(x, y, marker=".", s=4)
    axe.plot(xlin, ylin)
    
axe.grid()

结果稳定（即使没有标准化）：

# 0 0.844119238551804   [-4.29585719e-07 -1.27397584e-04] [6.27676896e-09 1.52481902e-06]
# 1 0.733147409548466   [ 8.58502916e-06 -3.05191197e-04] [8.18939540e-07 5.56526052e-05]
# 2 0.43378140123287146 [-1.62503164e-07  4.12090730e-04] [6.33832903e-09 1.69385354e-06]
# 3 0.6562435542810321  [1.18757579e-05 5.22226928e-04]   [1.37632797e-06 4.97901997e-05]

如您所见，斜率约为 1e-6 而不是 1e4，不确定性很高，但对于您的数据集来说是可以接受的。