octave_packages/data-smoothing-1.3.0/rgdtsmcore.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{v}] =} rgdtsmcore (@var{x}, @var{y}, @var{d}, @var{lambda}, [@var{options}])
  18 ##
  19 ## Smooths @var{y} vs. @var{x} values by Tikhonov regularization.
  20 ## Although this function can be used directly, the more feature rich
  21 ## function "regdatasmooth" should be used instead.  In addition to
  22 ## @var{x} and @var{y}, required input includes the smoothing derivative
  23 ## @var{d} and the regularization parameter @var{lambda}.  The smooth
  24 ## y-values are returned as @var{yhat}.  The generalized cross
  25 ## validation variance @var{v} may also be returned.
  26 ##
  27 ## Note:  the options have changed!
  28 ## Currently supported input options are (multiple options are allowed):
  29 ##
  30 ## @table @code
  31 ## @item "xhat", @var{vector}
  32 ##  A vector of x-values to use for the smooth curve; must be
  33 ##  monotonically increasing and must at least span the data
  34 ## @item "weights", @var{vector}
  35 ##  A vector of weighting values for fitting each point in the data.
  36 ## @item "relative"
  37 ##  use relative differences for the goodnes of fit term.  Conflicts
  38 ##  with the "weights" option.
  39 ## @item "midpointrule"
  40 ##  use the midpoint rule for the integration terms rather than a direct
  41 ##  sum; this option conflicts with the option "xhat"
  42 ## @end table
  43 ##
  44 ## References:  Anal. Chem. (2003) 75, 3631; AIChE J. (2006) 52, 325
  45 ## @seealso{regdatasmooth}
  46 ## @end deftypefn
  47
  48
  49 function [yhat, v] = rgdtsmcore (x, y, d, lambda, varargin)
  50
  51   if (nargin < 4)
  52     print_usage;
  53   endif
  54
  55   ## Defaults if not provided
  56   xhatprov = 0;
  57   xhat = x;
  58   weights = 0;
  59   relative = 0;
  60   midpr = 0;
  61
  62   ## parse the provided options
  63   if ( length(varargin) )
  64     for i = 1:length(varargin)
  65       arg = varargin{i};
  66       if ischar(arg)
  67         switch arg
  68           case "xhat"
  69             xhatprov = 1;
  70             xhat = varargin{i+1};
  71           case "weights"
  72             weights = 1;
  73             weightv = varargin{i+1};
  74           case "relative"
  75             relative = 1;
  76           case "midpointrule"
  77             midpr = 1;
  78           otherwise
  79             printf("Option '%s' is not implemented;\n", arg)
  80         endswitch
  81       endif
  82     endfor
  83   endif
  84   if (xhatprov && midpr)
  85     warning("midpointrule is currently not used if xhat is provided (since x,y may be scattered)")
  86     midpr = 0;
  87   endif
  88   if (weights && relative)
  89     warning("relative differences is not used if a weighting vector is provided")
  90   endif
  91
  92   N = length(x);
  93   Nhat = length(xhat);
  94
  95   ## test that xhat is increasing
  96   if !all(diff(xhat)>0)
  97     if xhatprov
  98       error("xhat must be monotonically increasing")
  99     else
 100       error("x must be monotonically increasing if xhat is not provided")
 101     endif
 102   endif
 103   ## test that xhat spans x
 104   if ( min(x) < min(xhat) || max(xhat) < max(x) )
 105     error("xhat must at least span the data")
 106   endif
 107
 108   ## construct M, D
 109   M = speye(Nhat);
 110   idx = interp1(xhat,1:Nhat,x,"nearest"); # works for unequally spaced xhat
 111   M = M(idx,:);
 112   D = ddmat(xhat,d);
 113
 114   ## construct "weighting" matrices W and U
 115   if (weights)
 116     ## use arbitrary weighting as provided
 117     W = diag(weightv);
 118   elseif (relative)
 119     ## use relative differences
 120     Yinv = sparse(diag(1./y));
 121     W = Yinv^2;
 122   else
 123     W = speye(N);
 124   endif
 125   ## use midpoint rule integration (rather than simple sums)
 126   if (midpr)
 127     Bhat = sparse(diag(-ones(N-1,1),-1)) + sparse(diag(ones(N-1,1),1));
 128     Bhat(1,1) = -1;
 129     Bhat(N,N) = 1;
 130     B = 1/2*sparse(diag(Bhat*x));
 131     if ( floor(d/2) == d/2 ) # test if d is even
 132       dh = d/2;
 133       Btilda = B(dh+1:N-dh,dh+1:N-dh);
 134     else # d is odd
 135       dh = ceil(d/2);
 136       Btilda = B(dh:N-dh,dh:N-dh);
 137     endif
 138     W = W*B;
 139     U = Btilda;
 140   else
 141     ## W = W*speye(Nhat);
 142     U = speye(Nhat-d);
 143   endif
 144
 145   ## Smoothing
 146   delta = trace(D'*D)/Nhat^(2+d);  # using "relative" or other weighting affects this!
 147   yhat = (M'*W*M + lambda*delta^(-1)*D'*U*D) \ M'*W*y;
 148   #[R,P,S] = splchol(M'*W*M + lambda*delta^(-1)*D'*U*D);
 149   #yhat = S*(R'\(R\(S'*M'*W*y)));
 150
 151   ## Computation of hat diagonal and cross-validation
 152   if (nargout > 1)
 153     ## from AIChE J. (2006) 52, 325
 154     ## note: chol factorization does not help speed up the computation of H;
 155     ## should implement Eiler's partial H computation if many point smoothing by GCV is needed
 156     ##H = M*(S*(R'\(R\(S'*M'*W))));
 157     H = M*((M'*W*M + lambda*delta^(-1)*D'*U*D)\M'*W);
 158     ## note: this is variance, squared of the standard error that Eilers uses
 159     v = (M*yhat - y)'*(M*yhat - y)/N / (1 - trace(H)/N)^2;
 160   endif
 161
 162   ## test mapping
 163   ##figure(5)
 164   ##plot(x,y,'o',x,M*yhat,'x')
 165
 166 endfunction