/usr/share/octave/packages/optim-1.3.0/polyconf.m is in octave-optim 1.3.0-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
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 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 | ## Author: Paul Kienzle <pkienzle@gmail.com>
## This program is granted to the public domain.
## [y,dy] = polyconf(p,x,s)
##
## Produce prediction intervals for the fitted y. The vector p
## and structure s are returned from polyfit or wpolyfit. The
## x values are where you want to compute the prediction interval.
##
## polyconf(...,['ci'|'pi'])
##
## Produce a confidence interval (range of likely values for the
## mean at x) or a prediction interval (range of likely values
## seen when measuring at x). The prediction interval tells
## you the width of the distribution at x. This should be the same
## regardless of the number of measurements you have for the value
## at x. The confidence interval tells you how well you know the
## mean at x. It should get smaller as you increase the number of
## measurements. Error bars in the physical sciences usually show
## a 1-alpha confidence value of erfc(1/sqrt(2)), representing
## one standandard deviation of uncertainty in the mean.
##
## polyconf(...,1-alpha)
##
## Control the width of the interval. If asking for the prediction
## interval 'pi', the default is .05 for the 95% prediction interval.
## If asking for the confidence interval 'ci', the default is
## erfc(1/sqrt(2)) for a one standard deviation confidence interval.
##
## Example:
## [p,s] = polyfit(x,y,1);
## xf = linspace(x(1),x(end),150);
## [yf,dyf] = polyconf(p,xf,s,'ci');
## plot(xf,yf,'g-;fit;',xf,yf+dyf,'g.;;',xf,yf-dyf,'g.;;',x,y,'xr;data;');
## plot(x,y-polyval(p,x),';residuals;',xf,dyf,'g-;;',xf,-dyf,'g-;;');
function [y,dy] = polyconf(p,x,varargin)
alpha = s = [];
typestr = 'pi';
for i=1:length(varargin)
v = varargin{i};
if isstruct(v), s = v;
elseif ischar(v), typestr = v;
elseif isscalar(v), alpha = v;
else s = [];
end
end
if (nargout>1 && (isempty(s)||nargin<3)) || nargin < 2
print_usage;
end
if isempty(s)
y = polyval(p,x);
else
## For a polynomial fit, x is the set of powers ( x^n ; ... ; 1 ).
n=length(p)-1;
k=length(x(:));
if columns(s.R) == n, ## fit through origin
A = (x(:) * ones (1, n)) .^ (ones (k, 1) * (n:-1:1));
p = p(1:n);
else
A = (x(:) * ones (1, n+1)) .^ (ones (k, 1) * (n:-1:0));
endif
y = dy = x;
[y(:),dy(:)] = confidence(A,p,s,alpha,typestr);
end
end
%!test
%! # data from Hocking, RR, "Methods and Applications of Linear Models"
%! temperature=[40;40;40;45;45;45;50;50;50;55;55;55;60;60;60;65;65;65];
%! strength=[66.3;64.84;64.36;69.70;66.26;72.06;73.23;71.4;68.85;75.78;72.57;76.64;78.87;77.37;75.94;78.82;77.13;77.09];
%! [p,s] = polyfit(temperature,strength,1);
%! [y,dy] = polyconf(p,40,s,0.05,'ci');
%! assert([y,dy],[66.15396825396826,1.71702862681486],200*eps);
%! [y,dy] = polyconf(p,40,s,0.05,'pi');
%! assert(dy,4.45345484470743,200*eps);
## [y,dy] = confidence(A,p,s)
##
## Produce prediction intervals for the fitted y. The vector p
## and structure s are returned from wsolve. The matrix A is
## the set of observation values at which to evaluate the
## confidence interval.
##
## confidence(...,['ci'|'pi'])
##
## Produce a confidence interval (range of likely values for the
## mean at x) or a prediction interval (range of likely values
## seen when measuring at x). The prediction interval tells
## you the width of the distribution at x. This should be the same
## regardless of the number of measurements you have for the value
## at x. The confidence interval tells you how well you know the
## mean at x. It should get smaller as you increase the number of
## measurements. Error bars in the physical sciences usually show
## a 1-alpha confidence value of erfc(1/sqrt(2)), representing
## one standandard deviation of uncertainty in the mean.
##
## confidence(...,1-alpha)
##
## Control the width of the interval. If asking for the prediction
## interval 'pi', the default is .05 for the 95% prediction interval.
## If asking for the confidence interval 'ci', the default is
## erfc(1/sqrt(2)) for a one standard deviation confidence interval.
##
## Confidence intervals for linear system are given by:
## x' p +/- sqrt( Finv(1-a,1,df) var(x' p) )
## where for confidence intervals,
## var(x' p) = sigma^2 (x' inv(A'A) x)
## and for prediction intervals,
## var(x' p) = sigma^2 (1 + x' inv(A'A) x)
##
## Rather than A'A we have R from the QR decomposition of A, but
## R'R equals A'A. Note that R is not upper triangular since we
## have already multiplied it by the permutation matrix, but it
## is invertible. Rather than forming the product R'R which is
## ill-conditioned, we can rewrite x' inv(A'A) x as the equivalent
## x' inv(R) inv(R') x = t t', for t = x' inv(R)
## Since x is a vector, t t' is the inner product sumsq(t).
## Note that LAPACK allows us to do this simultaneously for many
## different x using sqrt(sumsq(X/R,2)), with each x on a different row.
##
## Note: sqrt(F(1-a;1,df)) = T(1-a/2;df)
##
## For non-linear systems, use x = dy/dp and ignore the y output.
function [y,dy] = confidence(A,p,S,alpha,typestr)
if nargin < 4, alpha = []; end
if nargin < 5, typestr = 'ci'; end
y = A*p(:);
switch typestr,
case 'ci', pred = 0; default_alpha=erfc(1/sqrt(2));
case 'pi', pred = 1; default_alpha=0.05;
otherwise, error("use 'ci' or 'pi' for interval type");
end
if isempty(alpha), alpha = default_alpha; end
s = tinv(1-alpha/2,S.df)*S.normr/sqrt(S.df);
dy = s*sqrt(pred+sumsq(A/S.R,2));
end
|