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

   1 #include <stdio.h>
   2 #include <stdlib.h>
   3 #include <math.h>
   4 #include <gsl/gsl_bspline.h>
   5 #include <gsl/gsl_multifit.h>
   6 #include <gsl/gsl_rng.h>
   7 #include <gsl/gsl_randist.h>
   8 #include <gsl/gsl_statistics.h>
   9
  10 /* number of data points to fit */
  11 #define N        200
  12
  13 /* number of fit coefficients */
  14 #define NCOEFFS  12
  15
  16 /* nbreak = ncoeffs + 2 - k = ncoeffs - 2 since k = 4 */
  17 #define NBREAK   (NCOEFFS - 2)
  18
  19 int
  20 main (void)
  21 {
  22   const size_t n = N;
  23   const size_t ncoeffs = NCOEFFS;
  24   const size_t nbreak = NBREAK;
  25   size_t i, j;
  26   gsl_bspline_workspace *bw;
  27   gsl_vector *B;
  28   double dy;
  29   gsl_rng *r;
  30   gsl_vector *c, *w;
  31   gsl_vector *x, *y;
  32   gsl_matrix *X, *cov;
  33   gsl_multifit_linear_workspace *mw;
  34   double chisq;
  35   double Rsq;
  36   double dof;
  37
  38   gsl_rng_env_setup();
  39   r = gsl_rng_alloc(gsl_rng_default);
  40
  41   /* allocate a cubic bspline workspace (k = 4) */
  42   bw = gsl_bspline_alloc(4, nbreak);
  43   B = gsl_vector_alloc(ncoeffs);
  44
  45   x = gsl_vector_alloc(n);
  46   y = gsl_vector_alloc(n);
  47   X = gsl_matrix_alloc(n, ncoeffs);
  48   c = gsl_vector_alloc(ncoeffs);
  49   w = gsl_vector_alloc(n);
  50   cov = gsl_matrix_alloc(ncoeffs, ncoeffs);
  51   mw = gsl_multifit_linear_alloc(n, ncoeffs);
  52
  53   printf("#m=0,S=0\n");
  54   /* this is the data to be fitted */
  55   for (i = 0; i < n; ++i)
  56     {
  57       double sigma;
  58       double xi = (15.0 / (N - 1)) * i;
  59       double yi = cos(xi) * exp(-0.1 * xi);
  60
  61       sigma = 0.1 * yi;
  62       dy = gsl_ran_gaussian(r, sigma);
  63       yi += dy;
  64
  65       gsl_vector_set(x, i, xi);
  66       gsl_vector_set(y, i, yi);
  67       gsl_vector_set(w, i, 1.0 / (sigma * sigma));
  68
  69       printf("%f %f\n", xi, yi);
  70     }
  71
  72   /* use uniform breakpoints on [0, 15] */
  73   gsl_bspline_knots_uniform(0.0, 15.0, bw);
  74
  75   /* construct the fit matrix X */
  76   for (i = 0; i < n; ++i)
  77     {
  78       double xi = gsl_vector_get(x, i);
  79
  80       /* compute B_j(xi) for all j */
  81       gsl_bspline_eval(xi, B, bw);
  82
  83       /* fill in row i of X */
  84       for (j = 0; j < ncoeffs; ++j)
  85         {
  86           double Bj = gsl_vector_get(B, j);
  87           gsl_matrix_set(X, i, j, Bj);
  88         }
  89     }
  90
  91   /* do the fit */
  92   gsl_multifit_wlinear(X, w, y, c, cov, &chisq, mw);
  93
  94   dof = n - ncoeffs;
  95   Rsq = 1.0 - chisq / gsl_stats_wtss(w->data, 1, y->data, 1, y->size);
  96
  97   fprintf(stderr, "chisq/dof = %e, Rsq = %f\n", chisq / dof, Rsq);
  98
  99   /* output the smoothed curve */
 100   {
 101     double xi, yi, yerr;
 102
 103     printf("#m=1,S=0\n");
 104     for (xi = 0.0; xi < 15.0; xi += 0.1)
 105       {
 106         gsl_bspline_eval(xi, B, bw);
 107         gsl_multifit_linear_est(B, c, cov, &yi, &yerr);
 108         printf("%f %f\n", xi, yi);
 109       }
 110   }
 111
 112   gsl_rng_free(r);
 113   gsl_bspline_free(bw);
 114   gsl_vector_free(B);
 115   gsl_vector_free(x);
 116   gsl_vector_free(y);
 117   gsl_matrix_free(X);
 118   gsl_vector_free(c);
 119   gsl_vector_free(w);
 120   gsl_matrix_free(cov);
 121   gsl_multifit_linear_free(mw);
 122
 123   return 0;
 124 } /* main() */