% CSAPSDEM	Demonstrate cubic smoothing spline.
% Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
                                         % latest change: September 28, 1990
clc;clg;echo on;home
% The M-file  csaps  provides the  s m o o t h i n g  spline. This is a 
% cubic spline which 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:
help csaps
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(1,length(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,'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:

hold on
nx=length(xi);yy=zeros(5,nx);
for j=1:5
   yy(j,:) = csaps(xi,ybad,.5+j/10,xi);
end
plot(xi,yy);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 
% abscissae. 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')
hold on
for j=1:5
   yy(j,:) = csaps(xi,ybad,1/(1+epsilon*3^(j-3)),xi);
end
plot(xi,yy);pause
hold off

% We are not disappointed.

pause                                             % Touch any key to finish
