e02bef

e02bef © Numerical Algorithms Group, 2002.

Purpose

E02BEF Least-squares cubic spline curve fit, automatic knot placement

Synopsis

[lamda,c,fp,wrk,iwrk,n,ifail] = e02bef(x,y,s,lamda<,start,w,wrk,iwrk,n,lwrk,...
ifail>)

Description

 
 This routine determines a smooth cubic spline approximation s(x) 
 to the set of data points (x ,y ), with weights w , for 
                             r  r                 r     
 r=1,2,...,m.
 
 The spline is given in the B-spline representation
 
                             n-4                             
                             --                              
                       s(x)= >  c N (x)                        (1)
                             --  i i                         
                             i=1                             
 
 where N (x) denotes the normalised cubic B-spline defined upon 
        i                                                      
 the knots (lambda) ,(lambda)   ,...,(lambda)   .
                   i         i+1             i+4
 
 The total number n of these knots and their values 
 (lambda) ,...,(lambda)  are chosen automatically by the routine. 
         1             n                                         
 The knots (lambda) ,...,(lambda)    are the interior knots; they 
                   5             n-4                             
 divide the approximation interval [x ,x ] into n-7 sub-intervals.
                                     1  m           
 The coefficients c ,c ,...,c    are then determined as the 
                   1  2      n-4                           
 solution of the following constrained minimization problem:
 
 minimize
 
                             n-4                             
                             --        2                     
                      (eta)= >  (delta)                        (2)
                             --        i                     
                             i=5                             
 
 subject to the constraint
 
                           m                                 
                           --          2                     
                  (theta)= >  (epsilon) <=S                    (3)
                           --          r                     
                           r=1                               
 
 where: (delta)    stands for the discontinuity jump in the third
               i   order derivative of s(x) at the interior knot
 		  (lambda) ,
                           i
 
        (epsilon)  denotes the weighted residual w (y -s(x )),
                 r                                r  r    r  
 
 and    S          is a non-negative number to be specified by 
                   the user.
 
 The quantity (eta) can be seen as a measure of the (lack of) 
 smoothness of s(x), while closeness of fit is measured through 
 (theta). By means of the parameter S, 'the smoothing factor', the
 user will then control the balance between these two (usually 
 conflicting) properties. If S is too large, the spline will be 
 too smooth and signal will be lost (underfit); if S is too small,
 the spline will pick up too much noise (overfit). In the extreme 
 cases the routine will return an interpolating spline ((theta)=0)
 if S is set to zero, and the weighted least-squares cubic 
 polynomial ((eta)=0) if S is set very large. Experimenting with S
 values between these two extremes should result in a good 
 compromise.
 
 Values of the computed spline, or of its derivatives or definite 
 integral, can subsequently be computed by calling E02BBF, E02BCF 
 or E02BDF.
 

Parameters

e02bef

Required Input Arguments:

x (:)                                 real
y (:)                                 real
s                                     real
lamda (:)                             real

Optional Input Arguments:                       <Default>

start (1)                             string   'c'
w (:)                                 real     ones(length(x),1)
wrk (lwrk)                            real     if lower(start)=='c';wrk=...
                                        ...    zeros(lwrk,1);end
iwrk (:)                              integer  if lower(start)=='c';iwrk=...
                                        ...    zeros(length(lamda),1);end
n                                     integer  if lower(start)=='c';n=0;end
lwrk                                  integer  e02bef13(length(x),...
                                        ...    length(lamda))
ifail                                 integer  -1

Output Arguments:

lamda (:)                             real
c (:)                                 real
fp                                    real
wrk (lwrk)                            real
iwrk (:)                              integer
n                                     integer
ifail                                 integer