htswanalysis/MACS/lib/gsl/gsl-1.11/doc/examples/nlfit.c

   1 #include <stdlib.h>
   2 #include <stdio.h>
   3 #include <gsl/gsl_rng.h>
   4 #include <gsl/gsl_randist.h>
   5 #include <gsl/gsl_vector.h>
   6 #include <gsl/gsl_blas.h>
   7 #include <gsl/gsl_multifit_nlin.h>
   8
   9 #include "expfit.c"
  10
  11 #define N 40
  12
  13 void print_state (size_t iter, gsl_multifit_fdfsolver * s);
  14
  15 int
  16 main (void)
  17 {
  18   const gsl_multifit_fdfsolver_type *T;
  19   gsl_multifit_fdfsolver *s;
  20   int status;
  21   unsigned int i, iter = 0;
  22   const size_t n = N;
  23   const size_t p = 3;
  24
  25   gsl_matrix *covar = gsl_matrix_alloc (p, p);
  26   double y[N], sigma[N];
  27   struct data d = { n, y, sigma};
  28   gsl_multifit_function_fdf f;
  29   double x_init[3] = { 1.0, 0.0, 0.0 };
  30   gsl_vector_view x = gsl_vector_view_array (x_init, p);
  31   const gsl_rng_type * type;
  32   gsl_rng * r;
  33
  34   gsl_rng_env_setup();
  35
  36   type = gsl_rng_default;
  37   r = gsl_rng_alloc (type);
  38
  39   f.f = &expb_f;
  40   f.df = &expb_df;
  41   f.fdf = &expb_fdf;
  42   f.n = n;
  43   f.p = p;
  44   f.params = &d;
  45
  46   /* This is the data to be fitted */
  47
  48   for (i = 0; i < n; i++)
  49     {
  50       double t = i;
  51       y[i] = 1.0 + 5 * exp (-0.1 * t)
  52                  + gsl_ran_gaussian (r, 0.1);
  53       sigma[i] = 0.1;
  54       printf ("data: %u %g %g\n", i, y[i], sigma[i]);
  55     };
  56
  57   T = gsl_multifit_fdfsolver_lmsder;
  58   s = gsl_multifit_fdfsolver_alloc (T, n, p);
  59   gsl_multifit_fdfsolver_set (s, &f, &x.vector);
  60
  61   print_state (iter, s);
  62
  63   do
  64     {
  65       iter++;
  66       status = gsl_multifit_fdfsolver_iterate (s);
  67
  68       printf ("status = %s\n", gsl_strerror (status));
  69
  70       print_state (iter, s);
  71
  72       if (status)
  73         break;
  74
  75       status = gsl_multifit_test_delta (s->dx, s->x,
  76                                         1e-4, 1e-4);
  77     }
  78   while (status == GSL_CONTINUE && iter < 500);
  79
  80   gsl_multifit_covar (s->J, 0.0, covar);
  81
  82 #define FIT(i) gsl_vector_get(s->x, i)
  83 #define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
  84
  85   {
  86     double chi = gsl_blas_dnrm2(s->f);
  87     double dof = n - p;
  88     double c = GSL_MAX_DBL(1, chi / sqrt(dof));
  89
  90     printf("chisq/dof = %g\n",  pow(chi, 2.0) / dof);
  91
  92     printf ("A      = %.5f +/- %.5f\n", FIT(0), c*ERR(0));
  93     printf ("lambda = %.5f +/- %.5f\n", FIT(1), c*ERR(1));
  94     printf ("b      = %.5f +/- %.5f\n", FIT(2), c*ERR(2));
  95   }
  96
  97   printf ("status = %s\n", gsl_strerror (status));
  98
  99   gsl_multifit_fdfsolver_free (s);
 100   gsl_matrix_free (covar);
 101   gsl_rng_free (r);
 102   return 0;
 103 }
 104
 105 void
 106 print_state (size_t iter, gsl_multifit_fdfsolver * s)
 107 {
 108   printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f "
 109           "|f(x)| = %g\n",
 110           iter,
 111           gsl_vector_get (s->x, 0),
 112           gsl_vector_get (s->x, 1),
 113           gsl_vector_get (s->x, 2),
 114           gsl_blas_dnrm2 (s->f));
 115 }