Figure 9.14:

Confidence intervals with known (solid line) and unknown (dashed line) error variance.

Code for Figure 9.14

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
 98
 99
100
101
102
103
104
105
% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License as
% published by the Free Software Foundation; either version 2, or (at
% your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; see the file COPYING.  If not, write to
% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
% MA 02111-1307, USA.

%
% arhenius.m
%
p     = 2;
ndata = 10;
chisq = 5.99;  % chi square value for 95% confidence level with 2 parameters
Fstat = 4.46; % F statistic for 95 % confidence level with 2 parameters and
				% 10 data points (or 10-2=8 degrees of
				% freedcom)
lnk0  = 1;
E     = 100;
Tmin  = 300;
Tmax  = 500;
Tmeas = linspace(Tmin,Tmax,ndata)';
X     = [ones(ndata,1) -1./Tmeas];
lnk   = X*[lnk0; E];
k     = exp(lnk);
measvar= 1e-3;
measstddev = sqrt(measvar);
%
% set seed for 'reproducible' random numbers
%
rng(0);
%
%just a few noisy points
%
nexpts = 1;
clear lnkmeas;
for i = 1:nexpts lnkmeas(:,i) = lnk + measstddev*randn(ndata,1); end
%
%center the data first
%
Tcenter = -1./Tmeas + 1/mean(Tmeas);
Xcenter=[ones(ndata,1) Tcenter];
thetacenter = inv(Xcenter'*Xcenter)*Xcenter' * lnkmeas;

thetacentertr = thetacenter';

residual = lnkmeas-Xcenter*thetacenter;
samplevar = residual'*residual/(ndata-p);

npts = 181;
amat = Xcenter'*Xcenter;
level = chisq*measvar;
[x, y, major, minor, bbox] = ellipse (amat, level, npts);
%
% shift the ellipse's center to the estimated parameters
%
lnkest = thetacenter(1);
Eest   = thetacenter(2);
x = x+lnkest;
y = y+Eest;
minor(:,1)=minor(:,1)+lnkest;
minor(:,2)=minor(:,2)+Eest;
major(:,1)=major(:,1)+lnkest;
major(:,2)=major(:,2)+Eest;
bbox(:,1)=bbox(:,1)+lnkest;
bbox(:,2)=bbox(:,2)+Eest;
bbox1 = bbox;

outline1 = [x, y];

level = Fstat*p*samplevar;
[x, y, major, minor, bbox] = ellipse (amat, level, npts);
%
% shift the ellipse's center to the estimated parameters
%
x = x+lnkest;
y = y+Eest;
minor(:,1)=minor(:,1)+lnkest;
minor(:,2)=minor(:,2)+Eest;
major(:,1)=major(:,1)+lnkest;
major(:,2)=major(:,2)+Eest;
bbox(:,1)=bbox(:,1)+lnkest;
bbox(:,2)=bbox(:,2)+Eest;
bbox2 = bbox;

outline2 = [x, y];

save parestellipse.dat thetacentertr bbox1 outline1 bbox2 outline2;

if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
plot (thetacentertr(:,1), thetacentertr(:,2), '+', ...
       bbox1(:,1), bbox1(:,2), outline1(:,1), outline1(:,2), ...
       bbox2(:,1), bbox2(:,2), outline2(:,1), outline2(:,2));
axis ([0.6, 1.0, 80, 220]);
% TITLE
end % PLOTTING