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

tkrgdatasmdd.m

%%ys = tkrgdatasmdd (x, y, lambda)
%%ys = tkrgdatasmdd (x, y, lambda, d)
%%[ys, v] = tkrgdatasmdd (...)
%%
%% Smooths the tkrgdata of y vs. x values by
%% Tikhonov regularization and divided differences (arbitrary spacing of x).
%% Input:
%%   x:      data series of sampling positions (must be increasing)
%%   y:      data series, assumed to be sampled at equal intervals
%%   lambda: smoothing parameter; large lambda gives smoother result
%%   d:      order of smoothing derivative
%% Output:
%%   ys:     smoothed data
%%   v:      generalized cross-validation variance
%%
%% References:  Anal. Chem. (2003) 75, 3631; AIChE J. (2006) 52, 325
%%

%% Addapted with permission from 'whitsmdd', 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 [ys, v] = tkrgdatasmdd(x, y, lambda, d)
  
  %% Defaults if not provided
  if nargin <= 3
    d = 2;
  endif

  %% find the average dx in order to scale lambda
  L = x(end) - x(1);
  N = length(x);
  dx = L/(N-1);
  
  %% form the matrices
  %% D is the derivative matrix
  D = ddmat(x,d);

  %% B and Btilda are total integration matrices
  Bhat = spdiag(-ones(N-1,1),-1) + spdiag(ones(N-1,1),1);
  Bhat(1,1) = -1;
  Bhat(N,N) = 1;
  B = 1/2*spdiag(Bhat*x);
  %%B = 1/dx*speye(N,N);  % force equal waiting, even for variable spaced x?
  if ( floor(d/2) == d/2 ) % test if d is even
    dh = d/2;
    Btilda = B(dh+1:N-dh,dh+1:N-dh);
  else % d is odd
    dh = ceil(d/2);
    Btilda = B(dh:N-dh,dh:N-dh);
  endif

  %% Smoothing
  delta = sqrt(trace(D'*D)); % a suitable invariant of D for scaling lambda
  %%f = (B + lambda*D'*Btilda*D) \ B*y;
  C = chol(B + lambda*delta^(-1)*D'*Btilda*D);
  ys = C \ (C' \ B*y);
  
  %% Computation of hat diagonal and cross-validation
  if (nargout > 1)
    %% from AIChE J. (2006) 52, 325
    %%H = (B + lambda*D'*Btilda*D) \ B;
    H = C \ (C' \ B);
    %% note: this is variance, squared of the standard error that Eilers uses
    v = (y - ys)'*B*(y - ys) / (1 - trace(H)/N)^2;
  endif
  
endfunction

Generated by  Doxygen 1.6.0   Back to index