fitting.texi

开放gsl矩阵运算

int
main (void)
@{
  int i, n = 4;
  double x[4] = @{ 1970, 1980, 1990, 2000 @};
  double y[4] = @{   12,   11,   14,   13 @};
  double w[4] = @{  0.1,  0.2,  0.3,  0.4 @};

  double c0, c1, cov00, cov01, cov11, chisq;

  gsl_fit_wlinear (x, 1, w, 1, y, 1, n, 
                   &c0, &c1, &cov00, &cov01, &cov11, 
                   &chisq);

  printf("# best fit: Y = %g + %g X\n", c0, c1);
  printf("# covariance matrix:\n");
  printf("# [ %g, %g\n#   %g, %g]\n", 
         cov00, cov01, cov01, cov11);
  printf("# chisq = %g\n", chisq);

  for (i = 0; i < n; i++)
    printf("data: %g %g %g\n", 
                  x[i], y[i], 1/sqrt(w[i]));

  printf("\n");

  for (i = -30; i < 130; i++)
    @{
      double xf = x[0] + (i/100.0) * (x[n-1] - x[0]);
      double yf, yf_err;

      gsl_fit_linear_est (xf, 
                          c0, c1, 
                          cov00, cov01, cov11, 
                          &yf, &yf_err);

      printf("fit: %g %g\n", xf, yf);
      printf("hi : %g %g\n", xf, yf + yf_err);
      printf("lo : %g %g\n", xf, yf - yf_err);
    @}
  return 0;
@}
@end example
@noindent
The following commands extract the data from the output of the program
and display it using the @sc{gnu} plotutils @code{graph} utility, 

@example
$ ./demo > tmp
$ more tmp
# best fit: Y = -106.6 + 0.06 X
# covariance matrix:
# [ 39602, -19.9
#   -19.9, 0.01]
# chisq = 0.8

$ for n in data fit hi lo ; 
   do 
     grep "^$n" tmp | cut -d: -f2 > $n ; 
   done
$ graph -T X -X x -Y y -y 0 20 -m 0 -S 2 -Ie data 
     -S 0 -I a -m 1 fit -m 2 hi -m 2 lo
@end example

@iftex
@sp 1
@center @image{fit-wlinear,4in}
@end iftex

The next program performs a quadratic fit @math{y = c_0 + c_1 x + c_2
x^2} to a weighted dataset using the generalised linear fitting function
@code{gsl_multifit_wlinear}.  The model matrix @math{X} for a quadratic
fit is given by,

@tex
\beforedisplay
$$
X=\pmatrix{1&x_0&x_0^2\cr
1&x_1&x_1^2\cr
1&x_2&x_2^2\cr
\dots&\dots&\dots\cr}
$$
\afterdisplay
@end tex
@ifinfo
@example
X = [ 1   , x_0  , x_0^2 ;
      1   , x_1  , x_1^2 ;
      1   , x_2  , x_2^2 ;
      ... , ...  , ...   ]
@end example
@end ifinfo
@noindent
where the column of ones corresponds to the constant term @math{c_0}.
The two remaining columns corresponds to the terms @math{c_1 x} and and
@math{c_2 x^2}.

The program reads @var{n} lines of data in the format (@var{x}, @var{y},
@var{err}) where @var{err} is the error (standard deviation) in the
value @var{y}.

@example
#include <stdio.h>
#include <gsl/gsl_multifit.h>

int
main (int argc, char **argv)
@{
  int i, n;
  double xi, yi, ei, chisq;
  gsl_matrix *X, *cov;
  gsl_vector *y, *w, *c;

  if (argc != 2)
    @{
      fprintf(stderr,"usage: fit n < data\n");
      exit (-1);
    @}

  n = atoi(argv[1]);

  X = gsl_matrix_alloc (n, 3);
  y = gsl_vector_alloc (n);
  w = gsl_vector_alloc (n);

  c = gsl_vector_alloc (3);
  cov = gsl_matrix_alloc (3, 3);

  for (i = 0; i < n; i++)
    @{
      int count = fscanf(stdin, "%lg %lg %lg",
                         &xi, &yi, &ei);

      if (count != 3)
        @{
          fprintf(stderr, "error reading file\n");
          exit(-1);
        @}

      printf("%g %g +/- %g\n", xi, yi, ei);
      
      gsl_matrix_set (X, i, 0, 1.0);
      gsl_matrix_set (X, i, 1, xi);
      gsl_matrix_set (X, i, 2, xi*xi);
      
      gsl_vector_set (y, i, yi);
      gsl_vector_set (w, i, 1.0/(ei*ei));
    @}

  @{
    gsl_multifit_linear_workspace * work 
      = gsl_multifit_linear_alloc (n, 3);
    gsl_multifit_wlinear (X, w, y, c, cov,
                          &chisq, work);
    gsl_multifit_linear_free (work);
  @}

#define C(i) (gsl_vector_get(c,(i)))
#define COV(i,j) (gsl_matrix_get(cov,(i),(j)))

  @{
    printf("# best fit: Y = %g + %g X + %g X^2\n", 
           C(0), C(1), C(2));

    printf("# covariance matrix:\n");
    printf("[ %+.5e, %+.5e, %+.5e  \n",
              COV(0,0), COV(0,1), COV(0,2));
    printf("  %+.5e, %+.5e, %+.5e  \n", 
              COV(1,0), COV(1,1), COV(1,2));
    printf("  %+.5e, %+.5e, %+.5e ]\n", 
              COV(2,0), COV(2,1), COV(2,2));
    printf("# chisq = %g\n", chisq);
  @}
  return 0;
@}
@end example
@noindent
A suitable set of data for fitting can be generated using the following
program.  It outputs a set of points with gaussian errors from the curve
@math{y = e^x} in the region @math{0 < x < 2}.

@example
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_randist.h>

int
main (void)
@{
  double x;
  const gsl_rng_type * T;
  gsl_rng * r;
  
  gsl_rng_env_setup();
  
  T = gsl_rng_default;
  r = gsl_rng_alloc(T);

  for (x = 0.1; x < 2; x+= 0.1)
    @{
      double y0 = exp(x);
      double sigma = 0.1*y0;
      double dy = gsl_ran_gaussian(r, sigma)

      printf("%g %g %g\n", x, y0 + dy, sigma);
    @}
  return 0;
@}
@end example
@noindent
The data can be prepared by running the resulting executable program,

@example
$ ./generate > exp.dat
$ more exp.dat
0.1 0.97935 0.110517
0.2 1.3359 0.12214
0.3 1.52573 0.134986
0.4 1.60318 0.149182
0.5 1.81731 0.164872
0.6 1.92475 0.182212
....
@end example
@noindent
To fit the data use the previous program, with the number of data points
given as the first argument.  In this case there are 19 data points.

@example
$ ./fit 19 < exp.dat
0.1 0.97935 +/- 0.110517
0.2 1.3359 +/- 0.12214
...
# best fit: Y = 1.02318 + 0.956201 X + 0.876796 X^2
# covariance matrix:
[ +1.25612e-02, -3.64387e-02, +1.94389e-02  
  -3.64387e-02, +1.42339e-01, -8.48761e-02  
  +1.94389e-02, -8.48761e-02, +5.60243e-02 ]
# chisq = 23.0987
@end example
@noindent
The parameters of the quadratic fit match the coefficients of the
expansion of @math{e^x}, taking into account the errors on the
parameters and the @math{O(x^3)} difference between the exponential and
quadratic functions for the larger values of @math{x}.  The errors on
the parameters are given by the square-root of the corresponding
diagonal elements of the covariance matrix.  The chi-squared per degree
of freedom is 1.4, indicating a reasonable fit to the data.

@iftex
@sp 1
@center @image{fit-wlinear2,4in}
@end iftex

@node Fitting References and Further Reading
@section References and Further Reading
@noindent
A summary of formulas and techniques for least squares fitting can be
found in the "Statistics" chapter of the Annual Review of Particle
Physics prepared by the Particle Data Group.

@itemize @asis
@item
@cite{Review of Particle Properties}
R.M. Barnett et al., Physical Review D54, 1 (1996)
@url{http://pdg.lbl.gov/}
@end itemize
@noindent
The Review of Particle Physics is available online at the website given
above.

@noindent
The tests used to prepare these routines are based on the NIST
Statistical Reference Datasets. The datasets and their documentation are
available from NIST at the following website,

@center @url{http://www.nist.gov/itl/div898/strd/index.html}.