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{lambda}] =} regdatasmooth (@var{x}, @var{y}, [@var{options}])
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.
24 ## Note: the options have changed!
25 ## Currently supported input options are (multiple options are allowed):
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
39 ## use generalized cross-validation to determine the optimal value for
40 ## lambda; if neither "lambda" nor "stdev" options are given, this
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.
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"
58 ## Please run the demos for example usage.
60 ## References: Anal. Chem. (2003) 75, 3631; AIChE J. (2006) 52, 325
61 ## @seealso{rgdtsmcorewrap, rgdtsmcore}
64 function [yhat, lambda] = regdatasmooth (x, y, varargin)
68 elseif ( length(x) != length(y) )
69 error("x and y must be equal length vectors")
71 if ( isrow(x) ) x = x'; endif
72 if ( isrow(y) ) y = y'; endif
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
93 lambda = varargin{i+1};
96 stdev = varargin{i+1};
101 guess = log10(varargin{i+1});
109 ## add warning if more than one gcv, lambda, or stdev options provided?
113 ## do nothing and use the provided lambda
115 ## find the "optimal" lambda
117 ## match standard deviation
118 fhandle = @(log10lambda) rgdtsmcorewrap (log10lambda, x, y, d, {"stdev", stdev}, options{:});
120 ## perform cross-validation
121 fhandle = @(log10lambda) rgdtsmcorewrap (log10lambda, x, y, d, {"cve"}, options{:});
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);
136 warning("Iteration limit of %i exceeded\n",maxiter)
138 lambda = 10^log10lambda;
141 yhat = rgdtsmcore (x, y, d, lambda, options{:});
147 %! x = linspace(0,2*pi,npts)';
148 %! x = x + 2*pi/npts*(rand(npts,1)-0.5);
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");
155 %! yhp = ddmat(x,1)*yh;
156 %! yh2p = ddmat(x,2)*yh;
159 %! plot(x,y,'o','markersize',5,x,yh,x,sin(x))
161 %! legend("noisy","smoothed","sin(x)","location","northeast");
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])
166 %! legend("noisy","smoothed","cos(x)","location","southeast");
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])
171 %! legend("noisy","smoothed","-sin(x)","location","southeast");
172 %! %--------------------------------------------------------
173 %! % smoothing of monotonic data, using "stdev" to determine the optimal lambda
177 %! x = rand(npts,1)*2*pi;
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);
185 %! plot(x,y,'o','markersize',10,xh,yh,xh,sin(xh))
187 %! legend("noisy","smoothed","sin(x)","location","northeast");
188 %! %--------------------------------------------------------
189 %! % smoothing of scattered data, using "gcv" to determine the optimal lambda