Non-linear least square fitting of time-series exploits the function ao/curvefit.
During this exercise we will:
Let us 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 the global minimum. So, let us 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!