Figure 9.22:

Monte Carlo evaluation of confidence intervals.

Code for Figure 9.22

Text of the GNU GPL.

main.m


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
%
% jbr,  4/8/18
%


ymodel = struct;
model.x = {'ca'};
model.p = {'k', 'n'};

model.d = {'m_ca'};

model.ode = @(t,x,p) {-p.k*x.ca^p.n};

model.lsq = @(t,x,p) {x.ca-x.m_ca};

tfinal = 5;
nts    = 100;
tout = linspace(0,tfinal,nts);
model.tout = tout;

pe = paresto(model);

kac   = 0.5;
ca0ac = 2;
nac   = 2.5;
thetaac = [kac; nac; ca0ac];

x_ac = ca0ac;
p_ac = [kac; nac];

y_ac = pe.simulate(0, x_ac, p_ac);

% add measurement noise
measvar = 1e-2;
measstddev = sqrt(measvar);
rng(0);
noise = measstddev*randn(1,nts);

y_noisy = y_ac + noise;


% Initial guess, upper and lower bounds for the estimated parameters
theta0 = thetaac;
small = 1e-3;
large  = 5;
np = numel(theta0);
lbtheta = small*ones(np,1);
ubtheta = large*ones(np,1);


% loop here to create nmonte data sets and estimate parameters for each
page_output_immediately(1);
more off

ndata = nts;
nmonte = 10; % change to 500 for figure in text
alph(1:nmonte) = NaN;

for i = 1: nmonte
  %create  noisy data
  noise = measstddev*randn(1, nts);
  y_noisy = y_ac + noise;
  %estimate parameters
  [est, y, p] = pe.optimize(y_noisy, theta0, lbtheta, ubtheta);
  est.d2f_dtheta2;
  i
  theta_conf = pe.confidence(est, 1:np, 0.95);
  [est.theta theta_conf]
  % store the estimate, contour distance of estimate from true parameters,
  % alpha value of this contour
  obj = est.f;
  theta = est.theta;
  H = est.d2f_dtheta2;
  samplevar = obj/(ndata - np);
  thetamont(:,i) = theta;
  diff(:,i) = est.theta - thetaac;
  cont = diff(:,i)'*H*diff(:,i)/(2*np*samplevar);
  alph(i) = fcdf(cont, np, ndata-np);
end

% post processing for confidence region plot
amin = 0.0;
amax = 0.99995;
nas  = 100;
alphavec=linspace(amin,amax,nas)';
expected = alphavec*nmonte;
actual(1:nas) = NaN;
for i = 1: nas
  actual(i) = sum(alph <= alphavec(i));
end

table = [alphavec expected actual(:)];
save nthmonte.dat table;
if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
plot (alphavec, expected, alphavec, actual, '+');
% TITLE
end % PLOTTING