Re: [Help-gsl] nonlinear fit

[I.I. Stefan]

Hey,

Thanks. I implemented that, but now the routine is not converging. This is 
the error I'm getting:

status = the iteration has not converged yet

And the routine is below.
Thank you for your time!

Best wishes,
Irina


//fits a Gaussian function f(SB) to the Fake_Gaus_data text files through 
the gsl non-linear least squares. Doesn't properly work though...

    #include<iostream>
    #include<sstream>
    #include<fstream>
    #include<cmath>
    #include<vector>
    #include <stdlib.h>
    #include <stdio.h>
    #include <string>
    #include <gsl/gsl_rng.h>
    #include <gsl/gsl_randist.h>
    #include <gsl/gsl_vector.h>
    #include <gsl/gsl_blas.h>
    #include <gsl/gsl_multifit_nlin.h>
    #include <gsl/gsl_math.h>
    
	using namespace std;
    

    struct data {
      size_t n;
      double * y;
      double * SB;
    };
    
    int expb_f (const gsl_vector * x, void *data, gsl_vector * f)
    { size_t n = ((struct data *)data)-> n;
      double *SB = ((struct data *)data)->SB;
      double *y = ((struct data *)data)->y;

      double sigma = gsl_vector_get(x,0);
      double mu = gsl_vector_get(x, 1);       

      size_t i;
    
      for (i = 0; i < n; i++)
        {
          /* Model Yi = 
(1/sqrt(2*M_PI*sigma*sigma))*exp(-(t-mu)*(t-mu)/(2*sigma*sigma)) */
          double t = SB[i];
          double Yi = 
(1/sqrt(2*M_PI*sigma*sigma))*exp(-(t-mu)*(t-mu)/(2*sigma*sigma)) ;
          gsl_vector_set (f, i, Yi-y[i]);
        }
    
      return GSL_SUCCESS;
    }
    
    int
    expb_df (const gsl_vector * x, void *data, gsl_matrix * J)
    {
      size_t n = ((struct data *)data)->n;
      double *SB = ((struct data *)data)->SB;
    
      double sigma = gsl_vector_get (x, 0);
      double mu = gsl_vector_get (x, 1);
    
      size_t i;
    
      for (i = 0; i < n; i++)
        {
          /* Jacobian matrix J(i,j) = dfi / dxj, */
          /* Yi = 
(1/sqrt(2*M_PI*sigma*sigma))*exp(-(t-mu)*(t-mu)/(2*sigma*sigma)) */
          /* and the xj are the parameters (sigma, mu) */
          double t = SB[i];
          double e = exp(-(t-mu)*(t-mu)/(2*sigma*sigma));
           double A = (1/sqrt(2*M_PI*sigma*sigma));
          gsl_matrix_set (J, i, 0, 
((A/sigma)*e*((t-mu)*(t-mu)/(sigma*sigma)-1)));
          gsl_matrix_set (J, i, 1, (-(A/(sigma*sigma) * e * (t-mu))));
        }
      return GSL_SUCCESS;
    }
    
    int
    expb_fdf (const gsl_vector * x, void *data, gsl_vector * f, gsl_matrix 
* J)
    {
      expb_f (x, data, f);
      expb_df (x, data, J);
    
      return GSL_SUCCESS;
    }

//The main part of the program sets up a Levenberg-Marquardt solver. The 
initial guess for the parameters is chosen as (0.8, -18).

    void print_state (size_t iter, gsl_multifit_fdfsolver * s);
    
    int
    main (void)
    {
     unsigned int i=0, iter = 0;
     string line, temp;

ifstream read;
read.open ("fake_Gaus_data.txt");
if(read.is_open())
        while(!read.eof()){
                getline(read,line);
                i++;}
        read.close();

    #define N 200

      const gsl_multifit_fdfsolver_type *T;
      gsl_multifit_fdfsolver *s;
      int status;
      const size_t n = N;
      const size_t p = 2;
      int j, k;
    
      gsl_matrix *covar = gsl_matrix_alloc (p, p);
      double y[N], SB[N];
      struct data d = { n, y, SB};
      gsl_multifit_function_fdf f;
      double x_init[2] = {2, 1};
      gsl_vector_view x = gsl_vector_view_array (x_init, p);

      f.f = &expb_f;
      f.df = &expb_df;
      f.fdf = &expb_fdf;
      f.n = n;
      f.p = p;
      f.params = &d;
    
      /* This is the data to be fitted */
        read.open ("fake_Gaus_data.txt");
        if(read.is_open()){
                        i=0;
                        getline(read,line);
                        while(!read.eof()){
                                j=0;
                                k=0;
                                getline(read, line);
                                if (line.length()!=0){
                                        j=line.find(" ", 0);                    
      
                                        temp=line.substr(0, j);
                                        SB[i]= strtod(temp.c_str(), NULL);
                                        k=line.length();
                                        temp=line.substr(j+1, (k-j));
                                        y[i]=strtod(temp.c_str(), NULL);
                                        i++;}}
        read.close();}
        else cout<<"File was not open"<<endl;
        cout<<endl<<(i-1)<<endl<<endl;

      for (i = 0; i < n; i++)
        {
          printf ("data: %g %g\n", SB[i], y[i]);
        }
    
      T = gsl_multifit_fdfsolver_lmsder;
      s = gsl_multifit_fdfsolver_alloc (T, n, p);
      gsl_multifit_fdfsolver_set (s, &f, &x.vector);
    
      print_state (iter, s);
    
      do
        {
          iter++;
          status = gsl_multifit_fdfsolver_iterate (s);
    
          printf ("status = %s\n", gsl_strerror (status));
    
          print_state (iter, s);
    
          if (status)
            break;
    
          status = gsl_multifit_test_delta (s->dx, s->x, 1e-4, 1e-4);
        }
      while (status == GSL_CONTINUE && iter < 500);
    
      gsl_multifit_covar (s->J, 0.0, covar);
    
    #define FIT(i) gsl_vector_get(s->x, i)
    #define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
    
      { 
        double chi = gsl_blas_dnrm2(s->f);
        double dof = n - p;
        double c = GSL_MAX_DBL(1, chi / sqrt(dof)); 
    
        printf("chisq/dof = %g\n",  pow(chi, 2.0) / dof);
    
        printf ("sigma  = %.5f +/- %.5f\n", FIT(0), c*ERR(0));
        printf ("mu     = %.5f +/- %.5f\n", FIT(1), c*ERR(1));
      }
    
      printf ("status = %s\n", gsl_strerror (status));
    
      gsl_multifit_fdfsolver_free (s);
      gsl_matrix_free (covar);
      return 0;
    }
    
    void
    print_state (size_t iter, gsl_multifit_fdfsolver * s)
    {
      printf ("iter: %3u x = % 15.8f % 15.8f "
              "|f(x)| = %g\n",
              iter,
              gsl_vector_get (s->x, 0), 
              gsl_vector_get (s->x, 1),
              gsl_blas_dnrm2 (s->f));
    }




On Mar 9 2011, Francesco Abbate wrote:

> 2011/3/8 Irina Stefan <address@hidden>:
> > Hi,
> >
> > This is probably a very easy question, but I'm not very comfortable 
> > with gsl, so here it goes. I'm writing a routine to fit a Gaussian to a 
> > set of data points. I'm pretty much modifying the example from the gsl 
> > manual. However, I want to perform an unweighted fit and I'm not sure 
> > what function to use as f_i.
>
> Hi,
>
> this is actually an easy question :-)
>
> If you want to perform an unweighted fit just set f_i to Y_i - y_i,
> i.e., to the difference between the obsered value and the predicted
> one. Just don't forget to do the same thing for f_i and the Jacobian.
>
> In the example they divides by sigma_i so you have just to omit the
> division by sigma_i both for f_i and the jacobian.