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