i'm going use gnu scientific library (gsl) solving polynomial curve fittings. here function polyfit - see "c++ code". if use below example data, got result below - see "output". i've tried verification if ok or not python - see "python code" , python output. don't know why gsl , python result different. trend of python result similar original data. however, gsl result different. why different ? helpful if me.
output :
> 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++ code :
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); (int = 0; < n; i++) { (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; (int = 0; < 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); (std::vector<lenstiltmap>::iterator = _vecdata.begin(); != _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 code :
see below example data python. use example below link.
import numpy np import matplotlib.pyplot 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 output :
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
example data :
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]
using orthogonal least square fit algorightm, get:
-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-002 x + 5.9382e+003
which matches python output. link .rtf document algorithm use, including example c code:
Comments
Post a Comment