3
<title>Example programs for Nonlinear Least-Squares Fitting - GNU Scientific Library -- Reference Manual</title>
4
<meta http-equiv="Content-Type" content="text/html">
5
<meta name="description" content="GNU Scientific Library -- Reference Manual">
6
<meta name="generator" content="makeinfo 4.8">
7
<link title="Top" rel="start" href="index.html#Top">
8
<link rel="up" href="Nonlinear-Least_002dSquares-Fitting.html#Nonlinear-Least_002dSquares-Fitting" title="Nonlinear Least-Squares Fitting">
9
<link rel="prev" href="Computing-the-covariance-matrix-of-best-fit-parameters.html#Computing-the-covariance-matrix-of-best-fit-parameters" title="Computing the covariance matrix of best fit parameters">
10
<link rel="next" href="References-and-Further-Reading-for-Nonlinear-Least_002dSquares-Fitting.html#References-and-Further-Reading-for-Nonlinear-Least_002dSquares-Fitting" title="References and Further Reading for Nonlinear Least-Squares Fitting">
11
<link href="http://www.gnu.org/software/texinfo/" rel="generator-home" title="Texinfo Homepage">
13
Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 The GSL Team.
15
Permission is granted to copy, distribute and/or modify this document
16
under the terms of the GNU Free Documentation License, Version 1.2 or
17
any later version published by the Free Software Foundation; with the
18
Invariant Sections being ``GNU General Public License'' and ``Free Software
19
Needs Free Documentation'', the Front-Cover text being ``A GNU Manual'',
20
and with the Back-Cover Text being (a) (see below). A copy of the
21
license is included in the section entitled ``GNU Free Documentation
24
(a) The Back-Cover Text is: ``You have freedom to copy and modify this
25
GNU Manual, like GNU software.''-->
26
<meta http-equiv="Content-Style-Type" content="text/css">
27
<style type="text/css"><!--
28
pre.display { font-family:inherit }
29
pre.format { font-family:inherit }
30
pre.smalldisplay { font-family:inherit; font-size:smaller }
31
pre.smallformat { font-family:inherit; font-size:smaller }
32
pre.smallexample { font-size:smaller }
33
pre.smalllisp { font-size:smaller }
34
span.sc { font-variant:small-caps }
35
span.roman { font-family:serif; font-weight:normal; }
36
span.sansserif { font-family:sans-serif; font-weight:normal; }
42
<a name="Example-programs-for-Nonlinear-Least-Squares-Fitting"></a>
43
<a name="Example-programs-for-Nonlinear-Least_002dSquares-Fitting"></a>
44
Next: <a rel="next" accesskey="n" href="References-and-Further-Reading-for-Nonlinear-Least_002dSquares-Fitting.html#References-and-Further-Reading-for-Nonlinear-Least_002dSquares-Fitting">References and Further Reading for Nonlinear Least-Squares Fitting</a>,
45
Previous: <a rel="previous" accesskey="p" href="Computing-the-covariance-matrix-of-best-fit-parameters.html#Computing-the-covariance-matrix-of-best-fit-parameters">Computing the covariance matrix of best fit parameters</a>,
46
Up: <a rel="up" accesskey="u" href="Nonlinear-Least_002dSquares-Fitting.html#Nonlinear-Least_002dSquares-Fitting">Nonlinear Least-Squares Fitting</a>
50
<h3 class="section">37.9 Examples</h3>
52
<p>The following example program fits a weighted exponential model with
53
background to experimental data, Y = A \exp(-\lambda t) + b. The
54
first part of the program sets up the functions <code>expb_f</code> and
55
<code>expb_df</code> to calculate the model and its Jacobian. The appropriate
56
fitting function is given by,
57
where we have chosen t_i = i. The Jacobian matrix J is
58
the derivative of these functions with respect to the three parameters
59
(A, \lambda, b). It is given by,
60
where x_0 = A, x_1 = \lambda and x_2 = b.
62
<pre class="example"><pre class="verbatim"> /* expfit.c -- model functions for exponential + background */
71
expb_f (const gsl_vector * x, void *data,
74
size_t n = ((struct data *)data)->n;
75
double *y = ((struct data *)data)->y;
76
double *sigma = ((struct data *) data)->sigma;
78
double A = gsl_vector_get (x, 0);
79
double lambda = gsl_vector_get (x, 1);
80
double b = gsl_vector_get (x, 2);
84
for (i = 0; i < n; i++)
86
/* Model Yi = A * exp(-lambda * i) + b */
88
double Yi = A * exp (-lambda * t) + b;
89
gsl_vector_set (f, i, (Yi - y[i])/sigma[i]);
96
expb_df (const gsl_vector * x, void *data,
99
size_t n = ((struct data *)data)->n;
100
double *sigma = ((struct data *) data)->sigma;
102
double A = gsl_vector_get (x, 0);
103
double lambda = gsl_vector_get (x, 1);
107
for (i = 0; i < n; i++)
109
/* Jacobian matrix J(i,j) = dfi / dxj, */
110
/* where fi = (Yi - yi)/sigma[i], */
111
/* Yi = A * exp(-lambda * i) + b */
112
/* and the xj are the parameters (A,lambda,b) */
115
double e = exp(-lambda * t);
116
gsl_matrix_set (J, i, 0, e/s);
117
gsl_matrix_set (J, i, 1, -t * A * e/s);
118
gsl_matrix_set (J, i, 2, 1/s);
124
expb_fdf (const gsl_vector * x, void *data,
125
gsl_vector * f, gsl_matrix * J)
128
expb_df (x, data, J);
133
<p class="noindent">The main part of the program sets up a Levenberg-Marquardt solver and
134
some simulated random data. The data uses the known parameters
135
(1.0,5.0,0.1) combined with gaussian noise (standard deviation = 0.1)
136
over a range of 40 timesteps. The initial guess for the parameters is
137
chosen as (0.0, 1.0, 0.0).
139
<pre class="example"><pre class="verbatim"> #include <stdlib.h>
140
#include <stdio.h>
141
#include <gsl/gsl_rng.h>
142
#include <gsl/gsl_randist.h>
143
#include <gsl/gsl_vector.h>
144
#include <gsl/gsl_blas.h>
145
#include <gsl/gsl_multifit_nlin.h>
151
void print_state (size_t iter, gsl_multifit_fdfsolver * s);
156
const gsl_multifit_fdfsolver_type *T;
157
gsl_multifit_fdfsolver *s;
165
gsl_matrix *covar = gsl_matrix_alloc (p, p);
166
double y[N], sigma[N];
167
struct data d = { n, y, sigma};
168
gsl_multifit_function_fdf f;
169
double x_init[3] = { 1.0, 0.0, 0.0 };
170
gsl_vector_view x = gsl_vector_view_array (x_init, p);
171
const gsl_rng_type * type;
176
type = gsl_rng_default;
177
r = gsl_rng_alloc (type);
181
f.fdf = &expb_fdf;
186
/* This is the data to be fitted */
188
for (i = 0; i < n; i++)
191
y[i] = 1.0 + 5 * exp (-0.1 * t)
192
+ gsl_ran_gaussian (r, 0.1);
194
printf ("data: %d %g %g\n", i, y[i], sigma[i]);
197
T = gsl_multifit_fdfsolver_lmsder;
198
s = gsl_multifit_fdfsolver_alloc (T, n, p);
199
gsl_multifit_fdfsolver_set (s, &f, &x.vector);
201
print_state (iter, s);
206
status = gsl_multifit_fdfsolver_iterate (s);
208
printf ("status = %s\n", gsl_strerror (status));
210
print_state (iter, s);
215
status = gsl_multifit_test_delta (s->dx, s->x,
218
while (status == GSL_CONTINUE && iter < 500);
220
gsl_multifit_covar (s->J, 0.0, covar);
222
#define FIT(i) gsl_vector_get(s->x, i)
223
#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
226
double chi = gsl_blas_dnrm2(s->f);
228
double c = GSL_MAX_DBL(1, chi / sqrt(dof));
230
printf("chisq/dof = %g\n", pow(chi, 2.0) / dof);
232
printf ("A = %.5f +/- %.5f\n", FIT(0), c*ERR(0));
233
printf ("lambda = %.5f +/- %.5f\n", FIT(1), c*ERR(1));
234
printf ("b = %.5f +/- %.5f\n", FIT(2), c*ERR(2));
237
printf ("status = %s\n", gsl_strerror (status));
239
gsl_multifit_fdfsolver_free (s);
244
print_state (size_t iter, gsl_multifit_fdfsolver * s)
246
printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f "
249
gsl_vector_get (s->x, 0),
250
gsl_vector_get (s->x, 1),
251
gsl_vector_get (s->x, 2),
252
gsl_blas_dnrm2 (s->f));
255
<p class="noindent">The iteration terminates when the change in x is smaller than 0.0001, as
256
both an absolute and relative change. Here are the results of running
259
<pre class="smallexample"> iter: 0 x=1.00000000 0.00000000 0.00000000 |f(x)|=117.349
261
iter: 1 x=1.64659312 0.01814772 0.64659312 |f(x)|=76.4578
263
iter: 2 x=2.85876037 0.08092095 1.44796363 |f(x)|=37.6838
265
iter: 3 x=4.94899512 0.11942928 1.09457665 |f(x)|=9.58079
267
iter: 4 x=5.02175572 0.10287787 1.03388354 |f(x)|=5.63049
269
iter: 5 x=5.04520433 0.10405523 1.01941607 |f(x)|=5.44398
271
iter: 6 x=5.04535782 0.10404906 1.01924871 |f(x)|=5.44397
273
A = 5.04536 +/- 0.06028
274
lambda = 0.10405 +/- 0.00316
275
b = 1.01925 +/- 0.03782
278
<p class="noindent">The approximate values of the parameters are found correctly, and the
279
chi-squared value indicates a good fit (the chi-squared per degree of
280
freedom is approximately 1). In this case the errors on the parameters
281
can be estimated from the square roots of the diagonal elements of the
284
<p>If the chi-squared value shows a poor fit (i.e. <!-- {$\chi^2/(n-p) \gg 1$} -->
285
chi^2/dof >> 1) then the error estimates obtained from the
286
covariance matrix will be too small. In the example program the error estimates
287
are multiplied by <!-- {$\sqrt{\chi^2/(n-p)}$} -->
288
\sqrt{\chi^2/dof} in this case, a common way of increasing the
289
errors for a poor fit. Note that a poor fit will result from the use
290
an inappropriate model, and the scaled error estimates may then
291
be outside the range of validity for gaussian errors.
added
removed
Example-programs-for-Nonlinear-Least_002dSquares-Fitting.html
