Profile Estimation Experiments - the Rossler Equations
This page provides code to run profiled parameter estimation on data generated by the Rossler equations. It follows the same format as FhNEx.html and has commentary has therefore been kept to a minimum.
The Rossler equations are given by
Contents
- RHS Functions
- Various parameters
- Observation times
- Create trajectories
- Set up observations
- Fitting parameters
- Profiling optimisation control
- Setting up functional data objects
- Smooth the data
- Re-smoothing with model-based penalty
- Perform the Profiled Estimation
- Plot Smooth with Profile-Estimated Parameters
- Comparison with Smooth Using True Parameters
- Squared Error Performance
- Calculate Sample Information and Variance-Covariance Matrices
RHS Functions
odefn = @rossfunode; % Function for ODE solver fn.fn = @rossfun; % RHS function fn.dfdx = @rossdfdx; % Derivative wrt inputs fn.dfdp = @rossdfdp; % Derviative wrt parameters fn.d2fdx2 = @rossd2fdx2; % Hessian wrt inputs fn.d2fdxdp = @rossd2fdxdp; % Hessian wrt inputs and parameters fn.d2fdp2 = @rossd2fdp2; % Hessian wrt parameters. fn.d3fdx3 = @rossd3fdx3; % Third derivative wrt inputs. fn.d3fdx2dp = @rossd3fdx2dp; % Third derivative wrt intputs, inputs and pars. fn.d3fdxdp2 = @rossd3fdxdp2; % Third derivative wrt inputs, pars and pars.
Various parameters
y0 = [1.13293; -1.74953; 0.02207]; % Initial conditions pars = [0.2; 0.2; 3]; % Parameters sigma = 0.5; % Noise Level jitter = 0.2; % Perturbation for starting startpars = pars + jitter*randn(length(pars),1); % parameter estimates disp(['Initial par. values: ',num2str(startpars')])
Initial par. values: 0.13571 -0.16501 3.2038
Observation times
tspan = 0:0.05:20; % Observation times obs_pts{1} = 1:length(tspan); % Which components are observed at obs_pts{2} = 1:length(tspan); % which observation times. obs_pts{3} = 1:length(tspan); tfine = 0:0.05:20; % Times to plot solutions
Create trajectories
odeopts = odeset('RelTol',1e-13);
[full_time,full_path] = ode45(odefn,tspan,y0,odeopts,pars);
[plot_time,plot_path] = ode45(odefn,tfine,y0,odeopts,pars);
Set up observations
Tcell = cell(1,size(full_path,2)); path = Tcell; for i = 1:length(obs_pts) Tcell{i} = full_time(obs_pts{i}); path{i} = full_path(obs_pts{i},i); end % add noise Ycell = path; for i = 1:length(path) Ycell{i} = path{i} + sigma*randn(size(path{i})); end % and set wts wts = []; if isempty(wts) % estimate wts if not given for i = 1:length(Ycell) if ~isempty(Ycell{i}) wts(i) = 1./sqrt(var(Ycell{i})); else wts(i) = 1; end end end
Fitting parameters
lambda = 1e4; % Smoothing for model-based penalty lambda = lambda*wts; lambda0 = 1; % Smoothing for 1st-derivative penalty nknots = 401; % Number of knots to use. nquad = 5; % No. between-knots quadrature points. norder = 3; % Order of B-spline approximation
Profiling optimisation control
lsopts_out = optimset('DerivativeCheck','off','Jacobian','on',... 'Display','iter','MaxIter',1000,'TolFun',1e-8,'TolX',1e-10); % Other observed optimiation control lsopts_other = optimset('DerivativeCheck','off','Jacobian','on',... 'Display','iter','MaxIter',1000,'TolFun',1e-14,'TolX',1e-14,... 'JacobMult',@SparseJMfun); % Optimiation control within profiling lsopts_in = optimset('DerivativeCheck','off','Jacobian','on',... 'Display','off','MaxIter',1000,'TolFun',1e-14,'TolX',1e-14,... 'JacobMult',@SparseJMfun);
Setting up functional data objects
% set up knots range = [min(full_time),max(full_time)]; % range of observations knots = linspace(range(1),range(2),nknots); % set up bases basis_cell = cell(1,length(path)); % Create cell arrays. Lfd_cell = cell(1,length(path)); nbasis = nknots+norder-2; quadvals = [knots' ones(nknots,1)/nknots]; bbasis = MakeBasis(range,nbasis,norder,knots,quadvals,1); basis_cell(:) = {bbasis}; Lfd_cell(:) = {fdPar(bbasis,1,lambda0)};
Smooth the data
DEfd = smoothfd_cell(Ycell,Tcell,Lfd_cell); coefs = getcellcoefs(DEfd); figure(1) devals = eval_fdcell(tfine,DEfd,0); for i = 1:length(path) subplot(length(path),1,i); plot(Tcell{i},Ycell{i},'b.'); hold on; plot(tfine,devals{i},'r','LineWidth',2); if i==1 ylabel('\fontsize{13} X') title(['\fontsize{13} Raw data (.), ', ... 'and smooth fit (r-)']) elseif i==2 ylabel('\fontsize{13} Y') else xlabel('\fontsize{13} t') ylabel('\fontsize{13} Z') end end
Re-smoothing with model-based penalty
% Call the Gauss-Newton solver [newcoefs,resnorm2] = lsqnonlin(@SplineCoefErr,coefs,[],[],lsopts_other,... basis_cell,Ycell,Tcell,wts,lambda,fn,[],startpars); tDEfd = Make_fdcell(newcoefs,basis_cell); % Plot results along with exact solution figure(2) devals = eval_fdcell(tfine,tDEfd,0); for i = 1:length(path) subplot(length(path),1,i); plot(tfine,devals{i},'r','LineWidth',2); hold on; plot(Tcell{i},Ycell{i},'b.'); plot(plot_time,plot_path(:,i),'c'); hold off if i==1 ylabel('\fontsize{13} X') title(['\fontsize{13} Raw data (.), ', ... 'exact solution (r-) and true path (g-)']) elseif i==2 ylabel('\fontsize{13} Y') else xlabel('\fontsize{13} t') ylabel('\fontsize{13} Z') end end
Norm of First-order
Iteration Func-count f(x) step optimality CG-iterations
0 1 10744.1 1.28e+003
1 2 1226.46 4.11518 91.6 17
2 3 540.306 10 49.2 138
3 4 420.263 7.88586 23.4 218
4 5 400.472 3.37831 8.33 172
5 6 396.714 2.83926 3.73 180
6 7 396.623 2.16356 6.16 245
7 8 395.193 0.540891 30 229
8 9 394.779 1.13383 1.28 193
9 10 394.752 0.959113 0.787 298
10 11 394.496 0.239778 11.6 295
11 12 394.471 0.494898 0.379 210
12 13 394.413 0.119889 3.44 323
13 14 394.41 0.196755 0.123 227
14 15 394.399 0.0491888 0.531 351
15 16 394.399 0.0994648 0.076 197
16 17 394.396 0.0245944 0.357 344
17 18 394.396 0.0381179 0.0234 191
18 19 394.396 0.00952948 0.0397 341
19 20 394.396 0.0180631 0.013 235
20 21 394.396 0.00451578 0.0893 348
21 22 394.396 0.0071157 0.00513 183
22 23 394.396 0.00177893 0.00678 313
23 24 394.396 0.00325136 0.00248 232
24 25 394.396 0.00081284 0.00726 348
25 26 394.396 0.00131989 0.000897 187
26 27 394.396 0.000329973 0.00126 322
27 28 394.396 0.000593669 0.000426 226
28 29 394.396 0.000148417 0.00268 349
29 30 394.396 0.000243717 0.000168 183
30 31 394.396 6.09293e-005 0.000222 305
31 32 394.396 0.00010617 7.85e-005 199
32 33 394.396 2.65424e-005 0.000241 349
33 34 394.396 4.28563e-005 3.07e-005 195
34 35 394.396 1.07141e-005 5.57e-005 316
35 36 394.396 1.85588e-005 1.39e-005 217
36 37 394.396 4.63969e-006 2.44e-005 352
37 38 394.396 7.61147e-006 5.26e-006 184
38 39 394.396 1.90287e-006 1.05e-005 303
39 40 394.396 3.2754e-006 2.43e-006 210
40 41 394.396 8.1885e-007 1.29e-005 349
Optimization terminated: relative function value
changing by less than OPTIONS.TolFun.
Perform the Profiled Estimation
[newpars,newDEfd_cell] = Profile_GausNewt(startpars,lsopts_out,DEfd,fn,... lambda,Ycell,Tcell,wts,[],lsopts_in); disp(['New parameter values: ',num2str(newpars')]);
Iteration steps Residual Improvement Grad-norm parameters
1 1 186.049 0.296138 17.7 0.19446 0.25024 2.8813
2 1 182.233 0.0205111 0.816 0.19443 0.17858 2.9136
3 1 182.23 1.42551e-005 0.0229 0.19334 0.17796 2.9105
4 1 182.23 1.80858e-008 0.000892 0.19338 0.17798 2.9107
5 9 182.23 -2.4344e-010 0.000889 0.19338 0.17798 2.9107
New parameter values: 0.19338 0.17798 2.9107
Plot Smooth with Profile-Estimated Parameters
devals = eval_fdcell(tfine,newDEfd_cell,0); figure(3) for i = 1:length(path) subplot(length(path),1,i) plot(tfine,devals{i},'r','LineWidth',2); hold on; plot(Tcell{i},Ycell{i},'b.'); plot(plot_time,plot_path(:,i),'c'); hold off if i==1 ylabel('\fontsize{13} X') title(['\fontsize{13} Raw data (.), ', ... 'profiled solution (r-) and true path (g-)']) elseif i==2 ylabel('\fontsize{13} Y') else xlabel('\fontsize{13} t') ylabel('\fontsize{13} Z') end end
Comparison with Smooth Using True Parameters
coefs = getcellcoefs(DEfd); % Starting coefficient estimates [truecoefs,resnorm4] = lsqnonlin(@SplineCoefErr,coefs,[],[],... lsopts_other,basis_cell,Ycell,Tcell,wts,lambda,fn,[],pars); trueDEfd_cell = Make_fdcell(truecoefs,basis_cell); figure(4) devals = eval_fdcell(tfine,trueDEfd_cell,0); for i = 1:length(path) subplot(length(path),1,i) plot(plot_time,plot_path(:,i),'c') plot(tfine,devals{i},'r','LineWidth',2); hold on; plot(plot_time,plot_path(:,i),'c'); plot(Tcell{i},Ycell{i},'b.'); hold off; if i==1 ylabel('\fontsize{13} X') title(['\fontsize{13} Raw data (.), ', ... 'exact solution (r-) and true path (g-)']) elseif i==2 ylabel('\fontsize{13} Y') else xlabel('\fontsize{13} t') ylabel('\fontsize{13} Z') end end
Norm of First-order
Iteration Func-count f(x) step optimality CG-iterations
0 1 8654.04 1.35e+003
1 2 421.891 3.55023 59 16
2 3 200.788 2.92019 9.9 60
3 4 184.62 4.12092 7.66 277
4 5 183.018 0.095115 0.831 49
5 6 182.927 0.300157 0.202 254
6 7 182.926 0.0313522 0.00446 310
7 8 182.926 0.00287643 0.000353 291
8 9 182.926 0.000414949 3.57e-005 342
9 10 182.926 1.72255e-005 3.36e-006 253
10 11 182.926 3.05275e-006 2.98e-007 311
11 12 182.926 2.05804e-007 2.98e-007 194
12 13 182.926 5.14511e-008 2.98e-007 0
13 14 182.926 1.28628e-008 2.98e-007 0
14 15 182.926 3.21569e-009 2.98e-007 0
15 16 182.926 8.03923e-010 6.52e-007 0
Optimization terminated: relative function value
changing by less than OPTIONS.TolFun.
Squared Error Performance
% Squared error for estimated parameters newpreds = eval_fdcell(Tcell,newDEfd_cell,0); new_err = cell(length(newpreds)); for i = 1:length(path) new_err{i} = wts(i)*(newpreds{i} - Ycell{i}).^2; end new_err = mean(cell2mat(new_err)); % Squared error for true parameters truepreds = eval_fdcell(Tcell,trueDEfd_cell,0); true_err = cell(length(truepreds)); for i = 1:length(path) true_err{i} = wts(i)*(truepreds{i} - Ycell{i}).^2; end true_err = mean(cell2mat(true_err)); % print out a comparsion disp(['Estimated sqrd error: ',num2str(new_err)]) disp(['True sqrd error: ',num2str(true_err)]);
Estimated sqrd error: 0.15148 True sqrd error: 0.15174
Calculate Sample Information and Variance-Covariance Matrices
% Hessian of squared error with respect to parameters d2Jdp2 = make_d2jdp2(newDEfd_cell,fn,Ycell,Tcell,lambda,newpars,[],wts); % Second derivatives with respect to parameters and observations d2JdpdY = make_d2jdpdy(newDEfd_cell,fn,Ycell,Tcell,lambda,newpars,[],wts); % Resulting derivative of parameters with respect to observations dpdY = -d2Jdp2\d2JdpdY; % Variance of observations: S = make_sigma(DEfd,Tcell,Ycell,0); % Resulting parameter covariance matrix: Cov = dpdY*S*dpdY'; % Standard errors StdDev = sqrt(diag(Cov)); % Correlations Corr = Cov./(StdDev*StdDev'); % Display these results disp('Approximate covariance matrix for parameters:') disp(num2str(Cov)) disp('Approximate standard errors of parameters:') disp(num2str(StdDev')) disp('Approximate correlation matrix for parameters:') disp(num2str(Corr))
Approximate covariance matrix for parameters:
2.5357e-005 2.6026e-005 7.9651e-005
2.6026e-005 0.0003967 0.0012325
7.9651e-005 0.0012325 0.0050578
Approximate standard errors of parameters:
0.0050355 0.019917 0.071119
Approximate correlation matrix for parameters:
1 0.2595 0.22241
0.2595 1 0.87013
0.22241 0.87013 1