echo off
%CSAPSDEM Cubic smoothing spline demo.

%   Copyright 1987-2000 by C. de Boor and The MathWorks, Inc.
%   latest change: 23 oct 95; 28sep96; 8aug97
%   $Revision: 1.10 $

home, close, echo on
% The M-file  csaps  provides the  s m o o t h i n g  spline. This is a
% cubic spline that more or less follows the presumed underlying trend in
% noisy data.  A smoothing parameter, to be chosen by you, determines just
% how closely the smoothing spline follows the given data.
% Here is the basic information:

%CSAPS Cubic smoothing spline.
%
%   VALUES  = CSAPS( X, Y, P [, XX [, W ]])
%
%   Returns the values at XX of the cubic smoothing spline for the
%   given data (X,Y)  and depending on the smoothing parameter P from
%   [0 .. 1] .  This smoothing spline  f  minimizes
%
%   P * sum_i W(i)(Y(i) - f(X(i)))^2  +  (1-P) * integral (D^2 f)^2
%
pause                                             % Touch any key to continue
%      Here are some trial runs.

%  We start with a simple cubic
xi=[0:.05:1];
yi=xi.*xi.*xi;
%  .. and contaminate the data with some random noise:
ybad=yi+.3*(rand(size(xi))-.5);

% The smoothing requires choice of a particular value for the smoothing
% parameter  p , between 0 and 1. We choose  .5 :
ys = csaps(xi,ybad,.5,xi);
% and plot clean data (.), contaminated data (X), and smoothed values (O):
plot(xi,yi,':',xi,ybad,'x',xi,ys,'ro'),grid off
title('clean data (:), noisy data (x), smoothed values (o)'),pause


%  The smoothing is way overdone here. By choosing the smoothing parameter
%  p  closer to 1, we obtain a smoothing spline closer to the given data.
% We try  p = .6, .7, .8, .9, 1 , and plot all the smoothed values:

xxi = [0:100]/100;
hold on
nx=length(xxi);yy=zeros(5,nx);
for j=1:5
   yy(j,:) = csaps(xi,ybad,.5+j/10,xxi);
end
plot(xxi,yy)
title('smoothing splines for various values of smoothing parameter')
pause, hold off

% We see that the smoothing spline is very sensitive to the choice of the
% smoothing parameter. Even for  p =.9  , the smoothing spline is still
% far from the underlying trend, while for p=1 we get the interpolant to the
% (noisy) data.

% In fact, the formulation from p.235ff of PGS used here is very sensitive to
% scaling of the independent variable. A simple analysis of the equations
% used shows that the sensitive range for  p  is  around  1/(1+ epsilon) ,
% with  epsilon := h^3/16 , and  h  the average difference between the given
% sites. Specifically, one would expect a close following of the data
% when  p = 1/(1+epsilon/10)  and some satisfactory smoothing when
%  p = 1/(1+epsilon*10) .

% In our case,
epsilon = max(diff(xi))^3/16
% Let's try values of  p  near this magic number  1 - epsilon.

pause                                             % Touch any key to continue

plot(xi,yi,':',xi,ybad,'x'), grid off
hold on
for j=1:5
   yy(j,:) = csaps(xi,ybad,1/(1+epsilon*3^(j-3)),xxi);
end
plot(xxi,yy)
title('smoothing splines for smoothing parameter near its ''magic'' value')
pause, hold off

% We are not disappointed.

pause                                             % Touch any key to finish
echo off, close
