1 ## Copyright (C) 2008 Jonathan Stickel <jonathan.stickel@nrel.gov>
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
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
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/>.
17 ##@deftypefn {Function File} {[@var{yhat}, @var{v}] =} rgdtsmcore (@var{x}, @var{y}, @var{d}, @var{lambda}, [@var{options}])
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.
27 ## Note: the options have changed!
28 ## Currently supported input options are (multiple options are allowed):
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.
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"
44 ## References: Anal. Chem. (2003) 75, 3631; AIChE J. (2006) 52, 325
45 ## @seealso{regdatasmooth}
49 function [yhat, v] = rgdtsmcore (x, y, d, lambda, varargin)
55 ## Defaults if not provided
62 ## parse the provided options
63 if ( length(varargin) )
64 for i = 1:length(varargin)
73 weightv = varargin{i+1};
79 printf("Option '%s' is not implemented;\n", arg)
84 if (xhatprov && midpr)
85 warning("midpointrule is currently not used if xhat is provided (since x,y may be scattered)")
88 if (weights && relative)
89 warning("relative differences is not used if a weighting vector is provided")
95 ## test that xhat is increasing
98 error("xhat must be monotonically increasing")
100 error("x must be monotonically increasing if xhat is not provided")
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")
110 idx = interp1(xhat,1:Nhat,x,"nearest"); # works for unequally spaced xhat
114 ## construct "weighting" matrices W and U
116 ## use arbitrary weighting as provided
119 ## use relative differences
120 Yinv = sparse(diag(1./y));
125 ## use midpoint rule integration (rather than simple sums)
127 Bhat = sparse(diag(-ones(N-1,1),-1)) + sparse(diag(ones(N-1,1),1));
130 B = 1/2*sparse(diag(Bhat*x));
131 if ( floor(d/2) == d/2 ) # test if d is even
133 Btilda = B(dh+1:N-dh,dh+1:N-dh);
136 Btilda = B(dh:N-dh,dh:N-dh);
141 ## W = W*speye(Nhat);
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)));
151 ## Computation of hat diagonal and cross-validation
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;
164 ##plot(x,y,'o',x,M*yhat,'x')