## [77d7c2]: inst / tkrgscatdatasmooth.m  Maximize  Restore  History

### 155 lines (140 with data), 5.1 kB

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154``` ```%%[x, y, lambda] = tkrgscatdatasmooth (xm, ym) %%[x, y, lambda] = tkrgscatdatasmooth (xm, ym, N, d) %%[x, y, lambda] = tkrgscatdatasmooth (xm, ym, N, d, range) %%[x, y, lambda] = tkrgscatdatasmooth (xm, ym, N, d, range, option, value) %% %% Determines a smooth curve that approximates the scattered (xm,ym) %% data values by Tikhonov regularization. The number of points 'N' for %% the smooth curve and the order of the smoothing derivative 'd' can be %% provided (defaults are 100 and 2 respectively). Additionally, the %% desired output range for x ([min,max]) can be given; if the provided %% range does not completely span the range of the data, the range %% defaults to the min and max of the data. The option-value pair %% should be either the regularizaiton parameter "lambda" or the %% standard deviation "stdev" of the randomness in the data. With no %% option supplied, generalized cross-validation is used to determine %% lambda. %% Reference: Anal. Chem. (2003) 75, 3631. %% See also: datasmooth %% %% 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 . function [x, y, lambda] = tkrgscatdatasmooth (xm, ym, N, d, range, option, value) if (length(xm)!=length(ym)) error("xm and ym must be equal length vectors") endif if ( size(xm)(1)==1 ) xm = xm'; endif if ( size(ym)(1)==1 ) ym = ym'; endif if (nargin < 3) N = 100; d = 2; endif if (nargin < 5) range = [min(xm),max(xm)] endif guess = 0; if (nargin > 5) if ( strcmp(option,"lambda") ) %% if lambda provided, use it directly lambda = value; elseif ( strcmp(option,"stdev") ) %% if stdev is provided, scale it and use it stdev = value; opt = optimset("TolFun",1e-6,"MaxFunEvals",20); log10lambda = fminunc ("tkrgdatasmscatwrap", guess, opt, xm, ym, N, d, range, "stdev", stdev); lambda = 10^log10lambda; else warning("option %s is not recognized; using cross-validation",option) endif else %% otherwise, perform cross-validation opt = optimset("TolFun",1e-4,"MaxFunEvals",20); log10lambda = fminunc ("tkrgdatasmscatwrap", guess, opt, xm, ym, N, d, range, "cve"); lambda = 10^log10lambda; endif [x,y] = tkrgdatasmscat (xm, ym, lambda, N, d, range); endfunction %!demo %! npts = 80; %! xm = linspace(0,1,npts)'; %! stdev = 1e-1; %! xm = xm + stdev*randn(npts,1); %! ym = sin(10*xm); %! ym = ym + stdev*randn(npts,1); %! ymp = ddmat(xm,1)*ym; %! ym2p = ddmat(xm,2)*ym; %! [x, y, lambda] = tkrgscatdatasmooth (xm,ym,500,4,[-0.15,1.15]); %! lambda %! yp = ddmat(x,1)*y; %! y2p = ddmat(x,2)*y; %! figure(1); %! plot(xm,ym,'o',x,y) %! figure(2); %! plot(xm(1:end-1),ymp,'o',x(1:end-1),yp) %! axis([min(x),max(x),min(yp)-abs(min(yp)),max(yp)*2]) %! figure(3) %! plot(xm(2:end-1),ym2p,'o',x(2:end-1),y2p) %! axis([min(x),max(x),min(y2p)-abs(min(y2p)),max(y2p)*2]) %! %-------------------------------------------------------- %! % this demo used generalized cross-validation to determine lambda %!demo %! npts = 80; %! xm = linspace(0,1,npts)'; %! stdev = 1e-1; %! xm = xm + stdev*randn(npts,1); %! ym = sin(10*xm); %! ym = ym + stdev*randn(npts,1); %! ymp = ddmat(xm,1)*ym; %! ym2p = ddmat(xm,2)*ym; %! [x, y, lambda] = tkrgscatdatasmooth (xm,ym,500,4,[-0.15,1.15],"stdev",stdev); %! lambda %! yp = ddmat(x,1)*y; %! y2p = ddmat(x,2)*y; %! figure(1); %! plot(xm,ym,'o',x,y) %! figure(2); %! plot(xm(1:end-1),ymp,'o',x(1:end-1),yp) %! axis([min(x),max(x),min(yp)-abs(min(yp)),max(yp)*2]) %! figure(3) %! plot(xm(2:end-1),ym2p,'o',x(2:end-1),y2p) %! axis([min(x),max(x),min(y2p)-abs(min(y2p)),max(y2p)*2]) %! %-------------------------------------------------------- %! % this demo used standard deviation to determine lambda %!demo %! npts = 80; %! xm = linspace(0,1,npts)'; %! stdev = 1e-1; %! xm = xm + stdev*randn(npts,1); %! ym = sin(10*xm); %! ym = ym + stdev*randn(npts,1); %! ymp = ddmat(xm,1)*ym; %! ym2p = ddmat(xm,2)*ym; %! [x, y, lambda] = tkrgscatdatasmooth (xm,ym,500,4,[-0.15,1.15],"lambda",10000); %! lambda %! yp = ddmat(x,1)*y; %! y2p = ddmat(x,2)*y; %! figure(1); %! plot(xm,ym,'o',x,y) %! figure(2); %! plot(xm(1:end-1),ymp,'o',x(1:end-1),yp) %! axis([min(x),max(x),min(yp)-abs(min(yp)),max(yp)*2]) %! figure(3) %! plot(xm(2:end-1),ym2p,'o',x(2:end-1),y2p) %! axis([min(x),max(x),min(y2p)-abs(min(y2p)),max(y2p)*2]) %! %-------------------------------------------------------- %! % this demo used a user specified lambda that was too large ```