X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?p=CreaPhase.git;a=blobdiff_plain;f=octave_packages%2Foptim-1.2.0%2Fwpolyfit.m;fp=octave_packages%2Foptim-1.2.0%2Fwpolyfit.m;h=a9563fa27a39feeed9b1d6be8c69f08aade32336;hp=0000000000000000000000000000000000000000;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hpb=1705066eceaaea976f010f669ce8e972f3734b05

diff --git a/octave_packages/optim-1.2.0/wpolyfit.m b/octave_packages/optim-1.2.0/wpolyfit.m
new file mode 100644
index 0000000..a9563fa
--- /dev/null
+++ b/octave_packages/optim-1.2.0/wpolyfit.m
@@ -0,0 +1,261 @@
+## Author: Paul Kienzle <pkienzle@gmail.com>
+## This program is granted to the public domain.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{p}, @var{s}] =} wpolyfit (@var{x}, @var{y}, @var{dy}, @var{n})
+## Return the coefficients of a polynomial @var{p}(@var{x}) of degree
+## @var{n} that minimizes
+## @iftex
+## @tex
+## $$
+## \sum_{i=1}^N (p(x_i) - y_i)^2
+## $$
+## @end tex
+## @end iftex
+## @ifinfo
+## @code{sumsq (p(x(i)) - y(i))},
+## @end ifinfo
+## to best fit the data in the least squares sense.  The standard error
+## on the observations @var{y} if present are given in @var{dy}.
+##
+## The returned value @var{p} contains the polynomial coefficients 
+## suitable for use in the function polyval.  The structure @var{s} returns
+## information necessary to compute uncertainty in the model.
+##
+## To compute the predicted values of y with uncertainty use
+## @example
+## [y,dy] = polyconf(p,x,s,'ci');
+## @end example
+## You can see the effects of different confidence intervals and
+## prediction intervals by calling the wpolyfit internal plot
+## function with your fit:
+## @example
+## feval('wpolyfit:plt',x,y,dy,p,s,0.05,'pi')
+## @end example
+## Use @var{dy}=[] if uncertainty is unknown.
+##
+## You can use a chi^2 test to reject the polynomial fit:
+## @example
+## p = 1-chi2cdf(s.normr^2,s.df);
+## @end example
+## p is the probability of seeing a chi^2 value higher than that which 
+## was observed assuming the data are normally distributed around the fit.
+## If p < 0.01, you can reject the fit at the 1% level.
+##
+## You can use an F test to determine if a higher order polynomial 
+## improves the fit:
+## @example
+## [poly1,S1] = wpolyfit(x,y,dy,n);
+## [poly2,S2] = wpolyfit(x,y,dy,n+1);
+## F = (S1.normr^2 - S2.normr^2)/(S1.df-S2.df)/(S2.normr^2/S2.df);
+## p = 1-f_cdf(F,S1.df-S2.df,S2.df);
+## @end example
+## p is the probability of observing the improvement in chi^2 obtained
+## by adding the extra parameter to the fit.  If p < 0.01, you can reject 
+## the lower order polynomial at the 1% level.
+##
+## You can estimate the uncertainty in the polynomial coefficients 
+## themselves using
+## @example
+## dp = sqrt(sumsq(inv(s.R'))'/s.df)*s.normr;
+## @end example
+## but the high degree of covariance amongst them makes this a questionable
+## operation.
+##
+## @deftypefnx {Function File} {[@var{p}, @var{s}, @var{mu}] =} wpolyfit (...)
+##
+## If an additional output @code{mu = [mean(x),std(x)]} is requested then 
+## the @var{x} values are centered and normalized prior to computing the fit.
+## This will give more stable numerical results.  To compute a predicted 
+## @var{y} from the returned model use
+## @code{y = polyval(p, (x-mu(1))/mu(2)}
+##
+## @deftypefnx {Function File} wpolyfit (...)
+##
+## If no output arguments are requested, then wpolyfit plots the data,
+## the fitted line and polynomials defining the standard error range.
+##
+## Example
+## @example
+## x = linspace(0,4,20);
+## dy = (1+rand(size(x)))/2;
+## y = polyval([2,3,1],x) + dy.*randn(size(x));
+## wpolyfit(x,y,dy,2);
+## @end example
+##
+## @deftypefnx {Function File} wpolyfit (..., 'origin')
+##
+## If 'origin' is specified, then the fitted polynomial will go through
+## the origin.  This is generally ill-advised.  Use with caution.
+##
+## Hocking, RR (2003). Methods and Applications of Linear Models.
+## New Jersey: John Wiley and Sons, Inc.
+##
+## @end deftypefn
+## @seealso{polyfit,polyconf}
+
+function [p_out, s, mu] = wpolyfit (varargin)
+
+  ## strip 'origin' of the end
+  args = length(varargin);
+  if args>0 && ischar(varargin{args})
+    origin = varargin{args};
+    args--;
+  else
+    origin='';
+  endif
+  ## strip polynomial order off the end
+  if args>0
+    n = varargin{args};
+    args--;
+  else
+    n = [];
+  end
+  ## interpret the remainder as x,y or x,y,dy or [x,y] or [x,y,dy]
+  if args == 3
+    x = varargin{1};
+    y = varargin{2};
+    dy = varargin{3};
+  elseif args == 2
+    x = varargin{1};
+    y = varargin{2};
+    dy = [];
+  elseif args == 1
+    A = varargin{1};
+    [nr,nc]=size(A);
+    if all(nc!=[2,3])
+      error("wpolyfit expects vectors x,y,dy or matrix [x,y,dy]");
+    endif
+    dy = [];
+    if nc == 3, dy = A(:,3); endif
+    y = A(:,2);
+    x = A(:,1);
+  else
+    usage ("wpolyfit (x, y [, dy], n [, 'origin'])");
+  end
+
+  if (length(origin) == 0)
+    through_origin = 0;
+  elseif strcmp(origin,'origin')
+    through_origin = 1;
+  else
+    error ("wpolyfit: expected 'origin' but found '%s'", origin)
+  endif
+
+  if any(size (x) != size (y))
+    error ("wpolyfit: x and y must be vectors of the same size");
+  endif
+  if length(dy)>1 && length(y) != length(dy)
+    error ("wpolyfit: dy must be a vector the same length as y");
+  endif
+
+  if (! (isscalar (n) && n >= 0 && ! isinf (n) && n == round (n)))
+    error ("wpolyfit: n must be a nonnegative integer");
+  endif
+
+  if nargout == 3
+    mu = [mean(x), std(x)];
+    x = (x - mu(1))/mu(2);
+  endif
+
+  k = length (x);
+
+  ## observation matrix
+  if through_origin
+    ## polynomial through the origin y = ax + bx^2 + cx^3 + ...
+    A = (x(:) * ones(1,n)) .^ (ones(k,1) * (n:-1:1));
+  else
+    ## polynomial least squares y = a + bx + cx^2 + dx^3 + ...
+    A = (x(:) * ones (1, n+1)) .^ (ones (k, 1) * (n:-1:0));
+  endif
+
+  [p,s] = wsolve(A,y(:),dy(:));
+
+  if through_origin
+    p(n+1) = 0;
+  endif
+
+  if nargout == 0
+    good_fit = 1-chi2cdf(s.normr^2,s.df);
+    printf("Polynomial: %s  [ p(chi^2>observed)=%.2f%% ]\n", polyout(p,'x'), good_fit*100);
+    plt(x,y,dy,p,s,'ci');
+  else
+    p_out = p';
+  endif
+
+function plt(x,y,dy,p,s,varargin)
+
+  if iscomplex(p)
+    # XXX FIXME XXX how to plot complex valued functions?
+    # Maybe using hue for phase and saturation for magnitude
+    # e.g., Frank Farris (Santa Cruz University) has this:
+    # http://www.maa.org/pubs/amm_complements/complex.html
+    # Could also look at the book
+    #   Visual Complex Analysis by Tristan Needham, Oxford Univ. Press
+    # but for now we punt
+    return
+  end
+
+  ## decorate the graph
+  grid('on');
+  xlabel('abscissa X'); ylabel('data Y');
+  title('Least-squares Polynomial Fit with Error Bounds');
+
+  ## draw fit with estimated error bounds
+  xf = linspace(min(x),max(x),150)';
+  [yf,dyf] = polyconf(p,xf,s,varargin{:});
+  plot(xf,yf+dyf,"g.;;", xf,yf-dyf,"g.;;", xf,yf,"g-;fit;");
+
+  ## plot the data
+  hold on;
+  if (isempty(dy))
+    plot(x,y,"x;data;");
+  else
+    if isscalar(dy), dy = ones(size(y))*dy; end
+    errorbar (x, y, dy, "~;data;");
+  endif
+  hold off;
+
+  if strcmp(deblank(input('See residuals? [y,n] ','s')),'y')
+    clf;
+    if (isempty(dy))
+      plot(x,y-polyval(p,x),"x;data;");
+    else
+      errorbar(x,y-polyval(p,x),dy, '~;data;');
+    endif
+    hold on;
+    grid on;
+    ylabel('Residuals');
+    xlabel('abscissa X'); 
+    plot(xf,dyf,'g.;;',xf,-dyf,'g.;;');
+    hold off;
+  endif
+
+%!demo % #1  
+%! x = linspace(0,4,20);
+%! dy = (1+rand(size(x)))/2;
+%! y = polyval([2,3,1],x) + dy.*randn(size(x));
+%! wpolyfit(x,y,dy,2);
+  
+%!demo % #2
+%! x = linspace(-i,+2i,20);
+%! noise = ( randn(size(x)) + i*randn(size(x)) )/10;
+%! P = [2-i,3,1+i];
+%! y = polyval(P,x) + noise;
+%! wpolyfit(x,y,2)
+
+%!demo
+%! pin = [3; -1; 2];
+%! x = -3:0.1:3;
+%! y = polyval (pin, x);
+%!
+%! ## Poisson weights
+%! # dy = sqrt (abs (y));
+%! ## Uniform weights in [0.5,1]
+%! dy = 0.5 + 0.5 * rand (size (y));
+%!
+%! y = y + randn (size (y)) .* dy;
+%! printf ("Original polynomial: %s\n", polyout (pin, 'x'));
+%! wpolyfit (x, y, dy, length (pin)-1);
+
+