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.
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.
% 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...
% 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)
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!