📄 csapsdem.html
字号:
<html xmlns:mwsh="http://www.mathworks.com/namespace/mcode/v1/syntaxhighlight.dtd"> <head> <meta http-equiv="Content-Type" content="text/html; charset=utf-8"> <!--This HTML is auto-generated from an M-file.To make changes, update the M-file and republish this document. --> <title>Cubic smoothing spline</title> <meta name="generator" content="MATLAB 7.1"> <meta name="date" content="2005-07-27"> <meta name="m-file" content="csapsdem"> <link rel="stylesheet" type="text/css" href="../../matlab/demos/private/style.css"> </head> <body> <div class="header"> <div class="left"><a href="matlab:edit csapsdem">Open csapsdem.m in the Editor</a></div> <div class="right"><a href="matlab:echodemo csapsdem">Run in the Command Window</a></div> </div> <div class="content"> <h1>Cubic smoothing spline</h1> <introduction> <p>Illustration of the use of CSAPS and SPAPS.</p> </introduction> <h2>Contents</h2> <div> <ul> <li><a href="#1">CSAPS</a></li> <li><a href="#2">Example: noisy data from a cubic polynomial</a></li> <li><a href="#8">SPAPS</a></li> <li><a href="#9">Tolerance vs smoothing parameter</a></li> <li><a href="#11">CSAPS vs SPAPS</a></li> </ul> </div> <h2>CSAPS<a name="1"></a></h2> <p>The spline toolbox command 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, an abbreviated version of the online help: </p><pre class="codeinput"><span class="comment">%CSAPS Cubic smoothing spline.</span><span class="comment">%</span><span class="comment">% VALUES = CSAPS( X, Y, P, XX)</span><span class="comment">%</span><span class="comment">% Returns the values at XX of the cubic smoothing spline for the</span><span class="comment">% given data (X,Y) and depending on the smoothing parameter P, chosen from</span><span class="comment">% the interval [0 .. 1] . This smoothing spline f minimizes</span><span class="comment">%</span><span class="comment">% P * sum_i W(i)(Y(i) - f(X(i)))^2 + (1-P) * integral (D^2 f)^2</span><span class="comment">%</span></pre><h2>Example: noisy data from a cubic polynomial<a name="2"></a></h2> <p>Here are some trial runs. We start with data from a simple cubic, q(x) := x^3, contaminate the values with some noise, and choose the value of the smoothing parameter to be .5 , and plot the resulting smoothed values (o), along with the underlying cubic (:), and the contaminated data (x): </p><pre class="codeinput">q = inline(<span class="string">'x.^3'</span>);xi = (0:.05:1); yi= q(xi);ybad = yi+.3*(rand(size(xi))-.5);ys = csaps(xi,ybad,.5,xi);plot(xi,yi,<span class="string">':'</span>,xi,ybad,<span class="string">'x'</span>,xi,ys,<span class="string">'ro'</span>),grid <span class="string">off</span>title(<span class="string">'Clean data (:), noisy data (x), smoothed values (o)'</span>)legend(<span class="string">'exact'</span>,<span class="string">'noisy'</span>,<span class="string">'smoothed'</span>,2)</pre><img vspace="5" hspace="5" src="csapsdem_01.png"> <p>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. </p> <p>We try p = .6, .7, .8, .9, 1 , and plot the resulting smoothing splines.</p> <p>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. </p><pre class="codeinput">xxi = (0:100)/100;hold <span class="string">on</span>nx=length(xxi);yy=zeros(5,nx);<span class="keyword">for</span> j=1:5 yy(j,:) = csaps(xi,ybad,.5+j/10,xxi);<span class="keyword">end</span>plot(xxi,yy)title(<span class="string">'Smoothing splines for various values of smoothing parameter'</span>)legend(<span class="string">'exact'</span>,<span class="string">'noisy'</span>,<span class="string">'p= .5'</span>,<span class="string">'p= .6'</span>,<span class="string">'p= .7'</span>,<span class="string">'p= .8'</span>,<span class="string">'p= .9'</span>,<span class="string">'p=1'</span>,2)hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="csapsdem_02.png"> <p>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 neighboring 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) . </p> <p>The next slide shows the smoothing spline for values of p near this magic number 1/(1+epsilon).</p> <p>For this case, it is more informative to look at 1-p since the magic number, 1/(1+epsilon), is very close to 1; see the last calculation below. </p><pre class="codeinput">epsilon = max(diff(xi))^3/16;plot(xi,yi,<span class="string">':'</span>,xi,ybad,<span class="string">'x'</span>), grid <span class="string">off</span>hold <span class="string">on</span>, labels = cell(1,5);<span class="keyword">for</span> j=1:5 p = 1/(1+epsilon*3^(j-3)); yy(j,:) = csaps(xi,ybad,p,xxi); labels{j} = [<span class="string">'1-p= '</span>,num2str(1-p)];<span class="keyword">end</span>plot(xxi,yy)title(<span class="string">'Smoothing splines for smoothing parameter near its ''magic'' value'</span>)legend(<span class="string">'exact'</span>,<span class="string">'noisy'</span>,labels{1},labels{2},labels{3},labels{4},labels{5},2)hold <span class="string">off</span>1 - 1/(1+epsilon)</pre><pre class="codeoutput">ans = 7.8124e-006</pre><img vspace="5" hspace="5" src="csapsdem_03.png"> <p>In this example, choosing the smoothing parameter near the magic number seems to give the best results. To be sure, the noise is so large here that we have no hope of recovering the exact data to any accuracy. </p> <p>You can also supply CSAPS with error weights, to pay more attention to some data points than others. Also, if you do not supply the evaluation sites XX, then CSAPS returns the ppform of the smoothing spline. </p> <p>Finally, CSAPS can also handle vector-valued data and even multivariable, gridded data.</p> <h2>SPAPS<a name="8"></a></h2> <p>The cubic smoothing spline provided by the spline toolbox command SPAPS differs from the one constructed in CSAPS only in the way it is selected. </p> <p>Here is an abbreviated version of the online help for SPAPS:</p><pre class="codeinput"><span class="comment">%SPAPS Smoothing spline.</span><span class="comment">%</span><span class="comment">% [SP,VALUES] = SPAPS(X,Y,TOL) returns the B-form and, if asked, the</span><span class="comment">% values at X , of a cubic smoothing spline f for the given data</span><span class="comment">% (X(i),Y(:,i)), i=1,2, ..., n .</span><span class="comment">%</span><span class="comment">% The smoothing spline f minimizes the roughness measure</span><span class="comment">%</span><span class="comment">% F(D^2 f) := integral ( D^2 f(t) )^2 dt on X(1) < t < X(n)</span><span class="comment">%</span><span class="comment">% over all functions f for which the error measure</span><span class="comment">%</span><span class="comment">% E(f) := sum_j { W(j)*( Y(:,j) - f(X(j)) )^2 : j=1,...,n }</span><span class="comment">%</span><span class="comment">% is no bigger than the given TOL.</span><span class="comment">% Here, D^M f denotes the M-th derivative of f .</span><span class="comment">% The weights W are chosen so that E(f) is the Composite Trapezoid Rule</span><span class="comment">% approximation for F(y-f).</span><span class="comment">%</span><span class="comment">% f is constructed as the unique minimizer of</span><span class="comment">%</span><span class="comment">% rho*E(f) + F(D^2 f) ,</span><span class="comment">%</span><span class="comment">% with the smoothing parameter RHO so chosen so that E(f) equals TOL .</span><span class="comment">% Hence, FN2FM(SP,'pp') should be (up to roundoff) the same as the output</span><span class="comment">% from CPAPS(X,Y,RHO/(1+RHO)).</span>cla, title(<span class="string">''</span>), legend <span class="string">off</span>, axis <span class="string">off</span></pre><img vspace="5" hspace="5" src="csapsdem_04.png"> <h2>Tolerance vs smoothing parameter<a name="9"></a></h2> <p>It may be easier to supply a suitable tolerance than the smoothing parameter P required by CSAPS. In our earlier example, we added uniformly distributed random noise from the interval 0.3*[-0.5 .. 0.5]. Hence we can estimate a reasonable value for TOL as the value of the error measure E at such noise. The graph shows the resulting smoothing spline constructed by SPAPS. Note that I specified the error weights to be uniform, which is their default value in CSAPS. </p><pre class="codeinput">tol = sum((.3*(rand(size(yi))-.5)).^2);[sp,ys] = spaps(xi,ybad,tol,ones(size(xi)));plot(xi,yi,<span class="string">':'</span>,xi,ybad,<span class="string">'x'</span>,xi,ys,<span class="string">'ro'</span>),grid <span class="string">off</span>title(<span class="string">'Clean data (:), noisy data (x), smoothed values (o)'</span>)legend(<span class="string">'exact'</span>,<span class="string">'noisy'</span>,<span class="string">'smoothed'</span>,2)</pre><img vspace="5" hspace="5" src="csapsdem_05.png"> <p>Here, in addition, is also the smoothing spline provided by CSAPS when not given a smoothing parameter, in which case the parameter is chosen by a certain ad hoc procedure that attempts to locate the region where the smoothing spline is most sensitive to the smoothing parameter. </p><pre class="codeinput">hold <span class="string">on</span>plot(xi,fnval(csaps(xi,ybad),xi),<span class="string">'*'</span>)legend(<span class="string">'exact'</span>,<span class="string">'noisy'</span>,<span class="string">'from SPAPS, specified TOL'</span>,<span class="string">'from CSAPS, default P'</span>,2)xlabel(<span class="string">' ... also the values (*) smoothed by CSAPS using default parameter.'</span>)hold <span class="string">off</span></pre><img vspace="5" hspace="5" src="csapsdem_06.png"> <h2>CSAPS vs SPAPS<a name="11"></a></h2> <p>In addition to the different ways, smoothing parameter vs tolerance, in which you specify a particular smoothing spline, these two commands differ also in that SPAPS, in addition to the cubic smoothing spline, can also provide a linear or a quintic smoothing spline. </p> <p>The quintic smoothing spline is better than the cubic smoothing spline in the situation when you would like the second derivative to move as little as possible. </p> <p class="footer">Copyright 1987-2005 C. de Boor and The MathWorks, Inc.<br> Published with MATLAB® 7.1<br></p> </div> <!--##### SOURCE BEGIN #####%% Cubic smoothing spline%% Illustration of the use of CSAPS and SPAPS.% Copyright 1987-2005 C. de Boor and The MathWorks, Inc.% $Revision: 1.18.4.2 $⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -