Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo

Close

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

Download this file

tkrgdatasmscat.m    90 lines (77 with data), 2.9 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
%%[x, y] = tkrgdatasmscat (xm, ym, lambda)
%%[x, y] = tkrgdatasmscat (xm, ym, lambda, N, d)
%%[x, y] = tkrgdatasmscat (xm, ym, lambda, N, d, range)
%%[x, y, v] = tkrgdatasmscat (...)
%%
%% Smooths the data of the scattered y vs. x values by
%% Tikhonov regularization.
%% Input:
%% xm data series x
%% ym data series y
%% lambda: smoothing parameter
%% N: number of equally spaced smoothed points (default = 100)
%% d: order of smoothing derivative (default = 2)
%% range: two element vector [min,max] giving the desired output range of x;
%% if the range is less than the data, defaults to [min,max] of the data
%% Output:
%% x: smoothed x
%% y: smoothed y
%% v: generalized cross-validation variance
%%
%% References: Anal. Chem. (2003) 75, 3631; AIChE J. (2006) 52, 325
%%
%% Some code addapted with permission from 'whitsmscat', Paul Eilers, 2003
%% 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/>.
function [x, y, v] = tkrgdatasmscat (xm, ym, lambda, N, d, range)
%% Defaults if not provided
if (nargin <= 3)
N = 100;
d = 2;
endif
%% sort the scattered data
[xm,idx] = sort (xm);
ym = ym(idx);
%% construct x
Nm = length(xm);
if (nargin > 5)
xmin = min(range(1),xm(1));
xmax = max(range(2),xm(end));
else
xmin = xm(1);
xmax = xm(end);
endif
L = xmax - xmin;
dx = L/(N-1);
x = (xmin:dx:xmax)';
%% Itilda is the mapping matrix y->ym
Itilda = speye(N);
idx = round((xm - xmin) / dx) + 1;
Itilda = Itilda(idx,:);
%% D is the derivative matrix
%%D = ddmat(x,d);
D = dx^(-d)*diff(speye(N),d); % equivalent but a little faster
%% Smoothing
delta = sqrt(trace(D'*D)); % a suitable invariant of D
%%y = (Itilda'*Itilda + lambda*D'*D) \ Itilda'*ym;
C = chol(Itilda'*Itilda + lambda*delta^(-1)*D'*D);
y = C \ (C' \ Itilda'*ym);
%% Computation of hat diagonal and cross-validation
if (nargout > 2)
%% from AIChE J. (2006) 52, 325
H = Itilda * ( C \ (C' \ Itilda') );
%% note: this is variance, squared of the standard error that Eilers uses
v = (ym - Itilda*y)'*(ym - Itilda*y)/Nm / (1 - trace(H)/Nm)^2;
endif
endfunction