octave_packages/optim-1.2.0/polyconf.m

   1 ## Author: Paul Kienzle <pkienzle@gmail.com>
   2 ## This program is granted to the public domain.
   3
   4 ## [y,dy] = polyconf(p,x,s)
   5 ##
   6 ##   Produce prediction intervals for the fitted y. The vector p
   7 ##   and structure s are returned from polyfit or wpolyfit. The
   8 ##   x values are where you want to compute the prediction interval.
   9 ##
  10 ## polyconf(...,['ci'|'pi'])
  11 ##
  12 ##   Produce a confidence interval (range of likely values for the
  13 ##   mean at x) or a prediction interval (range of likely values
  14 ##   seen when measuring at x).  The prediction interval tells
  15 ##   you the width of the distribution at x.  This should be the same
  16 ##   regardless of the number of measurements you have for the value
  17 ##   at x.  The confidence interval tells you how well you know the
  18 ##   mean at x.  It should get smaller as you increase the number of
  19 ##   measurements.  Error bars in the physical sciences usually show
  20 ##   a 1-alpha confidence value of erfc(1/sqrt(2)), representing
  21 ##   one standandard deviation of uncertainty in the mean.
  22 ##
  23 ## polyconf(...,1-alpha)
  24 ##
  25 ##   Control the width of the interval. If asking for the prediction
  26 ##   interval 'pi', the default is .05 for the 95% prediction interval.
  27 ##   If asking for the confidence interval 'ci', the default is
  28 ##   erfc(1/sqrt(2)) for a one standard deviation confidence interval.
  29 ##
  30 ## Example:
  31 ##  [p,s] = polyfit(x,y,1);
  32 ##  xf = linspace(x(1),x(end),150);
  33 ##  [yf,dyf] = polyconf(p,xf,s,'ci');
  34 ##  plot(xf,yf,'g-;fit;',xf,yf+dyf,'g.;;',xf,yf-dyf,'g.;;',x,y,'xr;data;');
  35 ##  plot(x,y-polyval(p,x),';residuals;',xf,dyf,'g-;;',xf,-dyf,'g-;;');
  36
  37 function [y,dy] = polyconf(p,x,varargin)
  38   alpha = s = [];
  39   typestr = 'pi';
  40   for i=1:length(varargin)
  41     v = varargin{i};
  42     if isstruct(v), s = v;
  43     elseif ischar(v), typestr = v;
  44     elseif isscalar(v), alpha = v;
  45     else s = [];
  46     end
  47   end
  48   if (nargout>1 && (isempty(s)||nargin<3)) || nargin < 2
  49     print_usage;
  50   end
  51
  52   if isempty(s)
  53     y = polyval(p,x);
  54
  55   else
  56     ## For a polynomial fit, x is the set of powers ( x^n ; ... ; 1 ).
  57     n=length(p)-1;
  58     k=length(x(:));
  59     if columns(s.R) == n, ## fit through origin
  60       A = (x(:) * ones (1, n)) .^ (ones (k, 1) * (n:-1:1));
  61       p = p(1:n);
  62     else
  63       A = (x(:) * ones (1, n+1)) .^ (ones (k, 1) * (n:-1:0));
  64     endif
  65     y = dy = x;
  66     [y(:),dy(:)] = confidence(A,p,s,alpha,typestr);
  67
  68   end
  69 end
  70
  71 %!test
  72 %! # data from Hocking, RR, "Methods and Applications of Linear Models"
  73 %! temperature=[40;40;40;45;45;45;50;50;50;55;55;55;60;60;60;65;65;65];
  74 %! 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];
  75 %! [p,s] = polyfit(temperature,strength,1);
  76 %! [y,dy] = polyconf(p,40,s,0.05,'ci');
  77 %! assert([y,dy],[66.15396825396826,1.71702862681486],200*eps);
  78 %! [y,dy] = polyconf(p,40,s,0.05,'pi');
  79 %! assert(dy,4.45345484470743,200*eps);
  80
  81 ## [y,dy] = confidence(A,p,s)
  82 ##
  83 ##   Produce prediction intervals for the fitted y. The vector p
  84 ##   and structure s are returned from wsolve. The matrix A is
  85 ##   the set of observation values at which to evaluate the
  86 ##   confidence interval.
  87 ##
  88 ## confidence(...,['ci'|'pi'])
  89 ##
  90 ##   Produce a confidence interval (range of likely values for the
  91 ##   mean at x) or a prediction interval (range of likely values
  92 ##   seen when measuring at x).  The prediction interval tells
  93 ##   you the width of the distribution at x.  This should be the same
  94 ##   regardless of the number of measurements you have for the value
  95 ##   at x.  The confidence interval tells you how well you know the
  96 ##   mean at x.  It should get smaller as you increase the number of
  97 ##   measurements.  Error bars in the physical sciences usually show
  98 ##   a 1-alpha confidence value of erfc(1/sqrt(2)), representing
  99 ##   one standandard deviation of uncertainty in the mean.
 100 ##
 101 ## confidence(...,1-alpha)
 102 ##
 103 ##   Control the width of the interval. If asking for the prediction
 104 ##   interval 'pi', the default is .05 for the 95% prediction interval.
 105 ##   If asking for the confidence interval 'ci', the default is
 106 ##   erfc(1/sqrt(2)) for a one standard deviation confidence interval.
 107 ##
 108 ## Confidence intervals for linear system are given by:
 109 ##    x' p +/- sqrt( Finv(1-a,1,df) var(x' p) )
 110 ## where for confidence intervals,
 111 ##    var(x' p) = sigma^2 (x' inv(A'A) x)
 112 ## and for prediction intervals,
 113 ##    var(x' p) = sigma^2 (1 + x' inv(A'A) x)
 114 ##
 115 ## Rather than A'A we have R from the QR decomposition of A, but
 116 ## R'R equals A'A.  Note that R is not upper triangular since we
 117 ## have already multiplied it by the permutation matrix, but it
 118 ## is invertible.  Rather than forming the product R'R which is
 119 ## ill-conditioned, we can rewrite x' inv(A'A) x as the equivalent
 120 ##    x' inv(R) inv(R') x = t t', for t = x' inv(R)
 121 ## Since x is a vector, t t' is the inner product sumsq(t).
 122 ## Note that LAPACK allows us to do this simultaneously for many
 123 ## different x using sqrt(sumsq(X/R,2)), with each x on a different row.
 124 ##
 125 ## Note: sqrt(F(1-a;1,df)) = T(1-a/2;df)
 126 ##
 127 ## For non-linear systems, use x = dy/dp and ignore the y output.
 128 function [y,dy] = confidence(A,p,S,alpha,typestr)
 129   if nargin < 4, alpha = []; end
 130   if nargin < 5, typestr = 'ci'; end
 131   y = A*p(:);
 132   switch typestr,
 133     case 'ci', pred = 0; default_alpha=erfc(1/sqrt(2));
 134     case 'pi', pred = 1; default_alpha=0.05;
 135     otherwise, error("use 'ci' or 'pi' for interval type");
 136   end
 137   if isempty(alpha), alpha = default_alpha; end
 138   s = tinv(1-alpha/2,S.df)*S.normr/sqrt(S.df);
 139   dy = s*sqrt(pred+sumsq(A/S.R,2));
 140 end