octave_packages/data-smoothing-1.3.0/regdatasmooth.m

   1 ## Copyright (C) 2008 Jonathan Stickel <jonathan.stickel@nrel.gov>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ##@deftypefn {Function File} {[@var{yhat}, @var{lambda}] =} regdatasmooth (@var{x}, @var{y}, [@var{options}])
  18 ##
  19 ## Smooths the @var{y} vs. @var{x} values of 1D data by Tikhonov
  20 ## regularization. The smooth y-values are returned as @var{yhat}. The
  21 ## regularization parameter @var{lambda} that was used for the smoothing
  22 ## may also be returned.
  23 ##
  24 ## Note:  the options have changed!
  25 ## Currently supported input options are (multiple options are allowed):
  26 ##
  27 ##@table @code
  28 ##@item "d", @var{value}
  29 ## the smoothing derivative to use (default = 2)
  30 ##@item "lambda", @var{value}
  31 ## the regularization paramater to use
  32 ##@item "stdev", @var{value}
  33 ## the standard deviation of the measurement of @var{y}; an optimal
  34 ## value for lambda will be determined by matching the provided
  35 ## @var{value} with the standard devation of @var{yhat}-@var{y};
  36 ## if the option "relative" is also used, then a relative standard
  37 ## deviation is inferred
  38 ##@item "gcv"
  39 ## use generalized cross-validation to determine the optimal value for
  40 ## lambda; if neither "lambda" nor "stdev" options are given, this
  41 ## option is implied
  42 ##@item "lguess", @var{value}
  43 ## the initial value for lambda to use in the iterative minimization
  44 ## algorithm to find the optimal value (default = 1)
  45 ## @item "xhat", @var{vector}
  46 ##  A vector of x-values to use for the smooth curve; must be
  47 ##  monotonically increasing and must at least span the data
  48 ## @item "weights", @var{vector}
  49 ##  A vector of weighting values for fitting each point in the data.
  50 ## @item "relative"
  51 ##  use relative differences for the goodnes of fit term.  Conflicts
  52 ##  with the "weights" option.
  53 ##@item "midpointrule"
  54 ## use the midpoint rule for the integration terms rather than a direct
  55 ## sum; this option conflicts with the option "xhat"
  56 ##@end table
  57 ##
  58 ## Please run the demos for example usage.
  59 ##
  60 ## References:  Anal. Chem. (2003) 75, 3631; AIChE J. (2006) 52, 325
  61 ## @seealso{rgdtsmcorewrap, rgdtsmcore}
  62 ## @end deftypefn
  63
  64 function [yhat, lambda] = regdatasmooth (x, y, varargin)
  65
  66   if (nargin < 2)
  67     print_usage;
  68   elseif ( length(x) != length(y) )
  69     error("x and y must be equal length vectors")
  70   endif
  71   if ( isrow(x) ) x = x'; endif
  72   if ( isrow(y) ) y = y'; endif
  73
  74   ## defaults
  75   d = 2;
  76   lambda = 0;
  77   stdev = 0;
  78   guess = 0;
  79
  80   ## parse options for d, lambda, stdev, gcv, lguess;
  81   ## remaining options (gridx, Nhat, range, relative, midpointrule)
  82   ## will be sent directly to the core function
  83   idx = [];
  84   if ( nargin > 2)
  85     for i = 1:nargin-2
  86       arg = varargin{i};
  87       if ischar(arg)
  88         switch arg
  89           case "d"
  90             d = varargin{i+1};
  91             idx = [idx,i,i+1];
  92           case "lambda"
  93             lambda = varargin{i+1};
  94             idx = [idx,i,i+1];
  95           case "stdev"
  96             stdev = varargin{i+1};
  97             idx = [idx,i,i+1];
  98           case "gcv"
  99             idx = [idx,i];
 100           case "lguess"
 101             guess = log10(varargin{i+1});
 102             idx = [idx,i,i+1];
 103         endswitch
 104       endif
 105     endfor
 106   endif
 107   varargin(idx) = [];
 108   options = varargin;
 109   ## add warning if more than one gcv, lambda, or stdev options provided?
 110
 111   maxiter = 50;
 112   if (lambda)
 113     ## do nothing and use the provided lambda
 114   else
 115     ## find the "optimal" lambda
 116     if ( stdev )
 117       ## match standard deviation
 118       fhandle = @(log10lambda) rgdtsmcorewrap (log10lambda, x, y, d, {"stdev", stdev}, options{:});
 119     else
 120       ## perform cross-validation
 121       fhandle = @(log10lambda) rgdtsmcorewrap (log10lambda, x, y, d, {"cve"}, options{:});
 122     endif
 123     ## "fminunc" works OK, but a derivative-free method (below) is better for this problem
 124     ##opt = optimset("TolFun",1e-6,"MaxFunEvals",maxiter);
 125     ##[log10lambda,fout,exitflag] = fminunc (fhandle, guess, opt);
 126     ##[log10lambda,fout,exitflag] = fminunc_compat (fhandle, guess, opt);
 127     ## derivative-free optimization; should use "fminsearch" for Matlab
 128     ## compatibility, but fminsearch needs updates to be more compatible itself
 129     [log10lambda, fout, niter] = nelder_mead_min (fhandle, guess, "ftol", 1e-6, "maxev", maxiter);
 130     if (niter > maxiter)
 131       exitflag = 0;
 132     else
 133       exitflag = 1;
 134     endif
 135     if (!exitflag)
 136       warning("Iteration limit of %i exceeded\n",maxiter)
 137     endif
 138     lambda = 10^log10lambda;
 139   endif
 140
 141   yhat = rgdtsmcore (x, y, d, lambda, options{:});
 142
 143 endfunction
 144
 145 %!demo
 146 %! npts = 100;
 147 %! x = linspace(0,2*pi,npts)';
 148 %! x = x + 2*pi/npts*(rand(npts,1)-0.5);
 149 %! y = sin(x);
 150 %! y = y + 1e-1*randn(npts,1);
 151 %! yp = ddmat(x,1)*y;
 152 %! y2p = ddmat(x,2)*y;
 153 %! [yh, lambda] = regdatasmooth (x, y, "d",4,"stdev",1e-1,"midpointrule");
 154 %! lambda
 155 %! yhp = ddmat(x,1)*yh;
 156 %! yh2p = ddmat(x,2)*yh;
 157 %! clf
 158 %! subplot(221)
 159 %! plot(x,y,'o','markersize',5,x,yh,x,sin(x))
 160 %! title("y(x)")
 161 %! legend("noisy","smoothed","sin(x)","location","northeast");
 162 %! subplot(222)
 163 %! plot(x(1:end-1),[yp,yhp,cos(x(1:end-1))])
 164 %! axis([min(x),max(x),min(yhp)-abs(min(yhp)),max(yhp)*2])
 165 %! title("y'(x)")
 166 %! legend("noisy","smoothed","cos(x)","location","southeast");
 167 %! subplot(223)
 168 %! plot(x(2:end-1),[y2p,yh2p,-sin(x(2:end-1))])
 169 %! axis([min(x),max(x),min(yh2p)-abs(min(yh2p)),max(yh2p)*2])
 170 %! title("y''(x)")
 171 %! legend("noisy","smoothed","-sin(x)","location","southeast");
 172 %! %--------------------------------------------------------
 173 %! % smoothing of monotonic data, using "stdev" to determine the optimal lambda
 174
 175 %!demo
 176 %! npts = 20;
 177 %! x = rand(npts,1)*2*pi;
 178 %! y = sin(x);
 179 %! y = y + 1e-1*randn(npts,1);
 180 %! xh = linspace(0,2*pi,200)';
 181 %! [yh, lambda] = regdatasmooth (x, y, "d", 3, "xhat", xh);
 182 %! lambda
 183 %! clf
 184 %! figure(1);
 185 %! plot(x,y,'o','markersize',10,xh,yh,xh,sin(xh))
 186 %! title("y(x)")
 187 %! legend("noisy","smoothed","sin(x)","location","northeast");
 188 %! %--------------------------------------------------------
 189 %! % smoothing of scattered data, using "gcv" to determine the optimal lambda