Non-linear least square fitting of time-series exploits the function ao/curvefit.

Non-linear least square fitting of time series - Exercise

During this exercise we will:

  1. Load time series data
  2. Fit data with ao/curvefit
  3. Check results
  4. Refine the fit
  5. Redo a full parametrized fit

Let's open a new editor window and load test data.

    a = ao(plist('filename', 'topic5/T5_Ex05_TestNoise.xml'));
    a.setName;
    iplot(a)

As can be seen this is a chirped sine wave with some noise.

We could now try to doing the fit. The first parameter to pass to curvefit is a fit model. In this case we assume that we are dealing with a linearly chirped sine wave. Then we need to specify a starting guess for the function parameters and, if we want, we can set lower and upper bounds for the fit parameters. Further parameters (unchanged by the fit) could be added as a cell-array in the parameter ADDP. In the present case we added a bias of 5 for our data.

The output of ao/curvefit is a cdata analysis objects containing fit parameters.

    plfit = plist('Function', 'ADDP{1} + P(1).*sin(2.*pi.*(P(2) + P(3).*Xdata).*Xdata + P(4))', ...
              'P0', [4 3e-5 5e-6 0.5], ...
              'LB', [1 0 0 -pi], ...
              'UB', [5 1 1 pi],...
              'ADDP', {5});
    params = curvefit(a, plfit);

Once the fit is done. We can evaluate our model to check fit results.

    pleval = plist('Function', 'ADDP{1} + P(1).*sin(2.*pi.*(P(2) + P(3).*Xdata).*Xdata + P(4))', ...
              'Xdata', a, ...
              'dtype', 'tsdata', ...
              'ADDP', {5});
    b = evaluateModel(params, pleval);
    b.setYunits(a.yunits);
    b.setXunits(a.xunits);
    b.setName;    
    iplot(a,b)

As you can see, the fit is not accurate. One of the great problems of non-linear least square methods is that they easily find a local minimum of the chi square function and stop there without finding a global minimum. So, let's try to look at the data and to refine step by step the fit.

Looking at data we could try to fix some parameters like the amplitude and phase. The y starting level is near 0.4 + bias corresponding to a phase angle around 0.4 rad; the amplitude is close to 3. Now run the fit again with a reduced set of parameters.

    % Do the fit again
    plfit = plist('Function', 'ADDP{1} + 3.*sin(2.*pi.*(P(1) + P(2).*Xdata).*Xdata + 0.4)', ...
                'P0', [7e-4 9e-6], ...
                'LB', [1e-7 1e-7], ...
                'UB', [1 1e-4],...
                'ADDP', {5});
    
    params = curvefit(a, plfit);
    
    % Evaluate the model
    pleval = plist('Function', 'ADDP{1} + 3.*sin(2.*pi.*(P(1) + P(2).*Xdata).*Xdata + 0.4)', ...
                  'Xdata', a, ...
                  'dtype', 'tsdata', ...
                  'ADDP', {5});
    b = evaluateModel(params, pleval);
    b.setYunits(a.yunits);
    b.setXunits(a.xunits);
    b.setName;
    
    iplot(a,b)

The fit now looks like better...

Actual parameters contained in params are now a good starting guess for the wave frequancy and chirp parameter.
We could try a new fit with the complete set of parameters...

    % Do the fit again
    plfit = plist('Function', 'ADDP{1} + P(1).*sin(2.*pi.*(P(2) + P(3).*Xdata).*Xdata + P(4))', ...
                  'P0', [3 5e-5 1e-5 0.4], ...
                  'LB', [2.8 1e-5 1e-6 0.2], ...
                  'UB', [3.2 5e-4 5e-4 0.5],...
                  'ADDP', {5});
    
    params = curvefit(a, plfit);
    
    % Evaluate the model
    pleval = plist('Function', 'ADDP{1} + P(1).*sin(2.*pi.*(P(2) + P(3).*Xdata).*Xdata + P(4))', ...
                  'Xdata', a, ...
                  'dtype', 'tsdata', ...
                  'ADDP', {5});
    b = evaluateModel(params, pleval);
    b.setYunits(a.yunits);
    b.setXunits(a.xunits);
    b.setName;
    
    iplot(a,b)

Let us compare fit results with nominal parameters.
Data were generated with the following set of parameters:

    ADDP = 5
    P(1) = 3
    P(2) = 1e-4
    P(3) = 1e-5
    P(4) = 0.3

Fitted parameters are instead:

    ADDP = 5
    P(1) = 2.993
    P(2) = 0.000121
    P(3) = 9.983e-006
    P(4) = 0.278

The correlation matrix of the parameters is stored in the procinfo (processing information) field of the AO. This field is a plist and is used to additional information that can be returned from algorithms. In this case it returns the following plist:
Key Description

COR

The covariance matrix for the parameters.

DPARAMS

A structure of additional information about the parameters, for example, upper and lower values at the 1-sigma significance level.

GOF

A 'goodness of fit' structure which contains fields like, the Chi^2 value and the degrees of freedom.

To extract the covariance matrix from the procinfo you can do

    params.procinfo.find('cor')
Columns 1 through 3

                         1        0.0220543553134523       0.00840698749447142
        0.0220543553134523                         1        -0.963274881180157
       0.00840698749447142        -0.963274881180157                         1
       -0.0911417933676055        -0.833580057704702         0.692767145487321

  Column 4

       -0.0911417933676055
        -0.833580057704702
         0.692767145487321
                         1
                       

Not so bad!