Numpy Polynomial 和 Matlab polyfit 生成的线性方程输出之间的差异
Dissimilarity between the output of linear equation produced by Numpy Polynomial and Matlab polyfit
objective就是求一个P(x)次1（一条直线）的多项式的系数最适合数据y的最小二乘法感，然后计算 P.
的根
在 Matlab 中，线性最小二乘拟合可以通过
计算
[p, S, mu] = polyfit(x, y, 1)
这会产生系数 p，它定义了多项式：
-1.5810877
6.0094824
y = p(1) * x + p(2)
mu表示x的mean和std。
在 Python 中，这可以按照建议 :
通过 Numpy 的 polynomial 来实现
import numpy as np
x = [
6210, 6211, 6212, 6213, 6214, 6215, 6216, 6217, 6218, 6219,
6220, 6221, 6222, 6223, 6224, 6225, 6226, 6227, 6228, 6229,
6230, 6231, 6232, 6233, 6234, 6235, 6236, 6237, 6238, 6239,
6240, 6241, 6242, 6243, 6244, 6245, 6246, 6247, 6248, 6249,
6250, 6251, 6252, 6253, 6254, 6255, 6256, 6257, 6258, 6259,
6260, 6261, 6262, 6263, 6264, 6265, 6266, 6267, 6268, 6269,
6270, 6271, 6272, 6273, 6274, 6275, 6276, 6277, 6278, 6279,
6280, 6281, 6282, 6283, 6284, 6285, 6286, 6287, 6288
]
y = [
7.8625913, 7.7713094, 7.6833806, 7.5997391, 7.5211883,
7.4483986, 7.3819046, 7.3221073, 7.2692747, 7.223547,
7.1849418, 7.1533613, 7.1286001, 7.1103559, 7.0982385,
7.0917811, 7.0904517, 7.0936642, 7.100791, 7.1111741,
7.124136, 7.1389918, 7.1550579, 7.1716633, 7.1881566,
7.2039142, 7.218349, 7.2309117, 7.2410989, 7.248455,
7.2525721, 7.2530937, 7.249711, 7.2421637, 7.2302341,
7.213747, 7.1925621, 7.1665707, 7.1356878, 7.0998487,
7.0590014, 7.0131001, 6.9621005, 6.9059525, 6.8445964,
6.7779589, 6.7059474, 6.6284504, 6.5453324, 6.4564347,
6.3615761, 6.2605534, 6.1531439, 6.0391097, 5.9182019,
5.7901659, 5.6547484, 5.5117044, 5.360805, 5.2018456,
5.034656, 4.8591075, 4.6751242, 4.4826899, 4.281858,
4.0727611, 3.8556159, 3.6307325, 3.3985188, 3.1594861,
2.9142516, 2.6635408, 2.4081881, 2.1491354, 1.8874279,
1.6242117,1.3607255,1.0982931,0.83831298
]
p = np.polynomial.polynomial.Polynomial.fit(x, y, 1, domain=[-1, 1])
print(p)
436.53467443432453 - 0.0688950539698132·x¹
mean和std都可以计算如下：
x_mean = np.mean(x)
x_std = np.std(x)
但是我注意到，Matlab 的 polyfit 和 Numpy 的 Polynomial 产生的系数不同：
MATLAB：[-1.5810877, 6.0094824] 对比 Python：[436.53467443432453, 0.0688950539698132]
这个区别对我的用例来说很重要，因为我想计算直线的根。
在 Matlab 中，这可以计算为
>> x_m = p(1);
>> x_mean = mu(1);
>> x_c = p(2);
>> x_std = mu(2);
>> x_intercept = (x_m * x_mean - x_c * x_std) ./ x_m
x_intercept =
6336.2266
而在 Python
>>> x_mean = np.mean(x)
>>> x_std = np.std(x)
>>> x_c, x_m = line.coef
>>> x_intercept = (x_m * x_mean - x_c * x_std) / x_m
>>> x_intercept
150737.19742902054
显然，这两者之间的区别很大。如何在 Python 中复制使用 Matlab 计算的 x_intercept？
为了数值精度，在计算拟合时将多项式数据映射到域 [-1, 1] 或一些类似的 well-bounded 区域很方便。 numpy docs regarding fitting Chebyshev polynomials 中给出了一个经典的例子，如果你缩小得足够远，它似乎会失去所有有趣的特征。该映射对于直线不太有趣，但对于高阶多项式非常有用，其中 x**n 可能会爆炸。特定的 window [-1, 1] 恰好用于防止具有大 n 的任意移位数据的此类爆炸（因为 +/-1**n 始终只是 +/-1，并且 x**n 对于 |x| < 1 总是有界的）。
让我们从查看数据开始（这总是好的第一步）：
我在这里添加了一条直线。仅凭眼睛，我可以看出 ~6340 处的根对于直线拟合是合理的。同样很明显，数据实际上更像立方体，“实际”根在 6290 左右，但这与问题无关。
MATLAB 采用映射 x-data 的方法，使得域的 1-sigma 部分适合 window [-1, 1]。这意味着您获得的多项式系数适用于 (x - mean(xRight)) / std(xRight)。如果您选择 scaling, by requesting three output arguments from polyfit，您将无法选择不同的比例：
>> [p_scaled, s, m] = polyfit(xRight, yRight, 1)
p_scaled =
-1.5811 6.0095
s =
struct with fields:
R: [2×2 double]
df: 77
normr: 8.5045
m =
1.0e+03 *
6.2490
0.0229
>> p_unscaled = polyfit(xRight, yRight, 1)
p_unscaled =
-0.0689 436.5347
您可以使用 polyval:
在 x = 6250 处计算两个拟合值
>> polyval(p_scaled, 6250, s, mu)
ans =
5.9406
>> polyval(p_unscaled, 6250)
ans =
5.9406
>> polyval(p_scaled, (6250 - mean(xRight)) / std(xRight))
ans =
5.9406
当然还有手动：
>> (6250 - m(1)) / m(2) * p_scaled(1) + p_scaled(2)
ans =
5.9406
>> 6250 * p_unscaled(1) + p_unscaled(2)
ans =
5.9406
Python 对 np.polynomial.Polynomial.fit 的 domain 和 window 参数做了类似的事情。正如 MATLAB 将 [-std(x), std(x)] + mean(x) 映射到 [-1, 1] 一样，domain 也映射到 window。最大的不同是domain和window都可以选择。以下是一些常用选项：
>>> p_nomap = np.polynomial.Polynomial.fit(xRight, yRight, 1, domain=[-1, 1]); print(p_nomap)
436.53467443435204 - 0.06889505396981752 x**1
>>> p_default = np.polynomial.Polynomial.fit(xRight, yRight, 1); print(p_default)
6.009482176962028 - 2.686907104822783 x**1
>>> p_minmax = np.polynomial.Polynomial.fit(xRight, yRight, 1, domain=[xRight_arr.min(), xRight_arr.max()]); print(p_minmax)
6.009482176962028 - 2.686907104822783 x**1
>>> p_matlab = np.polynomial.Polynomial.fit(xRight, yRight, 1, domain=np.std(xRight) * np.arange(-1, 2, 2) + np.mean(xRight)); print(p_matlab)
6.009482176962061 - 1.571048948945243 x**1
您可以看到以下等价物：
- Numpy 的
p_matlab <=> MATLAB 的 p_scaled
- Numpy 的
p_nomap <=> MATLAB 的 p_unscaled
- Numpy 的
p_minmax <=> Numpy 的 p_default
您可以使用与 MATLAB 中类似的测试进行验证，除了您可以直接 call a Polynomial 对象而不是 polyval：
>>> p_nomap(6250)
5.940587122992554
>>> p_default(6250)
5.940587122992223
>>> p_minmax(6250)
5.940587122992223
>>> p_matlab(6250)
5.94058712299223
Python 在这种情况下有更好的语法，并且由于您使用的是一个对象，它会在内部为您跟踪所有比例和偏移量。
既然您了解了 domain 的不同选择是如何工作的，那么您可以通过多种方式估算线的根。在 MATLAB 中：
>> roots(p_unscaled)
ans =
6.3362e+03
>> roots(p_scaled) * m(2) + m(1)
ans =
6.3362e+03
>> -p_unscaled(2) / p_unscaled(1)
ans =
6.3362e+03
>> -p_scaled(2) / p_scaled(1) * m(2) + m(1)
ans =
6.3362e+03
你可以在 numpy 中做同样的事情：
>>> p_nomap.roots()
array([6336.22661252])
>>> p_default.roots()
array([6336.22661252])
>>> p_minmax.roots()
array([6336.22661252])
>>> p_matlab.roots()
array([6336.22661252])
手动求值，我们可以使用Polynomial.mapparms，大致相当于MATLAB中的输出参数m。主要区别在于 mapparms 输出 [-np.mean(xRight) / np.std(xRight), 1 / np.std(xRight)]，这简化了评估过程。
>>> -p_nomap.coef[0] / p_nomap.coef[1]
6336.22661251699
>>> (-p_default.coef[0] / p_default.coef[1] - p_default.mapparms()[0]) / p_default.mapparms()[1]
6336.226612516987
>>> (-p_minmax.coef[0] / p_minmax.coef[1] - p_minmax.mapparms()[0]) / p_minmax.mapparms()[1]
6336.226612516987
>>> (-p_matlab.coef[0] / p_matlab.coef[1] - p_matlab.mapparms()[0]) / p_matlab.mapparms()[1]
6336.226612516987
因此，获得等效结果的关键是要么选择相同的参数（并理解它们的对应关系），要么使用提供的函数来计算多项式及其根。总的来说我推荐后者而不考虑前者。提供这些功能的原因是您可以从 domain 的任何选择中获得准确的结果，只要您始终如一地传递数据（对于 polyval 手动传递，对于 Polynomial.__call__ 自动传递）
objective就是求一个P(x)次1（一条直线）的多项式的系数最适合数据y的最小二乘法感，然后计算 P.
在 Matlab 中，线性最小二乘拟合可以通过
计算[p, S, mu] = polyfit(x, y, 1)
这会产生系数 p，它定义了多项式：
-1.5810877
6.0094824
y = p(1) * x + p(2)
mu表示x的mean和std。
在 Python 中，这可以按照建议
polynomial 来实现
import numpy as np
x = [
6210, 6211, 6212, 6213, 6214, 6215, 6216, 6217, 6218, 6219,
6220, 6221, 6222, 6223, 6224, 6225, 6226, 6227, 6228, 6229,
6230, 6231, 6232, 6233, 6234, 6235, 6236, 6237, 6238, 6239,
6240, 6241, 6242, 6243, 6244, 6245, 6246, 6247, 6248, 6249,
6250, 6251, 6252, 6253, 6254, 6255, 6256, 6257, 6258, 6259,
6260, 6261, 6262, 6263, 6264, 6265, 6266, 6267, 6268, 6269,
6270, 6271, 6272, 6273, 6274, 6275, 6276, 6277, 6278, 6279,
6280, 6281, 6282, 6283, 6284, 6285, 6286, 6287, 6288
]
y = [
7.8625913, 7.7713094, 7.6833806, 7.5997391, 7.5211883,
7.4483986, 7.3819046, 7.3221073, 7.2692747, 7.223547,
7.1849418, 7.1533613, 7.1286001, 7.1103559, 7.0982385,
7.0917811, 7.0904517, 7.0936642, 7.100791, 7.1111741,
7.124136, 7.1389918, 7.1550579, 7.1716633, 7.1881566,
7.2039142, 7.218349, 7.2309117, 7.2410989, 7.248455,
7.2525721, 7.2530937, 7.249711, 7.2421637, 7.2302341,
7.213747, 7.1925621, 7.1665707, 7.1356878, 7.0998487,
7.0590014, 7.0131001, 6.9621005, 6.9059525, 6.8445964,
6.7779589, 6.7059474, 6.6284504, 6.5453324, 6.4564347,
6.3615761, 6.2605534, 6.1531439, 6.0391097, 5.9182019,
5.7901659, 5.6547484, 5.5117044, 5.360805, 5.2018456,
5.034656, 4.8591075, 4.6751242, 4.4826899, 4.281858,
4.0727611, 3.8556159, 3.6307325, 3.3985188, 3.1594861,
2.9142516, 2.6635408, 2.4081881, 2.1491354, 1.8874279,
1.6242117,1.3607255,1.0982931,0.83831298
]
p = np.polynomial.polynomial.Polynomial.fit(x, y, 1, domain=[-1, 1])
print(p)
436.53467443432453 - 0.0688950539698132·x¹
mean和std都可以计算如下：
x_mean = np.mean(x)
x_std = np.std(x)
但是我注意到，Matlab 的 polyfit 和 Numpy 的 Polynomial 产生的系数不同：
MATLAB：[-1.5810877, 6.0094824] 对比 Python：[436.53467443432453, 0.0688950539698132]
这个区别对我的用例来说很重要，因为我想计算直线的根。
在 Matlab 中，这可以计算为
>> x_m = p(1);
>> x_mean = mu(1);
>> x_c = p(2);
>> x_std = mu(2);
>> x_intercept = (x_m * x_mean - x_c * x_std) ./ x_m
x_intercept =
6336.2266
而在 Python
>>> x_mean = np.mean(x)
>>> x_std = np.std(x)
>>> x_c, x_m = line.coef
>>> x_intercept = (x_m * x_mean - x_c * x_std) / x_m
>>> x_intercept
150737.19742902054
显然，这两者之间的区别很大。如何在 Python 中复制使用 Matlab 计算的 x_intercept？
为了数值精度，在计算拟合时将多项式数据映射到域 [-1, 1] 或一些类似的 well-bounded 区域很方便。 numpy docs regarding fitting Chebyshev polynomials 中给出了一个经典的例子，如果你缩小得足够远，它似乎会失去所有有趣的特征。该映射对于直线不太有趣，但对于高阶多项式非常有用，其中 x**n 可能会爆炸。特定的 window [-1, 1] 恰好用于防止具有大 n 的任意移位数据的此类爆炸（因为 +/-1**n 始终只是 +/-1，并且 x**n 对于 |x| < 1 总是有界的）。
让我们从查看数据开始（这总是好的第一步）：
我在这里添加了一条直线。仅凭眼睛，我可以看出 ~6340 处的根对于直线拟合是合理的。同样很明显，数据实际上更像立方体，“实际”根在 6290 左右，但这与问题无关。
MATLAB 采用映射 x-data 的方法，使得域的 1-sigma 部分适合 window [-1, 1]。这意味着您获得的多项式系数适用于 (x - mean(xRight)) / std(xRight)。如果您选择 scaling, by requesting three output arguments from polyfit，您将无法选择不同的比例：
>> [p_scaled, s, m] = polyfit(xRight, yRight, 1)
p_scaled =
-1.5811 6.0095
s =
struct with fields:
R: [2×2 double]
df: 77
normr: 8.5045
m =
1.0e+03 *
6.2490
0.0229
>> p_unscaled = polyfit(xRight, yRight, 1)
p_unscaled =
-0.0689 436.5347
您可以使用 polyval:
x = 6250 处计算两个拟合值
>> polyval(p_scaled, 6250, s, mu)
ans =
5.9406
>> polyval(p_unscaled, 6250)
ans =
5.9406
>> polyval(p_scaled, (6250 - mean(xRight)) / std(xRight))
ans =
5.9406
当然还有手动：
>> (6250 - m(1)) / m(2) * p_scaled(1) + p_scaled(2)
ans =
5.9406
>> 6250 * p_unscaled(1) + p_unscaled(2)
ans =
5.9406
Python 对 np.polynomial.Polynomial.fit 的 domain 和 window 参数做了类似的事情。正如 MATLAB 将 [-std(x), std(x)] + mean(x) 映射到 [-1, 1] 一样，domain 也映射到 window。最大的不同是domain和window都可以选择。以下是一些常用选项：
>>> p_nomap = np.polynomial.Polynomial.fit(xRight, yRight, 1, domain=[-1, 1]); print(p_nomap)
436.53467443435204 - 0.06889505396981752 x**1
>>> p_default = np.polynomial.Polynomial.fit(xRight, yRight, 1); print(p_default)
6.009482176962028 - 2.686907104822783 x**1
>>> p_minmax = np.polynomial.Polynomial.fit(xRight, yRight, 1, domain=[xRight_arr.min(), xRight_arr.max()]); print(p_minmax)
6.009482176962028 - 2.686907104822783 x**1
>>> p_matlab = np.polynomial.Polynomial.fit(xRight, yRight, 1, domain=np.std(xRight) * np.arange(-1, 2, 2) + np.mean(xRight)); print(p_matlab)
6.009482176962061 - 1.571048948945243 x**1
您可以看到以下等价物：
- Numpy 的
p_matlab<=> MATLAB 的p_scaled - Numpy 的
p_nomap<=> MATLAB 的p_unscaled - Numpy 的
p_minmax<=> Numpy 的p_default
您可以使用与 MATLAB 中类似的测试进行验证，除了您可以直接 call a Polynomial 对象而不是 polyval：
>>> p_nomap(6250)
5.940587122992554
>>> p_default(6250)
5.940587122992223
>>> p_minmax(6250)
5.940587122992223
>>> p_matlab(6250)
5.94058712299223
Python 在这种情况下有更好的语法，并且由于您使用的是一个对象，它会在内部为您跟踪所有比例和偏移量。
既然您了解了 domain 的不同选择是如何工作的，那么您可以通过多种方式估算线的根。在 MATLAB 中：
>> roots(p_unscaled)
ans =
6.3362e+03
>> roots(p_scaled) * m(2) + m(1)
ans =
6.3362e+03
>> -p_unscaled(2) / p_unscaled(1)
ans =
6.3362e+03
>> -p_scaled(2) / p_scaled(1) * m(2) + m(1)
ans =
6.3362e+03
你可以在 numpy 中做同样的事情：
>>> p_nomap.roots()
array([6336.22661252])
>>> p_default.roots()
array([6336.22661252])
>>> p_minmax.roots()
array([6336.22661252])
>>> p_matlab.roots()
array([6336.22661252])
手动求值，我们可以使用Polynomial.mapparms，大致相当于MATLAB中的输出参数m。主要区别在于 mapparms 输出 [-np.mean(xRight) / np.std(xRight), 1 / np.std(xRight)]，这简化了评估过程。
>>> -p_nomap.coef[0] / p_nomap.coef[1]
6336.22661251699
>>> (-p_default.coef[0] / p_default.coef[1] - p_default.mapparms()[0]) / p_default.mapparms()[1]
6336.226612516987
>>> (-p_minmax.coef[0] / p_minmax.coef[1] - p_minmax.mapparms()[0]) / p_minmax.mapparms()[1]
6336.226612516987
>>> (-p_matlab.coef[0] / p_matlab.coef[1] - p_matlab.mapparms()[0]) / p_matlab.mapparms()[1]
6336.226612516987
因此，获得等效结果的关键是要么选择相同的参数（并理解它们的对应关系），要么使用提供的函数来计算多项式及其根。总的来说我推荐后者而不考虑前者。提供这些功能的原因是您可以从 domain 的任何选择中获得准确的结果，只要您始终如一地传递数据（对于 polyval 手动传递，对于 Polynomial.__call__ 自动传递）