/*  SPLINE.SRC  - Cubic Spline Smoothing routines

**  David Baird   AgResearch PO Box 60 Lincoln New Zealand   25/7/1995

**  This code is freely available for public non-commercial use. Please

**  cite if used in published research. No performance guarantees are given.

**  Format                Purpose                                    Line

** =======================================================================

** c = csspline(x,y,w,p); Fits weighted cubic smoothing spline         34

** y = fspline(c,x);      y on cubic spline given x                    99

** x = inv_spl(c,y);      x on monotonic cubic spline given y         119

** d = dspline(c,x,k);    kth derivative of cubic spline              156

** a = ispline(c,x);      Integral of cubic spline                    186

** s = press(x,y,w,p);    PRESS statistic                             222

** s = xvsc(x,y,w,p);     Cross Validation (CV) statistic             307

** s = gxvsc(x,y,w,p);    Generalised CV (GCV) statistic              339

** s = agxvsc(x,y,w,p);   Aprroximate GCV (AGCV) statistic            409

** l = spl_clim(x,y,w,p); Confidence limits for spline                445

** p = min_prs(x,y,w);    Smoothing parameter which minimizes PRESS   537

** p = min_cv(x,y,w);     Smoothing parameter which minimizes CV      557

** p = min_gcv(x,y,w);    Smoothing parameter which minimizes GCV     582

** p = min_agcv(x,y,w);   Smoothing parameter which minimizes AGCV    607

** c = prspline(x,y,w);   Cubic spline which minimizes PRESS          632

** c = cvspline(x,y,w);   Cubic spline which minimizes CV             652

** c = gcspline(x,y,w);   Cubic spline which minimizes GCV            675

** c = agspline(x,y,w);   Cubic spline which minimizes AGCV           695

** call splset            Resets globals used in spline procedures    716

** ====================================================================*/

/* Declarations of Globals used in Spline Routines */

declare matrix _spl_c   != 0 ;    @ Approximation constants used by AGCV @

declare matrix _spl_r   != 0 ;    @ Linear regression design matrix @

declare matrix _spl_q   != 0 ;    @ Cubic spline design matrix @

declare matrix _spl_x   != 0 ;    @ Passes knots   to proc _spl_fh @

declare matrix _spl_w   != 0 ;    @ Passes weights to proc _spl_fh @

declare matrix _spl_tol != 1e-6 ; @ Controls precision of approx. constants @

/* CSSPLINE - Fit Cubic Smoothing Spline

**

** Usage:   c = csspline(x,y,w,p) ;

**

** Input:   x - N x 1 vector of x values of the set of points;

**              x must be sorted into ascending order and must not contain

**              duplicate values.

**          y - N x 1 vector of y values of the set of points.

**          w - N x 1 vector of weights for the set of points;

**              for equal weighting set w to a scalar for efficiency.

**          p - Scalar giving the smoothing parameter. This defines the

**              relative weighting of the residual sum of squares and the

**              integrated squared second derivative (a roughness penalty)

**              of the curve; p = 0 gives an interpolating cubic spline and

**              p = 1 gives a linear regression fit (0 <= p <= 1).

**

** Output:  c - N x 5 matrix of coefficients of cubic smoothing spline;

**              c[.,1] = knots,

**              c[.,2..5] = function value and first 3 derivatives at knots.

**

** Reference: C. de Boor 1978, "A Practical Guide to Splines", Springer-Verlag.

*/

proc csspline(x,y,w,p) ;

     local n,n1,n2,dx,dy,dxi,r,q,qw,qy,t,u,c0,c1,c2,c3,oldtrap ;

     n   = rows(x) ;

     w   = w./meanc(w) ;

     t   = 1 - p ;

     n1  = n - 1 ;

     n2  = n - 2 ;

     dx  = x[2:n] - x[1:n1] ;

     if not (dx > 0) ;

        errorlog("ERROR: CSSPLINE - x sorted into ascending order and must"\

                 " not contain duplicates") ;

        retp(error(1)) ;

     endif ;

     dy  = (y[2:n] - y[1:n1])./dx ;

     dxi = 1./dx ;

     @ Create R matrix which fits linear trend to y @

     r = rotater((0|dx[2:n2])~2.*(dx[1:n2]+dx[2:n1])~(dx[2:n2]|0)~

                    zeros(n2,n-5),seqa(-1,1,n2))./6 ;

     @ Create Q matrix which fits cubic spline to y @

     q = rotater(dxi[1:n2]~(-dxi[1:n2]-dxi[2:n1])~dxi[2:n1]~zeros(n2,n-3),

                 seqa(0,1,n2))' ;

     @ Solve for U where ((1-p)6(Q'W^(-1)Q + pR)U = Q'y @

     qw = q./sqrt(w) ;

     qw = t.*moment(qw,0) + p.*r ;

     qy = dy[2:n1] - dy[1:n2] ;

     oldtrap = trapchk(0xffff) ;

     trap 1 ;

     u = bandsolpd(qy,band(qw,2));

     trap oldtrap,0xffff ;

     if scalmiss(u) ;

        errorlog("ERROR: CSSPLINE - spline design matrix failed to invert") ;

        u = bandsolpd(qy,band(r,2)) ;

     endif ;

     @ Find derivatives, c0-c3, of smoothing spline at x @

     c2  = 0 | (p.*u) | 0 ; @ f''(x) @

     u   = ((u|0) - (0|u))./dx ;

     c0  = y - (t./w).*((u|0) - (0|u)) ; @ f(x) @

     c3  = (c2[2:n] - c2[1:n1])./dx | 0 ; @ f'''(x) @

     c1  = (c0[2:n] - c0[1:n1])./dx - dx.*(c2[1:n1]./3 + c2[2:n]./6) ;

     c1  = c1 | (c1[n1] + dx[n1].*(c2[n1] + c2[n])./2) ; @ f'(x) @

     retp(x~c0~c1~c2~c3) ;

endp ;



/* FSPLINE - Calculates fitted values for a cubic spline

**

** Usage:  y = fspline(c,x) ;

**

** Input:  c - N x 5 matrix of spline interval start points and coefficients

**         x - N x 1 vector of x values for the spline is to be evaluated at

**

** Output: y - N x 1 vector of fitted values corresponding to the x values

*/



proc fspline(c,x) ;

     local i,z ;

     @ I locates which interval each value of x is in @

     i = 1 + sumc(x' .> c[2:rows(c),1]) ;

     @ Calculate distance from start of interval that x falls in @

     z = x - c[i,1] ;

     @ Return cubic for x @

     retp(c[i,2] + z.*(c[i,3] + (z.>0).*z.*(c[i,4]./2 + z.*c[i,5]./6))) ;

endp ;



/* INV_SPL - Calculates inverse values for a monotonic cubic spline

**

** Usage:  x = inv_spl(c,y) ;

**

** Input:  c   - N x 5 matrix of spline interval start points and coefficients

**         y   - M x 1 vector of y values for the spline is to be evaluated at

**

** Output: x - M x 1 vector of x values such that y = f(x) where f is a

**             monotonically increasing or decreasing cubic spline function

**

** Globals: _spl_tol - controls the accuracy to which y is estimated

*/



proc inv_spl(c,y) ;

     local i,k,n,converge,xi,yi,ci,cd,xt,yt,xn ;

     n = rows(c) ;

     @ i locates which interval each value of y is in @

     if c[1,2] > c[2,2] ;

        i = 1 + sumc(y' .< c[2:n,2]) ;

     else ;

        i = 1 + sumc(y' .> c[2:n,2]) ;

     endif ;

     xi = c[i,1] ;

     yi = y - c[i,2] ;

     ci = c[i,3:5]./(1~2~6) ;

     cd = ci.*(1~2~3) ;

     clear xt,converge,k ;

     do until(converge or k > 20) ;

        yt = xt.*(ci[.,1] + (xt.>0).*xt.*(ci[.,2] + xt.*ci[.,3])) ;

        xn = xt + (yi - yt)./(cd[.,1] + (xt.>0).*xt.*(cd[.,2] + xt.*cd[.,3])) ;

        converge = abs(xt-xn) < _spl_tol and abs(yi-yt) < _spl_tol ;

        xt = xn ;

        k = k + 1 ;

     endo ;

     retp(xi + xt) ;

endp ;



/* DSPLINE - Calculates fitted values for a differentials of cubic spline

**

** Usage:  d = dspline(c,x,k) ;

**

** Input:  c - N x 5 matrix of spline interval start points and coefficients

**         x - N x 1 vector of x values for the spline is to be evaluated at

**         k - scalar giving the degree of the derivative

**

** Output: d - N x 1 vector of derivatives corresponding to the x values

*/



proc dspline(c,x,k) ;

     local i,z ;

     i = 1 + sumc(x' .> c[2:rows(c),1]) ;

     z = x - c[i,1] ;

     if    k == -1 ; @ Return integral from start of interval @

        retp(z.*(c[i,2]+z.*(c[i,3]./2+(z.>=0).*z.*(c[i,4]./6+z.*c[i,5]./24))));

     elseif k == 0 ; @ No derivative - just return cubic @

        retp(c[i,2] + z.*(c[i,3] + (z.>=0).*z.*(c[i,4]./2 + z.*c[i,5]./6))) ;

     elseif k == 1 ; @ First derivative - quadratic @

        retp(c[i,3] + (z.>=0).*z.*(c[i,4] + z.*c[i,5]./2)) ;

     elseif k == 2 ; @ Second derivative - linear @

        retp((z.>=0).*(c[i,4] + z.*c[i,5])) ;

     elseif k == 3 ; @ Third derivative - constant @

        retp((z.>=0).*c[i,5]) ;

     else ; @ All higher derivatives - zero @

        retp(zeros(rows(x),1)) ;

     endif ;

endp ;



/* ISPLINE - Calculates integral of cubic spline

**

** Usage:  area = ispline(c,x) ;

**

** Input:  c   - N x 5 matrix of spline interval start points and coefficients

**         x   - 2 x M vector of x limits for the integral

**               xl[1,.] = upper limits of x , xl[2,.] = lower limits of x

**

** Output: area - M x 1 matrix of integrals the cubic spline curve

*/



proc ispline(c,xl) ;

    local x,i,nc,l,u,in,area,fuzz ;

    if rows(xl) /= 2 ;

       errorlog "Error: ISPLINE - limits matrix has incorrect dimensions" ;

       retp(error(1)) ;

    endif ;

    x = c[.,1] ;

    nc = cols(xl) ;

    fuzz = meanc(abs(x)).*1e-12 ;

    area = zeros(nc,1) ;

    i = 1 ;

    do until(i > nc) ;

       u  = xl[1,i] ;

       l  = xl[2,i] ;

       in = selif(x,(x .> l) .and (x .< u)) ;

       if scalmiss(in) ;

          area[i] = dspline(c,u,-1) - dspline(c,l,-1) ;

       else ;

          area[i] = sumc(dspline(c,in-fuzz|u,-1)) - dspline(c,l+fuzz,-1) ;

       endif ;

       i = i + 1 ;

    endo ;

    retp(area) ;

endp ;



/* PRESS - Calculates PRESS statistic for spline fit using first principles.

**         The PRESS statistic is obtained by deleting each observation in turn

**         and obtaining the prediction residual for that point based on the

**         remaining observations. The PRESS statistic is the sum of squares

**         of the prediction residuals. Minimising the PRESS statistic for

**         the cubic spline over the range of smoothing parameters gives a

**         "best" choice of the smoothing parameter. This method of choosing

**         the smoothing parameter is known as Cross-Validation (CV).

**

** Usage:  z = press(x,y,w,p) ;

**

** Input:   x - N x 1 vector of x values of the points to be fitted

**              x must be sorted into ascending order

**          y - N x 1 vector of y values of the points to be fitted

**          w - N x 1 vector of weights for each point

**              For equal weighting set w to a scalar for efficiency

**          p - scalar giving the smoothing parameter

**

** Output: z - scalar giving the PRESS statistic.

**

** Note: The direct approach to calculating this statistic is very expensive in

**       computer time, but is more numerically stable as the p parameter tends

**       towards 0.

*/



proc press(x,y,w,p) ;

     local n,rdel,j,del,wd ;

     n = rows(x) ;

     rdel = zeros(n,1) ;

     j = 1 ;

     do until(j > n) ; @ Drop out one observation at a time @

        if j == 1 ;

           del = seqa(2,1,n-1) ;

        elseif j == n ;

           del = seqa(1,1,n-1) ;

        else ;

           del = seqa(1,1,j-1) | seqa(j+1,1,n-j) ;

        endif ;

        if rows(w) > 1 ;

           wd = w[del] ;

        else ;

           wd = w ;

        endif ;

        @ Obtain residual from prediction of deleted value @

        rdel[j] = y[j] - fspline(csspline(x[del],y[del],wd,p),x[j]) ;

        j = j + 1 ;

     endo ;

     @ Form PRESS statistic as mean of squared deletion residuals @

     retp(meanc(w.*rdel^2)) ;

endp ;



/* Internal proc to calculate the projection matrix A

** used in procs XVSC and GXVSC */



proc _spl_A(x,w,p) ; @ Creates globals _spl_q and _spl_r to speed iterations @

     local n,n1,n2,dx,dxi,t,qw,v,oldtrap ;

     w = w./meanc(w) ;

     n = rows(x) ;

     if rows(_spl_q /= n) ;

        n1 = n - 1 ;

        n2 = n - 2 ;

        dx = x[2:n] - x[1:n1] ;

        dxi = 1./dx ;

        @ Create R matrix which fits linear trend to y @

        _spl_r = rotater((0|dx[2:n2])~2.*(dx[1:n2]+dx[2:n1])~(dx[2:n2]|0)~

                     zeros(n2,n-5),seqa(-1,1,n2))./6 ;

        @ Create Q matrix which fits cubic spline to y @

        _spl_q = rotater(dxi[1:n2]~(-dxi[1:n2]-dxi[2:n1])~dxi[2:n1]

                 ~zeros(n2,n-3),seqa(0,1,n2))' ;

     endif ;

     t = (1 - p) ;

     qw = _spl_q./sqrt(w) ;

     qw = t.*moment(qw,0) + p.*_spl_r ;

     oldtrap = trapchk(0xffff) ;

     trap 1 ;

     v = invpd(qw) ;

     trap oldtrap,0xffff ;

     if scalmiss(v) ;

        errorlog("ERROR: _SPL_A - spline design matrix failed to invert") ;

        retp(eye(n)-x'invswp(moment(x./sqrt(w),0))*x) ;

     else ;

        retp(((t./w).*_spl_q)*v*_spl_q') ;

     endif ;

endp ;



/* XVSC - Calculates the Cross-Validation (CV) statistic

**

** Usage:  s = xvsc(x,y,w,p) ;

**

** Input:   x - N x 1 vector of x values of the points to be fitted

**              x must be sorted into ascending order

**          y - N x 1 vector of y values of the points to be fitted

**          w - N x 1 vector of weights for each point

**              For equal weighting set w to a scalar for efficiency

**          p - scalar giving the smoothing parameter

**

** Output:  s - scalar giving the CV statistic.

**

** Global:  _spl_A - proc to calculate the projection matrix A

**

** Note: The CV statistic is just the PRESS statistic but this procedure

**       provides more efficient estimation of it. Where p is close to 1

**       the PRESS statistic is calculated directly for numerical stability.

*/



proc xvsc(x,y,w,p) ;

     local A,res,cv ;

     if p < .9999 ;

        A = _spl_A(x,w,p) ;

        res = A*y ;

        cv = meanc(w.*(res./diag(A))^2) ;

     else ;

        cv = press(x,y,w,p) ;

     endif ;

     retp(cv) ;

endp ;



/* GXVSC - Calculates the Generalized Cross-Validation (GCV) statistic

**

** Usage:  s = gxvsc(x,y,w,p) ;

**

** Input:   x - N x 1 vector of x values of the points to be fitted

**              x must be sorted into ascending order

**          y - N x 1 vector of y values of the points to be fitted

**          w - N x 1 vector of weights for each point

**              For equal weighting set w to a scalar for efficiency

**          p - scalar giving the smoothing parameter

**

** Output:  s - scalar giving the GCV statistic.

**

** Global:  _spl_A - proc to calculate the projection matrix A

**

** Note: The GCV statistic reweights the contributions of the deletion

**       residuals to the CV statistic. Where p is close to 1 an

**       approximation to the GCV is calculated for numerical stability.

*/



proc gxvsc(x,y,w,p) ;

     local A,res,gcv ;

     if p < .9999 ;

        A = _spl_A(x,w,p) ;

        res = A*y ;

        gcv = meanc(w.*(res)^2)./(minc(meanc(diag(A))|rows(x))^2) ;

     else ;

        gcv = agxvsc(x,y,w,p) ;

     endif ;

     retp(gcv) ;

endp ;



/* Internal Procedures used in AGCV:

** _spl_fh, _spl_fh4, _spl_c_f, and _spl_trA */



proc _spl_fh(t) ; @ Kernel smooth of knots used by _spl_fh4 and spl_clim @

     local h,mx ;

     mx = meanc(_spl_w.*_spl_x) ;

     h = 1.06.*sqrt(meanc(_spl_w.*(_spl_x - mx)^2))./(rows(_spl_x)^(1/5)) ;

     if h == 0 ;

        errorlog("ERROR: _SPL_FH - _spl_w or _spl_x NOT set") ;

        retp(error(1)) ;

     endif ;

     retp(meanc(_spl_w.*pdfn((t'-_spl_x)./h))./h) ;

endp ;



proc _spl_fh4(t) ; @ Fourth root of _spl_fh integrated by proc _spl_c_f @

     retp(_spl_fh(t)^(1/4)) ;

endp ;



proc _spl_c_f(x,w) ; @ Evaluation of c0 constants for design points @ ;

     local n,xr,xl,xu,c0 ;

     n = rows(x) ;

     xr  = x[n] - x[1] ;

     xl  = x[1] - xr./(2.*n) ;

     xu  = x[n] + xr./(2.*n) ;

     _spl_x = x ; _spl_w = w ;

     c0 = intsimp(&_spl_fh4,xu|xl,_spl_tol) ;

     c0 = (((pi./c0).*seqa(1.5,1,n-2))^4)./n ;

     clearg _spl_x,_spl_w ;

     retp(c0) ;

endp ;



proc _spl_trA(x,w,p) ; @ Approximation of trace of A (the projection matrix) @

     if rows(_spl_c) /= rows(x) - 2 ;

        _spl_c = _spl_c_f(x,w./meanc(w)) ;

     endif ;

     retp(rows(x) - 2 - sumc(1./(1 + _spl_c.*((1-p)./p)))) ;

endp ;



/* AGXVSC - Calculates the Approximate Generalized Cross-Validation (AGCV)

**          statistic

**

** Usage:  s = agxvsc(x,y,w,p) ;

**

** Input:   x - N x 1 vector of x values of the points to be fitted

**              x must be sorted into ascending order

**          y - N x 1 vector of y values of the points to be fitted

**          w - N x 1 vector of weights for each point

**              For equal weighting set w to scalar 1 for efficiency

**          p - scalar giving the smoothing parameter

**

** Output:  s - scalar giving the AGCV.

**

** Globals: _spl_trA - proc - approximation of trace of projection matrix A

**          _spl_tol - scalar controls the accuracy of the approx.

**

** Note: The AGCV statistic is an approximation of the GCV statistic, but

**       provides a more efficient estimation procedure. Where p is very

**       close to 1 the PRESS statistic is calculated directly for numerical

**       stability. For p very close to 1, AGCV and PRESS converge.

*/



proc agxvsc(x,y,w,p) ;

     local n,trA,res,agcv ;

     n  =  rows(x) ;

     if p < .9999 ;

        trA = _spl_trA(x,w,p) ;

        res = y - fspline(csspline(x,y,w,p),x) ;

        agcv = meanc(w.*res^2)./(trA./n)^2 ;

     else ; @ Use limit of AGCV = PRESS for stability @

        agcv = press(x,y,w,p) ;

     endif ;

     retp(agcv) ;

endp ;



/* SPL_CLIM - Calculates the confidence limits of the mean for a spline.

**

** Usage:   l = spl_clim(x,y,w,p,pr) ;

**

** Input:   x  - N x 1 vector of x values of the points to be fitted

**               x must be sorted into ascending order

**          y  - N x 1 vector of y values of the points to be fitted

**          w  - N x 1 vector of weights for each point

**               For equal weighting set w to a scalar for efficiency

**          p  - scalar giving the smoothing parameter

**          pr - scalar giving probability of confidence limits (eg .95)

**               (0 < pr < 1)

**

** Output:  l  - N x 2 matrix giving the upper and lower pr confidence

**               limits of the mean.

**

** Globals: _spl_trA - proc - approximation of trace of projection matrix A

**          _spl_tol - scalar controls the accuracy of the approx.

**          _spl_A   - proc to calculate the projection matrix A

**

** Note: The GCV statistic reweights the contributions of the deletion

**       residuals to the CV statistic. Where p is close to 1 an

**       approximation to the GCV is calculated for numerical stability.

*/



proc spl_clim(x,y,w,p,pr) ;

     local k,n,A,dA,trA,res,se,l4,ui,fh ;

     if not (pr > 0 and pr < 1) ;

        errorlog "ERROR: SPL_CLIM - Invalid probability value in arguments" ;

        retp(error(1)) ;

     endif ;

     w = w./meanc(w) ;

     k = invcdfn((1 + pr)./2) ;

     n = rows(x) ;

     if p < .9 ;

        A = _spl_A(x,w,p) ;

        res = A*y ;

        dA  = diag(A) ;

        trA = sumc(dA) ;

        se  = sumc(w.*res^2)./(trA) ;

     elseif p <.9999 ;

        res = y - fspline(csspline(x,y,w,p),x) ;

        clearg _spl_c ;

        trA = _spl_trA(x,w,p) ;

        _spl_x = x ; _spl_w = w ;

        fh = _spl_fh(x) ;

        clearg _spl_x,_spl_w,_spl_c ;

        l4 = ((1 - p)./(n.*p))^(1/4) ;

        ui = minc(abs(x'- x[1 n])).*(sqrt(2)./l4).*fh^(1/4) ;

        dA = w.*(1+(2+cos(ui)-sin(ui)).*exp(-ui))./(l4.*n.*(4.*fh)^(3/4));

        dA = 1 - dA ;

        se = sumc(w.*res^2)./(trA) ;

     else ;

        se = press(x,y,w,p) ;

        dA = zeros(n,1) ;

     endif ;

     retp((y-res) + (-k'~k').*sqrt(se.*(1 - minc(dA'|ones(1,n)))./w)) ;

endp ;



/* Utility procedure used to minimise statistics using golden ratio search */

proc _spl_min(&stat,x,y,w) ;

    local stat:proc,r,c,x0,x1,x2,x3,f1,f2,try,maxtry,tol ;

    tol = 1e-4 ;

    maxtry = 30 ;

    r = (sqrt(5) - 1)/2;

    c = 1-r ;

    x0 = 0 ;

    x1 = r^2 ;

    x2 = r ;

    x3 = 1 ;

    f1 = stat(x,y,w,x1) ;

    f2 = stat(x,y,w,x2) ;

    try = 0 ;

    do until( abs(x3-x0) < tol or try > maxtry) ;

       if f2 < f1 ;

          x0 = x1 ;

          x1 = x2 ;

          x2 = r*x1 + c*x3 ;

          f1 = f2 ;

          f2 = stat(x,y,w,x2);

       else;

          x3 = x2 ;

          x2 = x1 ;

          x1 = r*x2 + c*x0 ;

          f2 = f1 ;

          f1 = stat(x,y,w,x1) ;

       endif;

       try = try + 1 ;

    endo ;

    if f1 < f2 ; retp(x1) ; else ; retp(x2) ; endif ;

endp ;



/* MIN_PRS - Calculate smoothing parameter which minimises the PRESS statistic.

**

** Usage:   p = min_prs(x,y,w) ;

**

** Input:   x - N x 1 vector of x values of the points to be fitted

**              x must be sorted into ascending order

**          y - N x 1 vector of y values of the points to be fitted

**          w - N x 1 vector of weights for each point to be fitted

**              For equal weighting set w to a scalar for efficiency

**

** Globals: _spl_min - golden ratio search proc to minimise PRESS

**

** Output:  p - scalar giving the smoothing parameter which minimizes

**              the PRESS statistic.

*/



proc min_prs(x,y,w) ;

     retp(_spl_min(&press,x,y,w)) ;

endp ;



/* MIN_CV - Calculate smoothing parameter which minimises the CV statistic.

**

** Usage:   p = min_cv(x,y,w) ;

**

** Input:   x - N x 1 vector of x values of the points to be fitted

**              x must be sorted into ascending order

**          y - N x 1 vector of y values of the points to be fitted

**          w - N x 1 vector of weights for each point to be fitted

**              For equal weighting set w to a scalar for efficiency

**

** Globals: _spl_min - golden ratio search proc to minimise CV

**          _spl_c   - approximation constant

**

** Output: p - scalar giving the smoothing parameter which minimizes

**             the CV statistic.

*/



proc min_cv(x,y,w) ;

     local p ;

     clearg _spl_c,_spl_r,_spl_q ;

     p = _spl_min(&xvsc,x,y,w) ;

     clearg _spl_c,_spl_r,_spl_q ;

     retp(p) ;

endp ;



/* MIN_GCV - Calculate smoothing parameter which minimises the GCV statistic.

**

** Usage:   p = min_gcv(x,y,w) ;

**

** Input:   x - N x 1 vector of x values of the points to be fitted

**              x must be sorted into ascending order

**          y - N x 1 vector of y values of the points to be fitted

**          d - N x 1 vector of weights for each point to be fitted

**              For equal weighting set w to a scalar for efficiency

**

** Globals: _spl_min - golden ratio search proc to minimise GCV

**          _spl_c   - approximation constant

**

** Output: p - scalar giving the smoothing parameter which minimizes

**             the GCV statistic.

*/



proc min_gcv(x,y,w) ;

     local p ;

     clearg _spl_c,_spl_r,_spl_q ;

     p = _spl_min(&gxvsc,x,y,w) ;

     clearg _spl_c,_spl_r,_spl_q ;

     retp(p) ;

endp ;



/* MIN_AGCV - Calculate smoothing parameter which minimises the AGCV statistic.

**

** Usage:   p = min_agcv(x,y,w) ;

**

** Input:   x - N x 1 vector of x values of the points to be fitted

**              x must be sorted into ascending order

**          y - N x 1 vector of y values of the points to be fitted

**          d - N x 1 vector of weights for each point to be fitted

**              For equal weighting set w to a scalar for efficiency

**

** Globals: _spl_min - golden ratio search proc to minimise AGCV

**          _spl_c   - approximation constant

**

** Output: p - scalar giving the smoothing parameter which minimizes

**             the AGCV statistic.

*/



proc min_agcv(x,y,w) ;

     local p ;

     clearg _spl_c ;

     p = _spl_min(&gxvsc,x,y,w) ;

     clearg _spl_c ;

     retp(p) ;

endp ;



/* PRSPLINE - Calculates cubic smoothing spline coefficients of the spline

**            which minimises the PRESS statistic.

**

** Usage:   c = prspline(x,y,w) ;

**

** Input:   x - N x 1 vector of x values of the points to be fitted

**              x must be sorted into ascending order

**          y - N x 1 vector of y values of the points to be fitted

**          d - N x 1 vector of weights for each point

**              For equal weighting set w to a scalar for efficiency

**

** Output:  c - N x 5 matrix of coefficents of smoothing spline

**              c[.,1] = knots,

**              c[.,2..5] = function value and first 3 derivatives at knots

*/



proc prspline(x,y,w) ;

     retp(csspline(x,y,w,min_prs(x,y,w))) ;

endp ;



/* CVSPLINE - Calculates cubic smoothing spline coefficients of the spline

**            which minimises the Cross-Validation statistic.

**

** Usage:   c = cvspline(x,y,w) ;

**

** Input:   x - N x 1 vector of x values of the points to be fitted

**              x must be sorted into ascending order

**          y - N x 1 vector of y values of the points to be fitted

**          w - N x 1 vector of weights for each point

**              For equal weighting set w to a scalar for efficiency

**

** Output:  c - N x 5 matrix of coefficents of smoothing spline

**              c[.,1] = knots,

**              c[.,2..5] = function value and first 3 derivatives at knots

**

** Note:     This uses CV rather than PRESS to calculate the CV statistic

**           and hence is computationally much faster than PRSPLINE

*/



proc cvspline(x,y,w) ;

     retp(csspline(x,y,w,min_cv(x,y,w))) ;

endp ;



/* GCSPLINE - Calculates cubic smoothing spline coefficients of the spline

**            which minimises the Generalized Cross-Validation statistic.

**

** Usage:   c = gcspline(x,y,w) ;

**

** Input:   x - N x 1 vector of x values of the points to be fitted

**              x must be sorted into ascending order

**          y - N x 1 vector of y values of the points to be fitted

**          w - N x 1 vector of weights for each point

**              For equal weighting set w to a scalar for efficiency

**

** Output:  c - N x 5 matrix of coefficents of smoothing spline

**              c[.,1] = knots,

**              c[.,2..5] = function value and first 3 derivatives at knots

*/



proc gcspline(x,y,w) ;

     retp(csspline(x,y,w,min_gcv(x,y,w))) ;

endp ;



/* AGSPLINE - Calculates cubic smoothing spline coefficients of the spline

**            which minimises the Approximate Generalized Cross-Validation

**            statistic.

**

** Usage:   c = agspline(x,y,w) ;

**

** Input:   x - N x 1 vector of x values of the points to be fitted

**              x must be sorted into ascending order

**          y - N x 1 vector of y values of the points to be fitted

**          w - N x 1 vector of weights for each point

**              For equal weighting set w to a scalar for efficiency

**

** Output:  c - N x 5 matrix of coefficents of smoothing spline

**              c[.,1] = knots,

**              c[.,2..5] = function value and first 3 derivatives at knots

*/



proc agspline(x,y,w) ;

     retp(csspline(x,y,w,min_agcv(x,y,w))) ;

endp ;



/* SPL_SET - Reset globals used in Spline Routines

**

** Usage:   call spl_set ;

**

** Globals: _spl_c,_spl_r,_spl_q,_spl_x,_spl_w,_spl_tol ;

*/



proc 0 = splset ;

     clearg _spl_c,_spl_r,_spl_q,_spl_x,_spl_w ;

     _spl_tol = 1e-6 ;

endp ;



/* INVCDFN - Calculates the inverse of the Normal cdf

**

** Usage:  z = invcdfn(p)

**

** Input:  p - N x M matrix of probabilities (0 < p < 1)

**

** Output: z - N x M matrix of normal scores s.t. prob(N < z) = p

**

** Note:   Used in SPL_CLIM so added to this archive for completeness

*/



proc invcdfn(p) ;

     local dz,tol,t,n,d,converge,k,z,f,df ;

     if not (p > 0 and p < 1) ;

        errorlog "ERROR: INVCDFN - Illegal probability value \G" ;

     endif ;

     dz = 1e-8 ;

     tol = 1e-6 ;

     /* Initial approximation of z */

     t = sqrt(-2*ln(abs((p.>0.5) - p)));

     n = 2.515517 + t.*(0.802853 + t.*0.010328) ;

     d = 1 + t.*(1.432788 + t.*(0.189269 + t.*0.001308)) ;

     z = t - (n./d) ;

     z = (p.>0.5).*z - (p.<=0.5).*z ;

     /* Newton-Raphson adjustment of approximation to give exact value */

     clear converge ;

     k = 1 ;

     do until(converge or k > 10) ;

         f = cdfn(z) ;

        df = (cdfn(z + dz) - f)./dz ;

         n = z - (f - p)./df ;

        converge = (abs(z - n) < tol) ;

        z = n ;

        k = k + 1 ;

     endo ;

     retp(z) ;

endp ;

