Non-linear least square fitting of time series exploits the function ao/curvefit.
During this exercise we will:
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.
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.
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...
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)
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!