Scipy.curve_fit() 与 Matlab fit() 加权非线性最小二乘法

Question

我有一个 Matlab 参考例程，我正试图将其转换为 numpy/scipy。我遇到了无法在 Python 中解决的曲线拟合问题。所以这里有一个简单的例子来说明这个问题。数据是完全合成的，不是问题的一部分。

假设我正在尝试拟合噪声数据的直线模型 -

x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
y = [0.1075, 1.3668, 1.5482, 3.1724, 4.0638, 4.7385, 5.9133, 7.0685, 8.7157, 9.5539]

对于 Matlab 中的未加权解决方案，我将编码

g = @(m, b, x)(m*x + b)
f = fittype(g)
bestfit = fit(x, y, g)

生成 bestfit.m = 1.048、bestfit.b = -0.09219

的解

运行此数据通过 scipy.optimize.curve_fit() 产生相同的结果。

如果相反，拟合使用衰减函数来减少数据点的影响

dw = [0.7290, 0.5120, 0.3430, 0.2160, 0.1250, 0.0640, 0.0270, 0.0080,  0.0010, 0]
weightedfit = fit(x, y, g, 'Weights', dw)

这会产生斜率为 0.944 且偏移为 0.1484。

我还没有想出如何使用 sigma 参数从 scipy.optimize.curve_fit 中得出这个结果。如果我按照提供给 Matlab 的方式传递权重，“0”会导致除以零异常。显然，Matlab 和 scipy 对底层优化程序中权重的含义的思考非常不同。是否有一种简单的方法可以在两者之间进行转换，使我能够提供产生相同结果的加权函数？

Answer 1

好的，经过进一步调查，我可以提供答案，至少对于这个简单的例子。

import numpy as np
import scipy as sp
import scipy.optimize

def modelFun(x, m, b):
    return m * x + b

def testFit():
    w = np.diag([1.0, 1/0.7290, 1/0.5120, 1/0.3430, 1/0.2160, 1/0.1250, 1/0.0640, 1/0.0270, 1/0.0080, 1/0.0010])
    x = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
    y = np.array([0.1075, 1.3668, 1.5482, 3.1724, 4.0638, 4.7385, 5.9133, 7.0685, 8.7157, 9.5539])

    popt = sp.optimize.curve_fit(modelFun, x, y, sigma=w)

    print(popt[0])
    print(popt[1])

这会产生所需的结果。

为了强制sp.optimize.curve_fit使用曲线拟合工具箱最小化与 Matlab 相同的 chisq 度量，您必须做两件事：

使用权重因子的倒数
根据新的权重因子创建对角矩阵。根据 scipy 参考：

sigma None or M-length sequence or MxM array, optional Determines the uncertainty in ydata. If we define residuals as r = ydata - f(xdata, *popt), then the interpretation of sigma depends on its number of dimensions:

A 1-d sigma should contain values of standard deviations of errors in ydata. In this case, the optimized function is chisq = sum((r / sigma) ** 2).

A 2-d sigma should contain the covariance matrix of errors in ydata. In this case, the optimized function is chisq = r.T @ inv(sigma) @ r.

New in version 0.19.

None (default) is equivalent of 1-d sigma filled with ones.

Scipy.curve_fit() 与 Matlab fit() 加权非线性最小二乘法

Scipy.curve_fit() vs. Matlab fit() weighted nonlinear least squares

matlab

scipy-optimize