0
我將使用GNU科學圖書館（GSL）解決多項式曲線配件。這裏是我的polyFit函數 - 參見「C++代碼」。如果我使用下面的示例數據，那麼我得到以下結果 - 請參閱「輸出」。我試圖驗證，如果它是好的或不與python - 請參閱「Python代碼」和Python輸出。我不知道爲什麼GSL和Python的結果是不同的。 Python結果的趨勢與原始數據類似。但是，GSL結果不同。爲什麼它不同？如果有人幫助我，這將非常有幫助。GNU科學圖書館（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代碼：
見下文實例數據用蟒。我使用下面的例子鏈接。
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]