User Tools

Site Tools


Curve fitting with lmcurve()

This example shows how to fit a data set y(t) with a parametric function f(t;p).

/* demo/curve1.c */
 
#include "lmcurve.h"
#include <stdio.h>
 
/* model function: a parabola */
 
double f( double t, const double *p )
{
    return p[0] + p[1]*t + p[2]*t*t;
}
 
int main()
{
    int n = 3; /* number of parameters in model function f */
    double par[3] = { 100, 0, -10 }; /* really bad starting value */
 
    /* data points: a slightly distorted standard parabola */
    int m = 9;
    int i;
    double t[9] = { -4., -3., -2., -1.,  0., 1.,  2.,  3.,  4. };
    double y[9] = { 16.6, 9.9, 4.4, 1.1, 0., 1.1, 4.2, 9.3, 16.4 };
 
    lm_control_struct control = lm_control_double;
    lm_princon_struct princon = lm_princon_std;
    lm_status_struct status;
 
    printf( "Fitting ...\n" );
    /* now the call to lmfit */
    lmcurve( n, par, m, t, y, f, &control, &princon, &status );
 
    printf( "Results:\n" );
    printf( "status after %d function evaluations:\n  %s\n",
            status.nfev, lm_infmsg[status.info] );
 
    printf("obtained parameters:\n");
    for ( i = 0; i < n; ++i)
        printf("  par[%i] = %12g\n", i, par[i]);
    printf("obtained norm:\n  %12g\n", status.fnorm );
 
    printf("fitting data as follows:\n");
    for ( i = 0; i < m; ++i)
        printf( "  t[%2d]=%4g y=%6g fit=%10g residue=%12g\n",
                i, t[i], y[i], f(t[i],par), y[i] - f(t[i],par) );
 
    return 0;
}

On a standard 32 or 64 bit PC, it produces about the following output:

    Fitting ...
    Results:
    status after 21 function evaluations:
      success (the relative error in the sum of squares is at most tol)
    obtained parameters:
      par[0] =      0.12987
      par[1] =        -0.05
      par[2] =      1.03052
    obtained norm:
          0.450685
    fitting data as follows:
      t[ 0]=  -4 y=  16.6 fit=   16.8182 residue=   -0.218182
      t[ 1]=  -3 y=   9.9 fit=   9.55455 residue=    0.345455
      t[ 2]=  -2 y=   4.4 fit=   4.35195 residue=    0.048052
      t[ 3]=  -1 y=   1.1 fit=   1.21039 residue=    -0.11039
      t[ 4]=   0 y=     0 fit=   0.12987 residue=    -0.12987
      t[ 5]=   1 y=   1.1 fit=   1.11039 residue=  -0.0103896
      t[ 6]=   2 y=   4.2 fit=   4.15195 residue=    0.048052
      t[ 7]=   3 y=   9.3 fit=   9.25455 residue=   0.0454545
      t[ 8]=   4 y=  16.4 fit=   16.4182 residue=  -0.0181818

To obtain progress messages during the fit, one needs to modify the control record princon. This must be done after initialization of princon, and before calling lmcurve:

    lm_princon_struct princon = lm_princon_std;
    princon.flags = 3; /* a decent level of verbosity */
    /* ... */
    lmcurve( n, par, m, t, y, f, &control, &princon, &status );    

To obtain less verbose progress messages, switch to expert mode:

    princon.form = 1;