home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 5.ddi / SPLINES.DI$ / CSAPSDEM.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.6 KB  |  69 lines
  1. % CSAPSDEM    Demonstrate cubic smoothing spline.
  2. % Copyright (c) 1990-92 by Carl de Boor and The MathWorks, Inc.
  3.                                          % latest change: September 28, 1990
  4. clc;clg;echo on;home
  5. % The M-file  csaps  provides the  s m o o t h i n g  spline. This is a 
  6. % cubic spline which more or less follows the presumed underlying trend in 
  7. % noisy data.  A smoothing parameter, to be chosen by you, determines just
  8. % how closely the smoothing spline follows the given data. 
  9. % Here is the basic information:
  10. help csaps
  11. pause                                             % Touch any key to continue
  12. %      Here are some trial runs.
  13.  
  14. %  We start with a simple cubic
  15. xi=[0:.05:1];
  16. yi=xi.*xi.*xi;
  17. %  .. and contaminate the data with some random noise:
  18. ybad=yi+.3*(rand(1,length(xi))-.5);
  19.  
  20. % The smoothing requires choice of a particular value for the smoothing 
  21. % parameter  p , between 0 and 1. We choose  .5 :
  22. ys = csaps(xi,ybad,.5,xi);
  23. % and plot clean data (.), contaminated data (X), and smoothed values (O):
  24. plot(xi,yi,':',xi,ybad,'x',xi,ys,'o');pause
  25.  
  26.  
  27. %  The smoothing is way overdone here. By choosing the smoothing parameter
  28. %  p  closer to 1, we obtain a smoothing spline closer to the given data.
  29. % We try  p = .6, .7, .8, .9, 1 , and plot all the smoothed values:
  30.  
  31. hold on
  32. nx=length(xi);yy=zeros(5,nx);
  33. for j=1:5
  34.    yy(j,:) = csaps(xi,ybad,.5+j/10,xi);
  35. end
  36. plot(xi,yy);pause
  37. hold off
  38.  
  39. % We see that the smoothing spline is very sensitive to the choice of the
  40. % smoothing parameter. Even for  p =.9  , the smoothing spline is still 
  41. % far from the underlying trend, while for p=1 we get the interpolant to the
  42. % (noisy) data. 
  43.  
  44. % In fact, the formulation from p.235ff of PGS used here is very sensitive to
  45. % scaling of the independent variable. A simple analysis of the equations
  46. % used shows that the sensitive range for  p  is  around  1/(1+ epsilon) ,
  47. % with  epsilon := h^3/16 , and  h  the average difference between the given 
  48. % abscissae. Specifically, one would expect a close following of the data 
  49. % when  p = 1/(1+epsilon/10)  and some satisfactory smoothing when  
  50. %  p = 1/(1+epsilon*10) .
  51.  
  52. % In our case,
  53. epsilon = max(diff(xi))^3/16
  54. % Let's try values of  p  near this magic number  1 - epsilon. 
  55.  
  56. pause                                             % Touch any key to continue
  57.  
  58. plot(xi,yi,':',xi,ybad,'x')
  59. hold on
  60. for j=1:5
  61.    yy(j,:) = csaps(xi,ybad,1/(1+epsilon*3^(j-3)),xi);
  62. end
  63. plot(xi,yy);pause
  64. hold off
  65.  
  66. % We are not disappointed.
  67.  
  68. pause                                             % Touch any key to finish
  69.