Non-linear least square fitting of time series


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 fit parameter. Further papameters (unchanged by the fit) could be added as ADDP . In the present case we added a bias of 5 for our data.
ao/curvefit output 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 global minimum. So, let's try to look at data and to refine step by step the fit.

Looking at data we could try to fix some parameters like wave amplitude and phase. y starting level is near 0.4 + bias corresponding to a phase angle around 0.4 rad wave amplitude is instead easily set to 3. Run againg the fit with a reduced set of parameters.

  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);

  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...

  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);

  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.99300840665851
  P(2) = 0.000120496292686661
  P(3) = 9.98313996200842e-006
  P(4) = 0.27762605166717

Not so bad!




©LTP Team