Figure 9.14 (page 523):

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

Code for Figure 9.14

Text of the GNU GPL.

main.m

%% 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
%
randn('seed',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];

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 ('Figure 9.14')

../../util/common/ellipse.m

%% Copyright (C) 2001, James B. Rawlings and John W. Eaton
%%
%% 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.

%% [x, y, major, minor, bbox] = ellipse (amat, level, n, shift)
%%
%% Given a 2x2 matrix, generate ellipse data for plotting.  The
%% arguments N and SHIFT are optional.  If N is an empty matrix, a
%% default value of 100 is used.

function [x, y, major, minor, bbox] = ellipse (amat, level, n, shift)

  if (nargin < 3)
    n = 100;
  end

  if (isempty (n))
    n = 100;
  end

  if (nargin < 4)
    shift = [0, 0];
  end

  ss = size (shift);

  if (any (ss ~= [1, 2]))
    if (ss == [2, 1])
      shift = shift';
    else
      error ('shift must be a 2-element row vector');
    end
  end

  if (nargin > 1)

    [v, l] = eig (amat / level);

    dl = diag(l);
    if (any (imag (dl)) || any (dl <= 0))
      error ('ellipse: amat must be positive definite');
    end

    %% Generate contour data.

    a = 1 / sqrt (l(1,1));
    b = 1 / sqrt (l(2,2));

    t = linspace (0, 2*pi, n)';

    xt = a * cos (t);
    yt = b * sin (t);

    %% Rotate the contours.

    ra = atan2 (v(2,1), v(1,1));

    cos_ra = cos (ra);
    sin_ra = sin (ra);

    x = xt * cos_ra - yt * sin_ra + shift(1);
    y = xt * sin_ra + yt * cos_ra + shift(2);

    %% Endpoints of the major and minor axes.

    minor = (v * diag ([a, b]))';
    major = minor;

    major(2,:) = -major(1,:);
    minor(1,:) = -minor(2,:);

    t = [1; 1] * shift;

    major = major + t;
    minor = minor + t;

    %% Bounding box for the ellipse using magic formula.

    ainv = inv (amat);
    xbox = sqrt (level * ainv(1,1));
    ybox = sqrt (level * ainv(2,2));

    bbox = [xbox ybox; xbox -ybox; -xbox -ybox; -xbox ybox; xbox ybox];

    t = [1; 1; 1; 1; 1] * shift;
    bbox = bbox + t;

  else
    error ('usage: ellipse (amat, level, n, shift)');
  end

%% endfunction