Home > wafo > misc > smooth.m

# smooth

## PURPOSE

Calculates a smoothing spline.

## SYNOPSIS

[yy,coefs]= smooth(x,y,p,xx,LinExtrap,d2)

## DESCRIPTION

``` SMOOTH Calculates a smoothing spline.

CALL:  yy = smooth(x,y,p,xx,def,v)

x  = x-coordinates of data. (vector)
y  = y-coordinates of data. (vector)
p  = [0...1] is a smoothing parameter:
0 -> LS-straight line
1 -> cubic spline interpolant
xx = the x-coordinates in which to calculate the smoothed function.
yy = the calculated y-coordinates of the smoothed function if
xx is given. If xx is not given then yy contain the
pp-form of the spline, for use with PPVAL.
def = 0 regular smoothing spline (default)
1 a smoothing spline with a constraint on the ends to
ensure linear extrapolation outside the range of the data
v = variance of each y(i) (default  ones(length(X),1))

Given the approximate values

y(i) = g(x(i))+e(i)

of some smooth function, g, where e(i) is the error. SMOOTH tries to
recover g from y by constructing a function, f, which  minimizes

p * sum (Y(i) - f(X(i)))^2/d2(i)  +  (1-p) * int (f'')^2

The call  pp = smooth(x,y,p)  gives the pp-form of the spline,
for use with PPVAL.

Example:%
x = linspace(0,1).';
y = exp(x)+1e-1*randn(size(x));
pp = smooth(x,y,.9);
plot(x,y,x,smooth(x,y,.99,x),'g',x,ppval(pp,x),'k',x,exp(x),'r')

## CROSS-REFERENCE INFORMATION

This function calls:
 diff Difference and approximate derivative. error Display message and abort function. mkpp Make piecewise polynomial. ppval Evaluate piecewise polynomial. spdiags Sparse matrix formed from diagonals.
This function is called by:
 cdf2tr Estimate transformation, g, from observed CDF. dat2tr Estimate transformation, g, from data. dist2dsmfun Smooths the conditional DIST2D distribution parameters. dist2dsmfun2 Smooths the conditional DIST2D distribution parameters. hermitefun Calculates the transformation by a Hermite polynomial. jhsnlpdf Joint (Scf,Hd) PDF for nonlinear waves with a JONSWAP spectra. jhspdf Joint (Scf,Hd) PDF for linear waves with JONSWAP spectra. jhvnlpdf Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum. jhvpdf Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum. lc2tr Estimate transformation, g, from observed crossing intensity. lc2tr2 Estimate transformation, g, from observed crossing intensity, version2. ohhgparfun Wave height, Hd, distribution parameters for Ochi-Hubble spectra. ohhspdf Joint (Scf,Hd) PDF for linear waves with Ochi-Hubble spectra. ohhsspdf Joint (Scf,Hd) PDF linear waves in space with Ochi-Hubble spectra. ohhvpdf Joint (Vcf,Hd) PDF for lineare waves with Ochi-Hubble spectra. th2vhpdf Transform joint T-H density to V-H density thsnlpdf Joint (Scf,Hd) PDF for nonlinear waves with Torsethaugen spectra. thspdf Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra. thspdf2 Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra. thsspdf Joint (Scf,Hd) PDF for linear waves in space with Torsethaugen spectra. thvpdf Joint (Vcf,Hd) PDF for linear waves with Torsethaugen spectra. wafofig5 Joint distribution (pdf) of crest front velocity and wave height:

## SOURCE CODE

```001 function [yy,coefs]= smooth(x,y,p,xx,LinExtrap,d2)
002 %SMOOTH Calculates a smoothing spline.
003 %
004 %  CALL:  yy = smooth(x,y,p,xx,def,v)
005 %
006 %         x  = x-coordinates of data. (vector)
007 %         y  = y-coordinates of data. (vector)
008 %         p  = [0...1] is a smoothing parameter:
009 %              0 -> LS-straight line
010 %              1 -> cubic spline interpolant
011 %         xx = the x-coordinates in which to calculate the smoothed function.
012 %         yy = the calculated y-coordinates of the smoothed function if
013 %              xx is given. If xx is not given then yy contain the
014 %              pp-form of the spline, for use with PPVAL.
015 %        def = 0 regular smoothing spline (default)
016 %              1 a smoothing spline with a constraint on the ends to
017 %                ensure linear extrapolation outside the range of the data
018 %          v = variance of each y(i) (default  ones(length(X),1))
019 %
020 %  Given the approximate values
021 %
022 %                y(i) = g(x(i))+e(i)
023 %
024 %  of some smooth function, g, where e(i) is the error. SMOOTH tries to
025 %  recover g from y by constructing a function, f, which  minimizes
026 %
027 %              p * sum (Y(i) - f(X(i)))^2/d2(i)  +  (1-p) * int (f'')^2
028 %
029 %  The call  pp = smooth(x,y,p)  gives the pp-form of the spline,
030 %  for use with PPVAL.
031 %
032 % Example:%
033 %  x = linspace(0,1).';
034 %  y = exp(x)+1e-1*randn(size(x));
035 %  pp = smooth(x,y,.9);
036 %  plot(x,y,x,smooth(x,y,.99,x),'g',x,ppval(pp,x),'k',x,exp(x),'r')
037 %
038 % See also  lc2tr, dat2tr, ppval
039
040
041 %References:
042 % Carl de Boor (1978)
043 % 'Practical Guide to Splines'
044 %  Springer Verlag
045 %  Uses EqXIV.6--9, pp 239
046
047 % tested on: Matlab 4.x 5.x
048 %History:
049 % revised pab  26.11.2000
050 % - added example
051 % - fixed a bug: n=2 is now possible
052 % revised by pab 21.09.99
053 %    - added d2
054 % revised by pab 29.08.99
055 %    -new extrapolation: ensuring that
056 %     the smoothed function has contionous derivatives
057 %     in the first and last knot
058 % modified by Per A. Brodtkorb  23.09.98
059 %  secret option forcing linear extrapolation outside the ends when p>0
060 %  used in lc2tr
061
062 if (nargin<5)|(isempty(LinExtrap)),
063   LinExtrap=0; %do not force linear extrapolation in the ends (default)
064 end
065
066
067 [xi,ind]=sort(x(:));
068 n = length(xi);
069 y = y(:);
070 if n<2,
071    error('There must be >=2 data points.')
072 elseif any(diff(xi)<=0),
073    error('Two consecutive values in x can not be equal.')
074 elseif n~=length(y),
075    error('x and y must have the same length.')
076 end
077
078 if nargin<6|isempty(d2),
079   d2 = ones(n,1);  %not implemented yet
080 else
081   d2=d2(:);
082 end
083
084
085 yi=y(ind); %yi=yi(:);
086
087 dx = diff(xi);
088 dydx = diff(yi)./dx;
089
090 if (n==2)  % straight line
091   coefs=[dydx yi(1)];
092 else
093   if LinExtrap & n==3, p = 0;end % Force LS-fit
094   dx1=1./dx;
095   Q = spdiags([dx1(1:n-2) -(dx1(1:n-2)+dx1(2:n-1)) dx1(2:n-1)],0:-1:-2,n,n-2);
096   D = spdiags(d2,0,n,n);  % The variance
097   R = spdiags([dx(1:n-2) 2*(dx(1:n-2)+dx(2:n-1)) dx(2:n-1)],-1:1,n-2,n-2);
098
099   QQ = (6*(1-p))*(Q.'*D*Q)+p*R;
100   % Make sure Matlab uses symmetric matrix solver
101   u  = 2*((QQ+QQ')\diff(dydx));                     % faster than u=QQ\(Q'*yi);
102   ai = yi-6*(1-p)*D*diff([0;diff([0;u;0]).*dx1;0]); % faster than yi-6*(1-p)*Q*u
103
104   % The piecewise polynominals are written as
105   % fi=ai+bi*(x-xi)+ci*(x-xi)^2+di*(x-xi)^3
106   % where the derivatives in the knots according to Carl de Boor are:
107   %    ddfi  = 6*p*[0;u] = 2*ci;
108   %    dddfi = 2*diff([ci;0])./dx = 6*di;
109   %    dfi   = diff(ai)./dx-(ci+di.*dx).*dx = bi;
110
111   ci = 3*p*[0;u];
112
113
114   if LinExtrap & p~=0& n>3, %Forcing linear extrapolation in the ends
115     ci([2,  end]) = 0;
116   % New call
117   % fixing the coefficients so that we have continous
118   % derivatives everywhere
119     ai(1) = -(ai(3)-ai(2))*dx(1)/dx(2) +ai(2)+ ci(3)*dx(1)*dx(2)/3;
120     ai(n) = (ai(n-1)-ai(n-2))*dx(n-1)/dx(n-2) +ai(n-1)+ ci(n-2)*dx(n-2)*dx(n-1)/3;
121   end
122
123   di    = diff([ci;0]).*dx1/3;
124   bi    = diff(ai).*dx1-(ci+di.*dx).*dx;
125   coefs = [di ci bi ai(1:n-1)];
126 end
127
128 pp = mkpp(xi,coefs);
129 if (nargin<4)|(isempty(xx)),
130   yy = pp;
131 else
132   yy = ppval(pp,xx);
133 end
134
135
136 return
137 % old call: linear extrapolation
138 % crude the derivatives are not necessarily continous
139 % where we force ci and di to zero
140 if LinExtrap, %forcing linear extrapolation in the ends
141   ci([2,  end])=0;  % could be done better
142   di([1 end])=0;
143 end
144 % Old call for LinExtrap
145 if 0 %LinExtrap & n>6 & (p~=0), %forcing linear extrapolation in the ends
146     Q = Q(2:end-1,2:end-1);
147     D = D(2:end-1,2:end-1);
148     R = R(2:end-1,2:end-1);
149     % new call : so ci(1:2)=ci(n-1:n)=di(1)=di(n)=0
150
151     QQ=(6*(1-p))*(Q.'*D*Q)+p*R;
152     % Make sure Matlab uses symmetric matrix solver
153     u=2*((QQ+QQ')\diff(dydx(2:n-2))); %faster than 2*(QQ+QQ.'')\(Q'*yi);
154
155     % The piecewise polynominals are written as
156     % fi=ai+bi*(x-xi)+ci*(x-xi)^2+di*(x-xi)^3
157     ai=yi(2:n-1)-(6*(1-p)*D*diff([0;diff([0;u;0])./dx(2:n-2);0]));%faster than yi-6*(1-p)*D*Q*u ;
158     ci=3*p*[0;0;u;0];
159     % fixing the coefficients so that we have continous
160     % derivatives everywhere
161     a1=-(ai(2)-ai(1))*dx(1)/dx(2) +ai(1)+ ci(3)*dx(1)*dx(2)/3;
162     an=(ai(n-2)-ai(n-3))*dx(n-1)/dx(n-2) +ai(n-2)+ ci(n-2)*dx(n-2)*dx(n-1)/3;
163     ai=[a1;ai; an];
164
165   else
166
167   end
168
169
170```

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group

Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003