SPLINE and SPLINEC Calls

Chapter Contents

Language Reference

SPLINE and SPLINEC Calls

provide cubic spline fits

CALL SPLINE( fitted, data<, smooth><, delta><, nout>

<, type><, slope>);

CALL SPLINEC( fitted, coeff, endval, data<, smooth><, delta>

<, nout><, type><, slope>);

The SPLINEC subroutine returns the following values:

fitted: is an n ×2 matrix of fitted values.
coeff: is an n ×5 (or n ×9) matrix of spline coefficients. The matrix contains the cubic polynomial coefficients for the spline for each interval. Column 1 is the left endpoint of the x-interval for the regular (nonparametric) spline or the left endpoint of the parameter for the parametric spline. Columns 2-5 are the constant, linear, quadratic, and cubic coefficients, respectively, for the x-component. If a parametric spline is used, then columns 6-9 are the constant, linear, quadratic, and cubic coefficients, respectively, for the y-component. The coefficients for each interval are with respect to the variable x-x_i where x_i is the left endpoint of the interval and x is the point of interest. The matrix coeff can be processed to yield the integral or the derivative of the spline. This, in turn, can be used with the SPLINEV function to evaluate the resulting curves. The SPLINEC call returns coeff.
endval: is a 1 ×2 matrix of endpoint values. Slopes of the two ends of the curve are reported as angles expressed in degrees. The SPLINEC call returns the endval argument.

The inputs to the SPLINEC subroutine are as follows:

data

specifies a n ×2 (or n ×3) matrix of (x,y) points on which the spline is to be fit. The optional third column is used to specify a weight for each data point. If smooth >0, the weight column is used in calculations. A weight $\leq 0$ causes the data point to be ignored in calculations.

delta

is an optional scalar specifying the resolution constant. If delta is specified, the fitted points are spaced by the amount delta on the scale of the first column of data if a regular spline is used or on the scale of the curve length if a parametric spline is used. If both nout and delta are specified, nout is used and delta is ignored.

nout

is an optional scalar specifying the number of fitted points to be computed. The default is nout=200. If nout is specified, then nout equally spaced points are returned. The nout argument overrides the delta argument.

slope

is an optional 1 ×2 matrix of endpoint slopes given as angles in degrees. If a parametric spline is used, the angle values are used modulo 360. If a nonparametric spline is used, the tangent of the angles is used to set the slopes (that is, the effective angles range from -90 to 90 degrees).

smooth

is an optional scalar specifying the degree of smoothing to be used. If smooth is omitted or set equal to 0, then a cubic interpolating spline is fit to the data. If smooth >0, then a cubic spline is used. Larger values of smooth generate more smoothing.

type

is an optional 1 ×1 (or 1 ×2) character matrix or quoted literal giving the type of spline to be used. The first element of type should be one of the following:

periodic, which requests periodic endpoints
zero, which sets second derivatives at endpoints to 0

The type argument controls the endpoint constraints unless the slope argument is specified. If periodic is specified, the response values at the beginning and end of column 2 of data must be the same unless the smoothing spline is being used. If the values are not the same, an error message is printed and no spline is fit. The default value is zero. The second element of type should be one of the following:

nonparametric, which requests a nonparametric spline
parametric, which requests a parametric spline

If parametric is specified, a parameter sequence {t_i} is formed as follows: t₁=0 and

$t_i=t_{i-1} + \sqrt{ (x_i - x_{i-1})^2 + (y_i - y_{i-1})^2 }$

Splines are then fit to both the first and second columns of data. The resulting splined values are paired to form the output. Changing the relative scaling of the first two columns of data changes the output because the sequence {t_i} assumes Euclidean distance.

Note that if the points are not arranged in ascending order by the first columns of data, then a parametric method must be used. An error message results if the nonparametric spline is requested.

Refer to Stoer and Bulirsch (1980), Reinsch (1967), and Pizer (1975) for descriptions of the methods used to fit the spline.

   proc iml;
      data = { 0 1, 1 2, 2 3, 3 4, 4 5, 5 6, 6 7, 7 8, 8 9, 9 10 };
      call splinec(fitted,coeff,endval,data,,,,'zero',{45 45});
      v =  splinev(coeff,{-1 1 2 3 3.5 4 20});
      print v;
      v =  splinev(coeff,,3);
      print v;

This example takes a function defined by discrete data and finds the integral and the first moment of the function.

   data func;
      input x @@;
      y = x+0.1*sin(x);
      cards;
   0   2   5   7   8  10
   ;

   proc iml;

      use func;
      read all into a;
      call splinec(fit,coeff,endval,a,,,,'zero');

      start fcheck(x) global(coeff,pow);
         /************************************/
         /* Note that the first column of v  */
         /* contains the points of           */
         /* evaluation and the second column */
         /* contains the evaluation.         */           
         /************************************/
         v = x##pow # splinev(coeff,x //x)[1,2];
         return(v);
      finish;

      start moment(po) global(coeff,pow);
         pow = po;
         call quad(z,'fcheck',coeff[,1]) eps = 1.e-10;
         v1  = sum(z);
         return(v1);
      finish;

      mass  = moment(0);   /* to compute the mass */
      mass  = mass //
             (moment(1)/mass) // /* to compute the mean */
             (moment(2)/mass) ;  /* to compute the meansquare */
      print mass;

      /*****************************************************/
      /* Use Gauss-Legendre integration: this is not     */
      /* adaptive, but it is good for moments up to maxng. */
      /*****************************************************/

      gauss = {
          -9.3246951420315205e-01
          -6.6120938646626448e-01
          -2.3861918608319743e-01
           2.3861918608319713e-01
           6.6120938646626459e-01
           9.3246951420315183e-01,
           1.713244923791701e-01
           3.607615730481388e-01
           4.679139345726905e-01
           4.679139345726904e-01
           3.607615730481389e-01
           1.713244923791707e-01  };
      ngauss = ncol(gauss);
      maxng  =  2*ngauss-4;
      start moment1(pow) global(coeff,gauss,ngauss,maxng);
         if pow < maxng  then do;
            nrow    = nrow(coeff);
            ncol    = ncol(coeff);
            left    = coeff[1:nrow-1,1];
            right   = coeff[2:nrow,1];
            mid     = 0.5*(left+right);
            interv  = 0.5*(right - left);
          /* scale the weights on each interval */
            wgts    = shape(interv*gauss[2,],1);
          /* scale the points on each interval */
            pts     = shape(interv*gauss[1,] + mid * repeat(1,1,ngauss),1) ;
          /* evaluate the function */
            eval    = splinev(coeff,pts)[,2]`;
            mat     = sum (wgts#pts##pow#eval);
         end;
         return(mat);
      finish;

      mass    = moment1(0);       /* to compute the mass */
      mass    = mass // (moment1(1)/mass) //  (moment1(2)/mass) ;
      print mass;

Chapter Contents
Previous
Next
Top