octave-data-smoothing 1.3.0-4 / usr / share / octave / packages / data-smoothing-1.3.0 / regdatasmooth.m

## Copyright (C) 2008 Jonathan Stickel <jonathan.stickel@nrel.gov>
##
## 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 3 of the License, 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; if not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
##@deftypefn {Function File} {[@var{yhat}, @var{lambda}] =} regdatasmooth (@var{x}, @var{y}, [@var{options}])
##
## Smooths the @var{y} vs. @var{x} values of 1D data by Tikhonov
## regularization. The smooth y-values are returned as @var{yhat}. The
## regularization parameter @var{lambda} that was used for the smoothing
## may also be returned.
##
## Note:  the options have changed!
## Currently supported input options are (multiple options are allowed):
##
##@table @code
##@item "d", @var{value}
## the smoothing derivative to use (default = 2)
##@item "lambda", @var{value}
## the regularization paramater to use
##@item "stdev", @var{value}
## the standard deviation of the measurement of @var{y}; an optimal
## value for lambda will be determined by matching the provided
## @var{value} with the standard devation of @var{yhat}-@var{y};
## if the option "relative" is also used, then a relative standard
## deviation is inferred
##@item "gcv"
## use generalized cross-validation to determine the optimal value for
## lambda; if neither "lambda" nor "stdev" options are given, this
## option is implied
##@item "lguess", @var{value}
## the initial value for lambda to use in the iterative minimization
## algorithm to find the optimal value (default = 1)
## @item "xhat", @var{vector}
##  A vector of x-values to use for the smooth curve; must be
##  monotonically increasing and must at least span the data
## @item "weights", @var{vector}
##  A vector of weighting values for fitting each point in the data.
## @item "relative"
##  use relative differences for the goodnes of fit term.  Conflicts
##  with the "weights" option.
##@item "midpointrule"
## use the midpoint rule for the integration terms rather than a direct
## sum; this option conflicts with the option "xhat"
##@end table
##
## Please run the demos for example usage.
##
## References:  Anal. Chem. (2003) 75, 3631; AIChE J. (2006) 52, 325
## @seealso{rgdtsmcorewrap, rgdtsmcore}
## @end deftypefn

function [yhat, lambda] = regdatasmooth (x, y, varargin)

  if (nargin < 2)
    print_usage;
  elseif ( length(x) != length(y) )
    error("x and y must be equal length vectors")
  endif
  if ( isrow(x) ) x = x'; endif
  if ( isrow(y) ) y = y'; endif

  ## defaults
  d = 2;
  lambda = 0;
  stdev = 0;
  guess = 0;

  ## parse options for d, lambda, stdev, gcv, lguess;
  ## remaining options (gridx, Nhat, range, relative, midpointrule)
  ## will be sent directly to the core function
  idx = [];
  if ( nargin > 2)
    for i = 1:nargin-2
      arg = varargin{i};
      if ischar(arg)
        switch arg
          case "d"
            d = varargin{i+1};
            idx = [idx,i,i+1];
          case "lambda"
            lambda = varargin{i+1};
            idx = [idx,i,i+1];
          case "stdev"
            stdev = varargin{i+1};
            idx = [idx,i,i+1];
          case "gcv"
            idx = [idx,i];
          case "lguess"
            guess = log10(varargin{i+1});
            idx = [idx,i,i+1];
        endswitch
      endif
    endfor
  endif
  varargin(idx) = [];
  options = varargin;
  ## add warning if more than one gcv, lambda, or stdev options provided?

  maxiter = 50;
  if (lambda)
    ## do nothing and use the provided lambda
  else
    ## find the "optimal" lambda
    if ( stdev )
      ## match standard deviation
      fhandle = @(log10lambda) rgdtsmcorewrap (log10lambda, x, y, d, {"stdev", stdev}, options{:});
    else
      ## perform cross-validation
      fhandle = @(log10lambda) rgdtsmcorewrap (log10lambda, x, y, d, {"cve"}, options{:});
    endif
    ## "fminunc" works OK, but a derivative-free method (below) is better for this problem
    ##opt = optimset("TolFun",1e-6,"MaxFunEvals",maxiter);
    ##[log10lambda,fout,exitflag] = fminunc (fhandle, guess, opt);
    ##[log10lambda,fout,exitflag] = fminunc_compat (fhandle, guess, opt);
    ## derivative-free optimization; should use "fminsearch" for Matlab
    ## compatibility, but fminsearch needs updates to be more compatible itself
    [log10lambda, fout, niter] = nelder_mead_min (fhandle, guess, "ftol", 1e-6, "maxev", maxiter);
    if (niter > maxiter)
      exitflag = 0;
    else
      exitflag = 1;
    endif
    if (!exitflag)
      warning("Iteration limit of %i exceeded\n",maxiter)
    endif
    lambda = 10^log10lambda;
  endif
  
  yhat = rgdtsmcore (x, y, d, lambda, options{:});
  
endfunction

%!demo
%! npts = 100;
%! x = linspace(0,2*pi,npts)';
%! x = x + 2*pi/npts*(rand(npts,1)-0.5);
%! y = sin(x);
%! y = y + 1e-1*randn(npts,1);
%! yp = ddmat(x,1)*y;
%! y2p = ddmat(x,2)*y;
%! [yh, lambda] = regdatasmooth (x, y, "d",4,"stdev",1e-1,"midpointrule");
%! lambda
%! yhp = ddmat(x,1)*yh;  
%! yh2p = ddmat(x,2)*yh;
%! clf
%! subplot(221)
%! plot(x,y,'o','markersize',5,x,yh,x,sin(x))
%! title("y(x)")
%! legend("noisy","smoothed","sin(x)","location","northeast");
%! subplot(222)
%! plot(x(1:end-1),[yp,yhp,cos(x(1:end-1))])
%! axis([min(x),max(x),min(yhp)-abs(min(yhp)),max(yhp)*2])
%! title("y'(x)")
%! legend("noisy","smoothed","cos(x)","location","southeast");
%! subplot(223)
%! plot(x(2:end-1),[y2p,yh2p,-sin(x(2:end-1))])
%! axis([min(x),max(x),min(yh2p)-abs(min(yh2p)),max(yh2p)*2])
%! title("y''(x)")
%! legend("noisy","smoothed","-sin(x)","location","southeast");
%! %--------------------------------------------------------
%! % smoothing of monotonic data, using "stdev" to determine the optimal lambda

%!demo
%! npts = 20;
%! x = rand(npts,1)*2*pi;
%! y = sin(x);
%! y = y + 1e-1*randn(npts,1);
%! xh = linspace(0,2*pi,200)';
%! [yh, lambda] = regdatasmooth (x, y, "d", 3, "xhat", xh);
%! lambda
%! clf
%! figure(1);
%! plot(x,y,'o','markersize',10,xh,yh,xh,sin(xh))
%! title("y(x)")
%! legend("noisy","smoothed","sin(x)","location","northeast");
%! %--------------------------------------------------------
%! % smoothing of scattered data, using "gcv" to determine the optimal lambda