This file is indexed.

/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