具有GNU科学库(GSL)的多项式曲线配件
Polynomial Curve fittings with GNU Scientific Library (GSL)
我将使用GNU科学库(GSL)来求解多项式曲线配件。这是我的polyFit功能 - 请参见" C 代码" 。如果我使用以下示例数据,则我得到以下结果 - 请参见" output" 。我尝试验证是否可以使用Python - 请参阅" Python代码" 和 Python输出。我不知道为什么GSL和Python结果是不同的。Python结果的趋势与原始数据相似。但是,GSL结果是不同的。为什么有所不同?如果有人帮助我,这将非常有帮助。
输出:
> Polynomial Curve Fittings : order = 6 th, Data Size = 200
> best fit: Y = 840810 + 2354.69 X + 1.66468 X^2 + -0.000321755 X^3 + -3.53193e- 007 X^4 + 1.06755e-010 X^5 + -8.11813e-015 X^6
> chisq = 4.76957e+008
C 代码:
struct LensTiltMap{
double xPos;
double yPos1;
double yCompPos1;
};
std::vector<LensTiltMap> polyFit(std::vector<LensTiltMap> _vecData)
{
std::vector<LensTiltMap> _vec;
if (_vecData.size() != 0)
{
double* data_x = &_vecData[0].xPos;
double* data_y = &_vecData[0].yPos1;
int n = _vecData.size();
int order = 6;
double chisq;
gsl_vector *y, *c;
gsl_matrix *X, *cov;
y = gsl_vector_alloc(n);
c = gsl_vector_alloc(order + 1);
X = gsl_matrix_alloc(n, order + 1);
cov = gsl_matrix_alloc(order + 1, order + 1);
for (int i = 0; i < n; i++) {
for (int j = 0; j < order + 1; j++) {
gsl_matrix_set(X, i, j, pow(data_x[i], j));
}
gsl_vector_set(y, i, data_y[i]);
}
gsl_multifit_linear_workspace * work = gsl_multifit_linear_alloc(n, order + 1);
gsl_multifit_linear(X, y, c, cov, &chisq, work);
gsl_multifit_linear_free(work);
std::vector<double> vc;
for (int i = 0; i < order + 1; i++) {
vc.push_back(gsl_vector_get(c, i));
}
cout << "# Polynomial Curve Fittings : order = " << order << " th, Data Size = " << n << endl;
cout << "# best fit: Y = " << vc[0] << " + " << vc[1] << " X + " << vc[2] << " X^2 + " << vc[3] << " X^3 + " << vc[4] << " X^4 + " << vc[5] << " X^5 + " << vc[6] << " X^6 " << endl;
cout << "# covariance matrix = n"
<< "[" << gsl_matrix_get(cov, (0), (0)) << ", " << gsl_matrix_get(cov, (0), (1)) << ", " << gsl_matrix_get(cov, (0), (2)) << ", " << gsl_matrix_get(cov, (0), (3)) << ", " << gsl_matrix_get(cov, (0), (4)) << ", " << gsl_matrix_get(cov, (0), (5)) << ", " << gsl_matrix_get(cov, (0), (6))
<< gsl_matrix_get(cov, (1), (0)) << ", " << gsl_matrix_get(cov, (1), (1)) << ", " << gsl_matrix_get(cov, (1), (2)) << ", " << gsl_matrix_get(cov, (1), (3)) << ", " << gsl_matrix_get(cov, (1), (4)) << ", " << gsl_matrix_get(cov, (1), (5)) << ", " << gsl_matrix_get(cov, (1), (6))
<< gsl_matrix_get(cov, (2), (0)) << ", " << gsl_matrix_get(cov, (2), (1)) << ", " << gsl_matrix_get(cov, (2), (2)) << ", " << gsl_matrix_get(cov, (2), (3)) << ", " << gsl_matrix_get(cov, (2), (4)) << ", " << gsl_matrix_get(cov, (2), (5)) << ", " << gsl_matrix_get(cov, (2), (6))
<< gsl_matrix_get(cov, (3), (0)) << ", " << gsl_matrix_get(cov, (3), (1)) << ", " << gsl_matrix_get(cov, (3), (2)) << ", " << gsl_matrix_get(cov, (3), (3)) << ", " << gsl_matrix_get(cov, (3), (4)) << ", " << gsl_matrix_get(cov, (3), (5)) << ", " << gsl_matrix_get(cov, (3), (6))
<< gsl_matrix_get(cov, (4), (0)) << ", " << gsl_matrix_get(cov, (4), (1)) << ", " << gsl_matrix_get(cov, (4), (2)) << ", " << gsl_matrix_get(cov, (4), (3)) << ", " << gsl_matrix_get(cov, (4), (4)) << ", " << gsl_matrix_get(cov, (4), (5)) << ", " << gsl_matrix_get(cov, (4), (6))
<< gsl_matrix_get(cov, (5), (0)) << ", " << gsl_matrix_get(cov, (5), (1)) << ", " << gsl_matrix_get(cov, (5), (2)) << ", " << gsl_matrix_get(cov, (5), (3)) << ", " << gsl_matrix_get(cov, (5), (4)) << ", " << gsl_matrix_get(cov, (5), (5)) << ", " << gsl_matrix_get(cov, (5), (6))
<< gsl_matrix_get(cov, (6), (0)) << ", " << gsl_matrix_get(cov, (6), (1)) << ", " << gsl_matrix_get(cov, (6), (2)) << ", " << gsl_matrix_get(cov, (6), (3)) << ", " << gsl_matrix_get(cov, (6), (4)) << ", " << gsl_matrix_get(cov, (6), (5)) << ", " << gsl_matrix_get(cov, (6), (6))
<< "]" << endl;
cout << "# chisq = " << chisq << endl;
gsl_vector_free(y);
gsl_vector_free(c);
gsl_matrix_free(X);
gsl_matrix_free(cov);
for (std::vector<LensTiltMap>::iterator it = _vecData.begin(); it != _vecData.end(); ++it) {
LensTiltMap _data = *it;
_data.yCompPos1 = vc[6] + vc[5] * pow(_data.xPos, 1) + vc[4] * pow(_data.xPos, 2) + vc[3] * pow(_data.xPos, 3) + vc[2] * pow(_data.xPos, 4) + vc[1] * pow(_data.xPos, 5) + vc[0] * pow(_data.xPos, 6);
_vec.push_back(_data);
}
}
else
return _vecData;
return _vec.size() != 0 ? _vec : _vecData;
}
Python代码:
请参见带有Python的示例数据。我使用下面的示例链接。
import numpy as np
import matplotlib.pyplot as plt
x = np.arange(-1000,1000,10)
y = np.array([ 5347.21, 5338.78, 5365.01, 5351.12, 5349.49, 5351.44,
5321.54, 5302.74, 5354.44, 5349.04, 5322.55, 5353.69,
5366.55, 5345.69, 5295.52, 5331.35, 5343.48, 5327.36,
5364.93, 5369.18, 5341.57, 5326.26, 5381.95, 5343.6 ,
5372.34, 5341.09, 5341.8 , 5319.17, 5357.89, 5366.52,
5372.47, 5405.77, 5335.64, 5375.94, 5334.32, 5408.44,
5345.63, 5388.27, 5407.22, 5415.23, 5402.14, 5401.65,
5425.57, 5370.68, 5418.62, 5476.2 , 5447.66, 5467.31,
5444.86, 5450.44, 5525.4 , 5489.32, 5494.43, 5457.14,
5504.57, 5555.23, 5520.92, 5513.36, 5585.96, 5621.79,
5558.42, 5608.05, 5596.97, 5599.98, 5583.34, 5610.35,
5679.16, 5666.85, 5695.01, 5693.84, 5722.46, 5726.53,
5714.61, 5722.61, 5733.16, 5699.93, 5753.52, 5754.43,
5745.86, 5828.79, 5772.72, 5825.61, 5819.32, 5852.81,
5876. , 5852.52, 5849.53, 5863.86, 5892.23, 5907.96,
5858.39, 5942.41, 5938.36, 5935.82, 5955.2 , 5910.05,
5958.88, 5995.05, 5923.07, 5968.93, 5933.05, 5920.94,
5930.83, 5993.96, 5919.47, 5956.48, 5948.48, 5966.21,
5990.58, 5996.2 , 5937.79, 5922.37, 5903.46, 5925.97,
5942.13, 5878.51, 5915.93, 5895.85, 5881.16, 5835.25,
5895.39, 5794.58, 5842.72, 5809.81, 5834.05, 5843.11,
5771.03, 5741.2 , 5763.68, 5738.31, 5756.64, 5686.59,
5686.05, 5711.26, 5680.77, 5678. , 5670.78, 5626.55,
5599.49, 5572.86, 5573.88, 5572.26, 5532.51, 5523.21,
5541.77, 5528.95, 5531.11, 5542.49, 5515.9 , 5509.62,
5485.16, 5488.85, 5495.59, 5465.52, 5434.44, 5507.97,
5459.17, 5421.25, 5419.23, 5416.85, 5396.44, 5410.29,
5430.09, 5385.02, 5361.95, 5391.7 , 5345.41, 5350.12,
5345.22, 5370.72, 5322.03, 5348.25, 5370.73, 5338.4 ,
5300.9 , 5325.29, 5323.3 , 5341.07, 5316.03, 5281.7 ,
5333.72, 5287.52, 5355.5 , 5313.96, 5315.16, 5314.75,
5293.81, 5313.89, 5317.65, 5289.2 , 5322.9 , 5275.23,
5273.53, 5278.15, 5291.24, 5260.07, 5290.77, 5272.02,
5284.21, 5317.56])
z=numpy.polyfit(x,y,6)
xp = np.linspace(-1000,1000,200)
p = np.poly1d(z)
_ = plt.plot(x, y, '.', xp, p(xp), '-')
plt.ylim(4000,7000)
plt.show()
Python输出:
z
Out[35]:
array([ -1.93861649e-15, 1.30945729e-13, 3.97750836e-09,
-2.32744951e-07, -2.69634938e-03, 7.39431395e-02,
5.93818062e+03])
http://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html
示例数据:
x = [-1000, -990, -980, -970, -960, -950, -940, -930, -920,
-910, -900, -890, -880, -870, -860, -850, -840, -830,
-820, -810, -800, -790, -780, -770, -760, -750, -740,
-730, -720, -710, -700, -690, -680, -670, -660, -650,
-640, -630, -620, -610, -600, -590, -580, -570, -560,
-550, -540, -530, -520, -510, -500, -490, -480, -470,
-460, -450, -440, -430, -420, -410, -400, -390, -380,
-370, -360, -350, -340, -330, -320, -310, -300, -290,
-280, -270, -260, -250, -240, -230, -220, -210, -200,
-190, -180, -170, -160, -150, -140, -130, -120, -110,
-100, -90, -80, -70, -60, -50, -40, -30, -20,
-10, 0, 10, 20, 30, 40, 50, 60, 70,
80, 90, 100, 110, 120, 130, 140, 150, 160,
170, 180, 190, 200, 210, 220, 230, 240, 250,
260, 270, 280, 290, 300, 310, 320, 330, 340,
350, 360, 370, 380, 390, 400, 410, 420, 430,
440, 450, 460, 470, 480, 490, 500, 510, 520,
530, 540, 550, 560, 570, 580, 590, 600, 610,
620, 630, 640, 650, 660, 670, 680, 690, 700,
710, 720, 730, 740, 750, 760, 770, 780, 790,
800, 810, 820, 830, 840, 850, 860, 870, 880,
890, 900, 910, 920, 930, 940, 950, 960, 970,
980, 990]
y = [ 5347.21, 5338.78, 5365.01, 5351.12, 5349.49, 5351.44,
5321.54, 5302.74, 5354.44, 5349.04, 5322.55, 5353.69,
5366.55, 5345.69, 5295.52, 5331.35, 5343.48, 5327.36,
5364.93, 5369.18, 5341.57, 5326.26, 5381.95, 5343.6 ,
5372.34, 5341.09, 5341.8 , 5319.17, 5357.89, 5366.52,
5372.47, 5405.77, 5335.64, 5375.94, 5334.32, 5408.44,
5345.63, 5388.27, 5407.22, 5415.23, 5402.14, 5401.65,
5425.57, 5370.68, 5418.62, 5476.2 , 5447.66, 5467.31,
5444.86, 5450.44, 5525.4 , 5489.32, 5494.43, 5457.14,
5504.57, 5555.23, 5520.92, 5513.36, 5585.96, 5621.79,
5558.42, 5608.05, 5596.97, 5599.98, 5583.34, 5610.35,
5679.16, 5666.85, 5695.01, 5693.84, 5722.46, 5726.53,
5714.61, 5722.61, 5733.16, 5699.93, 5753.52, 5754.43,
5745.86, 5828.79, 5772.72, 5825.61, 5819.32, 5852.81,
5876. , 5852.52, 5849.53, 5863.86, 5892.23, 5907.96,
5858.39, 5942.41, 5938.36, 5935.82, 5955.2 , 5910.05,
5958.88, 5995.05, 5923.07, 5968.93, 5933.05, 5920.94,
5930.83, 5993.96, 5919.47, 5956.48, 5948.48, 5966.21,
5990.58, 5996.2 , 5937.79, 5922.37, 5903.46, 5925.97,
5942.13, 5878.51, 5915.93, 5895.85, 5881.16, 5835.25,
5895.39, 5794.58, 5842.72, 5809.81, 5834.05, 5843.11,
5771.03, 5741.2 , 5763.68, 5738.31, 5756.64, 5686.59,
5686.05, 5711.26, 5680.77, 5678. , 5670.78, 5626.55,
5599.49, 5572.86, 5573.88, 5572.26, 5532.51, 5523.21,
5541.77, 5528.95, 5531.11, 5542.49, 5515.9 , 5509.62,
5485.16, 5488.85, 5495.59, 5465.52, 5434.44, 5507.97,
5459.17, 5421.25, 5419.23, 5416.85, 5396.44, 5410.29,
5430.09, 5385.02, 5361.95, 5391.7 , 5345.41, 5350.12,
5345.22, 5370.72, 5322.03, 5348.25, 5370.73, 5338.4 ,
5300.9 , 5325.29, 5323.3 , 5341.07, 5316.03, 5281.7 ,
5333.72, 5287.52, 5355.5 , 5313.96, 5315.16, 5314.75,
5293.81, 5313.89, 5317.65, 5289.2 , 5322.9 , 5275.23,
5273.53, 5278.15, 5291.24, 5260.07, 5290.77, 5272.02,
5284.21, 5317.56]
使用正交最小成方拟合algorightm,我得到:
-1.9386E-015 x^6 1.3095e-013 x^5 3.9775e-009 x^4 -2.3274e-007 x^3 -2.6963e-003 x^2 7.3943e-002X 5.9382e 003
与Python输出匹配。链接到我使用的算法的.RTF文档,包括示例C代码:
http://rcgldr.net/misc/opls.rtf
相关文章:
- c++模板来表示多项式
- CGAL:如何创建填充边界曲线的曲面网格?
- 使用链表数据结构打印多项式
- 使用 SFML 绘制曲线会给出错误的形状
- 需要帮助重载多项式类运算符C++
- 使用 GSL 库制作样条曲线并使用它们进行集成
- 如何修复此教科书程序以在C++中添加多项式?
- 如何对两个 4 位数字进行乘法,将它们视为 C++ 中的多项式
- 扩充矩阵的行缩减-三维样条曲线计算
- 实现伪多项式DP子集和
- 向心Catmull-Rom样条曲线插值alpha参数
- 递归程序获得勒让德多项式
- 使用 C++ 在 Qt5 中显示曲线
- 在 NTL 中构造多项式的标准方法是什么?
- 最小二乘多项式拟合仅适用于偶数个坐标
- 多项式系数代码始终提供相同的答案
- 进入第二个多项式后如何修复分割错误?
- 用霍纳方法进行多项式求值的C++ constexpr
- 具有GNU科学库(GSL)的多项式曲线配件
- 图像曲线拟合的多项式最小二乘法