## Copyright (C) 2008 Jonathan Stickel ## 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 2 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{v}] =} rgdtsmcore (@var{x}, @var{y}, @var{d}, @var{lambda}, [@var{options}]) ## ## Smooths @var{y} vs. @var{x} values by Tikhonov regularization. ## Although this function can be used directly, the more feature rich ## function "regdatasmooth" should be used instead. In addition to ## @var{x} and @var{y}, required input includes the smoothing derivative ## @var{d} and the regularization parameter @var{lambda}. The smooth ## y-values are returned as @var{yhat}. The generalized cross ## validation variance @var{v} may also be returned. ## ## Note: the options have changed! ## Currently supported input options are (multiple options are allowed): ## ##@table @code ##@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 ## ## References: Anal. Chem. (2003) 75, 3631; AIChE J. (2006) 52, 325 ##@end deftypefn function [yhat, v] = rgdtsmcore (x, y, d, lambda, varargin) ## Defaults if not provided xhatprov = 0; xhat = x; weights = 0; relative = 0; midpr = 0; ## parse the provided options if ( length(varargin) ) for i = 1:length(varargin) arg = varargin{i}; if ischar(arg) switch arg case "xhat" xhatprov = 1; xhat = varargin{i+1}; case "weights" weights = 1; weightv = varargin{i+1}; case "relative" relative = 1; case "midpointrule" midpr = 1; otherwise printf("Option '%s' is not implemented;\n", arg) endswitch endif endfor endif if (xhatprov & midpr) warning("midpointrule is currently not used if xhat is provided (since x,y may be scattered)") midpr = 0; endif if (weights & relative) warning("relative differences is not used if a weighting vector is provided") endif N = length(x); Nhat = length(xhat); ## test that xhat is increasing if !all(diff(xhat)>0) if xhatprov error("xhat must be monotonically increasing") else error("x must be monotonically increasing if xhat is not provided") endif endif ## test that xhat spans x if ( min(x) < min(xhat) | max(xhat) < max(x) ) error("xhat must at least span the data") endif ## construct M, D M = speye(Nhat); idx = interp1(xhat,1:Nhat,x,"nearest"); # works for unequal spaced xhat M = M(idx,:); D = ddmat(xhat,d); ## construct "weighting" matrices W and U if (weights) ## use arbitrary weighting as provided W = diag(weightv); elseif (relative) ## use relative differences Yinv = spdiag(1./y); W = Yinv^2; else W = speye(N); endif ## use midpoint rule integration (rather than simple sums) if (midpr) Bhat = spdiag(-ones(N-1,1),-1) + spdiag(ones(N-1,1),1); Bhat(1,1) = -1; Bhat(N,N) = 1; B = 1/2*spdiag(Bhat*x); if ( floor(d/2) == d/2 ) # test if d is even dh = d/2; Btilda = B(dh+1:N-dh,dh+1:N-dh); else # d is odd dh = ceil(d/2); Btilda = B(dh:N-dh,dh:N-dh); endif W = W*B; U = Btilda; else ## W = W*speye(Nhat); U = speye(Nhat-d); endif ## Smoothing delta = trace(D'*D)/Nhat^(2+d); # using "relative" or other weighting affects this! yhat = (M'*W*M + lambda*delta^(-1)*D'*U*D) \ M'*W*y; #[R,P,S] = splchol(M'*W*M + lambda*delta^(-1)*D'*U*D); #yhat = S*(R'\(R\(S'*M'*W*y))); ## Computation of hat diagonal and cross-validation if (nargout > 1) ## from AIChE J. (2006) 52, 325 ## note: chol factorization does not help speed up the computation of H; ## should implement Eiler's partial H computation if many point smoothing by GCV is needed ##H = M*(S*(R'\(R\(S'*M'*W)))); H = M*((M'*W*M + lambda*delta^(-1)*D'*U*D)\M'*W); ## note: this is variance, squared of the standard error that Eilers uses v = (M*yhat - y)'*(M*yhat - y)/N / (1 - trace(H)/N)^2; endif endfunction

Generated by Doxygen 1.6.0 Back to index