sslvrg.f

   1 C Output from Public domain Ratfor, version 1.0
   2       subroutine sslvrg(penalt,dofoff,x,y,w,ssw, n, knot,nk,coef,
   3      *     sz,lev, crit,icrit, lambda, xwy, hs0,hs1,hs2,hs3,
   4      *     sg0,sg1,sg2,sg3, abd,p1ip,p2ip,ld4,ldnk,info)
   5
   6 C Purpose :
   7 C       Compute smoothing spline for smoothing parameter lambda
   8 C       and compute one of three `criteria' (OCV , GCV , "df match").
   9 C See comments in ./sbart.f from which this is called
  10
  11       integer n,nk,icrit,ld4,ldnk,info
  12       DOUBLE precision penalt,dofoff,x(n),y(n),w(n),ssw,
  13      &     knot(nk+4), coef(nk),sz(n),lev(n), crit, lambda,
  14      *     xwy(nk), hs0(nk),hs1(nk),hs2(nk),hs3(nk),
  15      *     sg0(nk),sg1(nk),sg2(nk),sg3(nk), abd(ld4,nk),
  16      &     p1ip(ld4,nk),p2ip(ldnk,nk)
  17
  18       EXTERNAL bvalue
  19       double precision bvalue
  20 C local variables
  21       double precision vnikx(4,1),work(16)
  22       integer i,ileft,j,mflag, lenkno, ilo
  23       double precision b0,b1,b2,b3,eps, xv,rss,df, sumw
  24 c
  25 c      integer interv
  26 c      external interv
  27
  28       lenkno = nk+4
  29       ileft = 1
  30       eps = 1d-11
  31           ilo = 1
  32
  33 C compute the coefficients coef() of estimated smooth
  34
  35       do 1 i=1,nk
  36          coef(i) = xwy(i)
  37          abd(4,i) = hs0(i)+lambda*sg0(i)
  38  1    continue
  39
  40       do 4 i=1,(nk-1)
  41          abd(3,i+1) = hs1(i)+lambda*sg1(i)
  42  4    continue
  43
  44       do 6 i=1,(nk-2)
  45  6       abd(2,i+2) = hs2(i)+lambda*sg2(i)
  46
  47       do 8 i=1,(nk-3)
  48  8       abd(1,i+3) = hs3(i)+lambda*sg3(i)
  49
  50 c     factorize banded matrix abd:
  51       call dpbfa(abd,ld4,nk,3,info)
  52       if(info.ne.0) then
  53 C        matrix could not be factorized -> ier := info
  54          return
  55       endif
  56 c     solve linear system (from factorize abd):
  57       call dpbsl(abd,ld4,nk,3,coef)
  58
  59 C     Value of smooth at the data points
  60       do 12 i=1,n
  61          xv = x(i)
  62  12      sz(i) = bvalue(knot,coef,nk,4,xv,0)
  63
  64 C     Compute the criterion function if requested
  65
  66       if(icrit .eq. 0)then
  67          return
  68       else
  69 C --- Ordinary or Generalized CV or "df match" ---
  70
  71 C     Get Leverages First
  72          call sinerp(abd,ld4,nk,p1ip,p2ip,ldnk,0)
  73          do 16 i=1,n
  74             xv = x(i)
  75             call interv(knot(1), nk+1, xv, ileft, mflag)
  76             if(mflag .eq. -1) then
  77                ileft = 4
  78                xv = knot(4)+eps
  79             else if(mflag .eq. 1) then
  80                ileft = nk
  81                xv = knot(nk+1) - eps
  82             endif
  83             j=ileft-3
  84 C           call bspvd(knot,4,1,xv,ileft,4,vnikx,work)
  85             call bsplvd(knot,lenkno,4,xv,ileft,work,vnikx,1)
  86             b0=vnikx(1,1)
  87             b1=vnikx(2,1)
  88             b2=vnikx(3,1)
  89             b3=vnikx(4,1)
  90             lev(i) = (
  91      &              p1ip(4,j)*b0**2   + 2.d0*p1ip(3,j)*b0*b1 +
  92      *           2.d0*p1ip(2,j)*b0*b2   + 2.d0*p1ip(1,j)*b0*b3 +
  93      *              p1ip(4,j+1)*b1**2 + 2.d0*p1ip(3,j+1)*b1*b2 +
  94      *           2.d0*p1ip(2,j+1)*b1*b3 +    p1ip(4,j+2)*b2**2 +
  95      &           2.d0*p1ip(3,j+2)*b2*b3 +    p1ip(4,j+3)*b3**2
  96      &           )*w(i)**2
  97  16      continue
  98
  99 C     Evaluate Criterion
 100
 101          if(icrit .eq. 1)then
 102 C     Generalized CV
 103             rss = ssw
 104             df = 0d0
 105             sumw = 0d0
 106 c       w(i) are sqrt( wt[i] ) weights scaled in ../R/smspline.R such
 107 c       that sumw =  number of observations with w(i) > 0
 108             do 24 i=1,n
 109                rss = rss + ((y(i)-sz(i))*w(i))**2
 110                df = df + lev(i)
 111                sumw = sumw + w(i)**2
 112  24         continue
 113
 114             crit = (rss/sumw)/((1d0-(dofoff + penalt*df)/sumw)**2)
 115 c            call dblepr("spar", 4, spar, 1)
 116 c            call dblepr("crit", 4, crit, 1)
 117
 118          else if(icrit .eq. 2) then
 119 C     Ordinary CV
 120                crit = 0d0
 121                do 30 i = 1,n
 122  30               crit = crit + (((y(i)-sz(i))*w(i))/(1-lev(i)))**2
 123                crit = crit/n
 124 c            call dblepr("spar", 4, spar, 1)
 125 c            call dblepr("crit", 4, crit, 1)
 126             else
 127 C     df matching
 128             crit = 0d0
 129             do 32 i=1,n
 130  32            crit = crit+lev(i)
 131             crit = 3 + (dofoff-crit)**2
 132          endif
 133          return
 134       endif
 135 C     Criterion evaluation
 136       end