Figure 9.11:

Reducing parameter correlation by centering the data.

Code for Figure 9.11

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
106
% 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
%
chisq = 5.99;  % chi square value for 95% confidence level with 2 parameters
lnk0  = 1;
E     = 100;
ndata = 10;
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 = 500;
clear lnkmeas;
for i = 1:nexpts lnkmeas(:,i) = lnk + measstddev*randn(ndata,1); end
%for i = 1:nexpts kmeas(:,i) = k + measstddev*randn(ndata,1); end
%lnkmeas = log(kmeas);
theta = (inv(X'*X)*X' * lnkmeas)';

%save_precision=5;
%save -ascii arrhenius_center.dat rest;
%
%center the data first
%
Tcenter = -1./Tmeas + 1/mean(Tmeas);
Xcenter=[ones(ndata,1) Tcenter];
thetacenter = (inv(Xcenter'*Xcenter)*Xcenter' * lnkmeas)';

%save -ascii arrhenius_center-theta.dat theta;
%save -ascii arrhenius_center-thetacenter.dat thetacenter;

npts = 181;
amat = X'*X/measvar;
level = chisq;
[x, y, major, minor, bbox] = ellipse (amat, level, npts);
%
% shift the ellipse's center to the optimal parameters
%
x = x+lnk0;
y = y+E;
minor(:,1)=minor(:,1)+lnk0;
minor(:,2)=minor(:,2)+E;
major(:,1)=major(:,1)+lnk0;
major(:,2)=major(:,2)+E;
bbox(:,1)=bbox(:,1)+lnk0;
bbox(:,2)=bbox(:,2)+E;
bbox1 = bbox;
%save -ascii arrhenius_center-bbox-1.dat bbox1;
outline1 = [x, y];
%save -ascii arrhenius_center-outline-1.dat outline1;

amat = Xcenter'*Xcenter/measvar;
[x, y, major, minor, bbox] = ellipse (amat, level, npts);
%
% shift the ellipse's center to the optimal parameters
%
lnkmean = lnk0 - E/mean(Tmeas);
x = x+lnkmean;
y = y+E;
minor(:,1)=minor(:,1)+lnkmean;
minor(:,2)=minor(:,2)+E;
major(:,1)=major(:,1)+lnkmean;
major(:,2)=major(:,2)+E;
bbox(:,1)=bbox(:,1)+lnkmean;
bbox(:,2)=bbox(:,2)+E;
bbox2 = bbox;
%save -ascii arrhenius_center-bbox-2.dat bbox2;
outline2 = [x, y];
%save -ascii arrhenius_center-outline-2.dat outline2;

save arrhenius_center.dat theta thetacenter bbox1 outline1 bbox2 outline2;

if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
plot (theta(:,1), theta(:,2), '+', ...
       thetacenter(:,1), thetacenter(:,2), 'x', ...
       bbox1(:,1), bbox1(:,2), outline1(:,1), outline1(:,2), ...
       bbox2(:,1), bbox2(:,2), outline2(:,1), outline2(:,2));
% TITLE
end % PLOTTING