Logo Search packages:      
Sourcecode: octave-data-smoothing version File versions  Download package

tkrgdatasmscat.m

%%[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

Generated by  Doxygen 1.6.0   Back to index