C/C++ code with GNU Scientific Library (GSL) gives different results to GNUPlot - possible floating point instabilities? -

sshort:

gnuplot gives better fit data gsl code does. why?

short:

i confused @ moment, question might not particularly worded... edit understanding improves.

the original title of question was: "g++ compiling code either -o1 -o2 or -o3 flags , floating point precision"

i believe code suffering numerical instabilities.

gnuplot gives better fit data gsl code does, surprising since believe gnuplot uses gsl libraries?

long:

i have written c / c++ code makes use of gnu scientific library (gsl). code fits non-linear functions sets of non-linear data. algorithms can highly sensitive order in floating point operations take place, due nature of numerical inaccuracies result accumulation of numerical round-off errors.

question: "could these caused effect of running 1 of optimization flags, -o1, -o2 or -o3?"

partial answer: turned off -on flags , recompiled code, results may have changed small amount, ie: delta_x / x ~= 1.0e-3. fit still poor in comparison gnuplot.

functions fitting:

i have provided these show numerical work occurring. suspect of prone numerical error.

typical values yi in range of 0.0 1.0. t typically in range 0.0 200.0. (but fit poor in first half of range.)

// block - weighted residuals  double t = time; // whatever may double s = sigma[ix]; // these errors, tried 1.0 , 0.001 no affect on parameter values obtained  double yi = (std::sqrt(rho) * r0) / std::sqrt((rho - r0*r0) * std::exp(-2.0 * t / tau) + r0*r0); // y value - model  gsl_vector_set(f, ix, (yi - y[ix])/sigma[ix]); // weighted residual  // block b - jacobian  double mem_a = std::exp(-2.0 * t / tau); // tried optimization here double mem_b = 1.0 - mem_a; // perhaps causing problem? double mem_c = rho - r0 * r0; double mem_d = mem_c * mem_a + r0*r0; double mem_e = std::pow(mem_d, 1.5);  double df_drho = (std::pow(r0, 3.0) * mem_b) / (2.0 * std::sqrt(rho) * mem_e); double df_dr0 = (std::pow(rho, 1.5) * mem_a) / mem_e; double df_dtau = -1.0 * (std::sqrt(rho) * r0 * mem_c * mem_a * t) / (tau * tau * mem_e);  gsl_matrix_set(j, ix, 0, df_drho / s); gsl_matrix_set(j, ix, 1, df_dr0 / s); gsl_matrix_set(j, ix, 2, df_dtau / s); 

here graph, isn't nice?

well here graph explains issue lot better have been able in words. may ignore green line, shows initial parameters given before fitting algorithm run changes parameters.

gnuplot fit results:

rhofit = 0.086173236829715 +- 2.61304934752193e-05  r0fit = 0.00395856812689133 +- 2.08898744280108e-05 taufit = 11.7694359189233 +- 0.016094629240588  // not sure how gnuplot calculates errors - not appear regular errors computed off diagonal elements of lm matrix after fitting. (if know little levenberg–marquardt algorithm.) 

c++ gsl fit results:

rho    = 0.08551510 +/- ... r0     = 0.00507645 +/- ... // nonsense errors due "not-real" errors on data points tau    = 12.99114719 +/- ... 

on careful inspection, see pink , blue lines not overlay each other quite considerable margin. pink line many describe "good fit". blue line comparison isn't particularly good.

i have tried making errorbars (although same size points - aren't "real" errorbars, artificial ones) smaller - doesn't help, changes chi-square value , associated errors of each parameter after fitting.

further random ideas:

• i built gsl wrong?
• gnuplot splits dataset small blocks keep numbers added same order of magnitude? (kind of how fft works.)

gsl fit output:

iter: 0 x = 0.1 0.001 10 |f(x)| = 12487.8 status = success iter: 1 x = 0.0854247 0.00323946 13.2064 |f(x)| = 10476.9 dx vector: -0.0145753, 0.00223946, 3.20642 status = success iter: 2 x = 0.0854309 0.00576809 13.7443 |f(x)| = 3670.4 dx vector: 6.18836e-06, 0.00252863, 0.537829 chisq/dof = 6746.03 rho    = 0.08543089 +/- 0.00013518 r0     = 0.00576809 +/- 0.00013165 tau    = 13.74425294 +/- 0.09012196 

i came across page because having exact same problem. needed fit function gsl, hadn't done before comparing results gnuplot's fitting routine. in case, fitting simple power law portion of galaxy power spectrum, , gsl giving me fits had chi^2/dof around 6.

in order solve this, figured out had been careless , x values data points didn't match x values fitting function being evaluated. easiest way fix create spline data values, , evaluate spline @ same x values fitting function evaluated. example:

#include <gsl/gsl_spline.h>     .     .     . std::vector< double > xvals; std::vector< double > yvals;  fin.open("somedatafile.dat", std::ios::in); while (!fin.eof()) {     double x, y;     fin >> x >> y;     xvals.push_back(x);     yvals.push_back(y); } fin.close();  gsl_spline *y = gsl_spline_alloc(gsl_interp_cspline, yvals.size()); gsl_interp_accel *acc = gsl_interp_accel_alloc();  gsl_spline_init(y, &xvals[0], &yvals[0], yvals.size());  double y[n];  (int = 0; < n; ++i) {     double x = xmin + i*dx; // xmin smallest x value , dx                             // (xmax-xmin)/n     y[i] = gsl_spline_eval(y, x, acc); } 

then in function difference between power law , data calculated, made sure use same xmin , dx x values same yi function in notation.

struct data {     size_t n;     double *y;     double xmin;     double xmax; };  int powerlaw(const gsl_vector *x, void *dat, gsl_vector *f) {     size_t n = ((data *) dat)->n;     double *y = ((data *) dat)->y;     double xmin = ((data *) dat)->xmin;     double xmax = ((data *) dat)->xmax;     double dx = (xmax-xmin)/double(n);      double = gsl_vector_get(x, 0);     double alpha = gsl_vector_get(x, 1);      (int = 0; < n; ++i) {         double xval = xmin + double(i)*dx;         double yi = a*pow(xval,alpha);         gsl_vector_set(f, i, yi - y[i]);     }      return gsl_success; } 

after that, values gnuplot , gsl agree quite nicely, gnuplot giving amplitude of 123.196 +/- 0.04484, , exponent of -1.13275 +/- 0.001903 , gsl giving 123.20464 +/- 0.98008 , -1.13272 +/- 0.00707. results of fits shown in graph below, fit gnuplot, , g(x) gsl (note: don't expect power laws precise match data, sufficient purposes). fits gnuplot , gsl virtually identical.

plot of data , fits gnuplot , gsl.

i have mentioned in comment question, have never had ask question here , have never answered one, didn't have enough rep.