multifit.texi

GNU Scientific Library,C语言开发的数值方面的函数库

the decrease in the function falls below machine precision

@item GSL_ETOLX
the change in the position vector falls below machine precision

@item GSL_ETOLG
the norm of the gradient, relative to the norm of the function, falls
below machine precision
@end table
@noindent
These error codes indicate that further iterations will be unlikely to
change the solution from its current value.

@end deffn

@deffn {Derivative Solver} gsl_multifit_fdfsolver_lmder
This is an unscaled version of the @sc{lmder} algorithm.  The elements of the
diagonal scaling matrix @math{D} are set to 1.  This algorithm may be
useful in circumstances where the scaled version of @sc{lmder} converges too
slowly, or the function is already scaled appropriately.
@end deffn

@node Minimization Algorithms without Derivatives
@section Minimization Algorithms without Derivatives

There are no algorithms implemented in this section at the moment.

@node Computing the covariance matrix of best fit parameters
@section Computing the covariance matrix of best fit parameters
@cindex covariance of best-fit parameters
@cindex best-fit parameters, covariance
@cindex least-squares, covariance of best-fit parameters

@deftypefun int gsl_multifit_covar (const gsl_matrix * @var{J}, double @var{epsrel}, gsl_matrix * @var{covar})
This function uses the Jacobian matrix @var{J} to compute the covariance
matrix of the best-fit parameters, @var{covar}.  The parameter
@var{epsrel} is used to remove linear-dependent columns when @var{J} is
rank deficient.

The covariance matrix is given by,

@tex
\beforedisplay
$$
C = (J^T J)^{-1}
$$
\afterdisplay
@end tex
@ifinfo
@example
covar = (J^T J)^@{-1@}
@end example
@end ifinfo
@noindent
and is computed by QR decomposition of J with column-pivoting.  Any
columns of @math{R} which satisfy 

@tex
\beforedisplay
$$
|R_{kk}| \leq epsrel |R_{11}|
$$
\afterdisplay
@end tex
@ifinfo
@example
|R_@{kk@}| <= epsrel |R_@{11@}|
@end example
@end ifinfo
@noindent
are considered linearly-dependent and are excluded from the covariance
matrix (the corresponding rows and columns of the covariance matrix are
set to zero).
@end deftypefun

@comment ============================================================

@node Example programs for Nonlinear Least-Squares Fitting
@section Examples

The following example program fits a weighted exponential model with
background to experimental data, @math{Y = A \exp(-\lambda t) + b}. The
first part of the program sets up the functions @code{expb_f} and
@code{expb_df} to calculate the model and its Jacobian.  The appropriate
fitting function is given by,

@tex
\beforedisplay
$$
f_i = ((A \exp(-\lambda t_i) + b) - y_i)/\sigma_i
$$
\afterdisplay
@end tex
@ifinfo
@example
f_i = ((A \exp(-\lambda t_i) + b) - y_i)/\sigma_i
@end example
@end ifinfo
@noindent
where we have chosen @math{t_i = i}.  The Jacobian matrix @math{J} is
the derivative of these functions with respect to the three parameters
(@math{A}, @math{\lambda}, @math{b}).  It is given by,

@tex
\beforedisplay
$$
J_{ij} = {\partial f_i \over \partial x_j}
$$
\afterdisplay
@end tex
@ifinfo
@example
J_@{ij@} = d f_i / d x_j
@end example
@end ifinfo
@noindent
where @math{x_0 = A}, @math{x_1 = \lambda} and @math{x_2 = b}.

@example
#include <stdlib.h>
#include <stdio.h>
#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>

struct data @{
  size_t n;
  double * y;
  double * sigma;
@};

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

  double A = gsl_vector_get (x, 0);
  double lambda = gsl_vector_get (x, 1);
  double b = gsl_vector_get (x, 2);

  size_t i;

  for (i = 0; i < n; i++)
    @{
      /* Model Yi = A * exp(-lambda * i) + b */
      double t = i;
      double Yi = A * exp (-lambda * t) + b;
      gsl_vector_set (f, i, (Yi - y[i])/sigma[i]);
    @}

  return GSL_SUCCESS;
@}

int
expb_df (const gsl_vector * x, void *params, 
         gsl_matrix * J)
@{
  size_t n = ((struct data *)params)->n;
  double *sigma = ((struct data *) params)->sigma;

  double A = gsl_vector_get (x, 0);
  double lambda = gsl_vector_get (x, 1);

  size_t i;

  for (i = 0; i < n; i++)
    @{
      /* Jacobian matrix J(i,j) = dfi / dxj, */
      /* where fi = (Yi - yi)/sigma[i],      */
      /*       Yi = A * exp(-lambda * i) + b  */
      /* and the xj are the parameters (A,lambda,b) */
      double t = i;
      double s = sigma[i];
      double e = exp(-lambda * t);
      gsl_matrix_set (J, i, 0, e/s); 
      gsl_matrix_set (J, i, 1, -t * A * e/s);
      gsl_matrix_set (J, i, 2, 1/s);


    @}
  return GSL_SUCCESS;
@}

int
expb_fdf (const gsl_vector * x, void *params,
          gsl_vector * f, gsl_matrix * J)
@{
  expb_f (x, params, f);
  expb_df (x, params, J);

  return GSL_SUCCESS;
@}
@end example
@noindent
The main part of the program sets up a Levenberg-Marquardt solver and
some simulated random data. The data uses the known parameters
(1.0,5.0,0.1) combined with gaussian noise (standard deviation = 0.1)
over a range of 40 timesteps. The initial guess for the parameters is
chosen as (0.0, 1.0, 0.0).

@example
#define N 40

int
main (void)
@{
  const gsl_multifit_fdfsolver_type *T;
  gsl_multifit_fdfsolver *s;

  int status;
  size_t i, iter = 0;

  const size_t n = N;
  const size_t p = 3;

  gsl_matrix *covar = gsl_matrix_alloc (p, p);

  double y[N], sigma[N];

  struct data d = @{ n, y, sigma@};
  
  gsl_multifit_function_fdf f;

  double x_init[3] = @{ 1.0, 0.0, 0.0 @};

  gsl_vector_view x = gsl_vector_view_array (x_init, p);

  const gsl_rng_type * type;
  gsl_rng * r;

  gsl_rng_env_setup();

  type = gsl_rng_default;
  r = gsl_rng_alloc (type);

  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 */

  for (i = 0; i < n; i++)
    @{
      double t = i;
      y[i] = 1.0 + 5 * exp (-0.1 * t) 
                 + gsl_ran_gaussian (r, 0.1);
      sigma[i] = 0.1;
      printf ("data: %d %g %g\n", i, y[i], sigma[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))

  printf ("A      = %.5f +/- %.5f\n", FIT(0), ERR(0));
  printf ("lambda = %.5f +/- %.5f\n", FIT(1), ERR(1));
  printf ("b      = %.5f +/- %.5f\n", FIT(2), ERR(2));

  @{ 
    double chi = gsl_blas_dnrm2(s->f);
    printf("chisq/dof = %g\n",  pow(chi, 2.0)/ (n - p));
  @}

  printf ("status = %s\n", gsl_strerror (status));

  gsl_multifit_fdfsolver_free (s);
  return 0;
@}

void
print_state (size_t iter, gsl_multifit_fdfsolver * s)
@{
  printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f "
          "|f(x)| = %g\n",
          iter,
          gsl_vector_get (s->x, 0), 
          gsl_vector_get (s->x, 1),
          gsl_vector_get (s->x, 2), 
          gsl_blas_dnrm2 (s->f));
@}
@end example
@noindent
The iteration terminates when the change in x is smaller than 0.0001, as
both an absolute and relative change.  Here are the results of running
the program:

@smallexample
iter: 0 x = 1.00000000 0.00000000 0.00000000 |f(x)| = 118.574
iter: 1 x = 1.64919392 0.01780040 0.64919392 |f(x)| = 77.2068
iter: 2 x = 2.86269020 0.08032198 1.45913464 |f(x)| = 38.0579
iter: 3 x = 4.97908864 0.11510525 1.06649948 |f(x)| = 10.1548
iter: 4 x = 5.03295496 0.09912462 1.00939075 |f(x)| = 6.4982
iter: 5 x = 5.05811477 0.10055914 0.99819876 |f(x)| = 6.33121
iter: 6 x = 5.05827645 0.10051697 0.99756444 |f(x)| = 6.33119
iter: 7 x = 5.05828006 0.10051819 0.99757710 |f(x)| = 6.33119

A      = 5.05828 +/- 0.05983
lambda = 0.10052 +/- 0.00309
b      = 0.99758 +/- 0.03944
chisq/dof = 1.08335
status = success
@end smallexample
@noindent
The approximate values of the parameters are found correctly, and the
chi-squared value indicates a good fit (the chi-squared per degree of
freedom is approximately 1).  In this case the errors on the parameters
can be estimated from the square roots of the diagonal elements of the
covariance matrix.  If the chi-squared value indicates a poor fit then
error estimates obtained from the covariance matrix are not valid, since
the Gaussian approximation would not apply.

@iftex
@sp 1
@center @image{fit-exp,3.4in}
@end iftex

@node References and Further Reading for Nonlinear Least-Squares Fitting
@section References and Further Reading

The @sc{minpack} algorithm is described in the following article,

@itemize @asis
@item
J.J. Mor@'e, @cite{The Levenberg-Marquardt Algorithm: Implementation and
Theory}, Lecture Notes in Mathematics, v630 (1978), ed G. Watson.
@end itemize

@noindent
The following paper is also relevant to the algorithms described in this
section,

@itemize @asis
@item 
J.J. Mor@'e, B.S. Garbow, K.E. Hillstrom, "Testing Unconstrained
Optimization Software", ACM Transactions on Mathematical Software, Vol
7, No 1 (1981), p 17-41.
@end itemize