| C Output from Public domain Ratfor, version 1.0 |
| |
| c Smoothing Spline LeVeRaGes = SSLVRG |
| c ----------------------------------- leverages = H_ii = diagonal entries of Hat matrix |
| subroutine sslvrg(penalt,dofoff,x,y,w,ssw, n, knot,nk,coef, |
| * sz,lev, crit,icrit, lambda, xwy, hs0,hs1,hs2,hs3, |
| * sg0,sg1,sg2,sg3, abd,p1ip,p2ip,ld4,ldnk,info) |
| |
| C Purpose : |
| C Compute smoothing spline for smoothing parameter lambda |
| C and compute one of three `criteria' (OCV , GCV , "df match"). |
| C See comments in ./sbart.c from which this is called |
| |
| integer n,nk,icrit,ld4,ldnk,info |
| DOUBLE precision penalt,dofoff,x(n),y(n),w(n),ssw, |
| & knot(nk+4), coef(nk),sz(n),lev(n), crit, lambda, |
| * xwy(nk), hs0(nk),hs1(nk),hs2(nk),hs3(nk), |
| * sg0(nk),sg1(nk),sg2(nk),sg3(nk), abd(ld4,nk), |
| & p1ip(ld4,nk), p2ip(ldnk,nk) |
| |
| EXTERNAL bvalue |
| double precision bvalue |
| C local variables |
| double precision vnikx(4,1),work(16) |
| integer i,ileft,j,mflag, lenkno |
| double precision b0,b1,b2,b3,eps, xv,rss,df, sumw |
| c |
| integer interv |
| external interv ! in ../../../appl/interv.c |
| |
| lenkno = nk+4 |
| ileft = 1 |
| eps = 1d-11 |
| |
| C compute the coefficients coef() of estimated smooth |
| |
| do i=1,nk |
| coef(i) = xwy(i) |
| abd(4,i) = hs0(i)+lambda*sg0(i) |
| end do |
| |
| do i=1,(nk-1) |
| abd(3,i+1) = hs1(i)+lambda*sg1(i) |
| end do |
| |
| do i=1,(nk-2) |
| abd(2,i+2) = hs2(i)+lambda*sg2(i) |
| end do |
| |
| do i=1,(nk-3) |
| abd(1,i+3) = hs3(i)+lambda*sg3(i) |
| end do |
| |
| c factorize banded matrix abd (into upper triangular): |
| call dpbfa(abd,ld4,nk,3,info) |
| if(info.ne.0) then |
| C matrix could not be factorized -> ier := info |
| return |
| endif |
| c solve linear system (from factorized abd): |
| call dpbsl(abd,ld4,nk,3,coef) |
| |
| C Value of smooth at the data points |
| do i=1,n |
| xv = x(i) |
| sz(i) = bvalue(knot,coef,nk,4,xv,0) |
| end do |
| |
| C Compute the criterion function if requested (icrit > 0) : |
| if(icrit .ge. 1) then |
| |
| C --- Ordinary or Generalized CV or "df match" --- |
| |
| C Get Leverages First |
| call sinerp(abd,ld4,nk,p1ip,p2ip,ldnk, 0) |
| do i=1,n |
| xv = x(i) |
| ileft = interv(knot(1), nk+1, xv, 0,0, ileft, mflag) |
| if(mflag .eq. -1) then |
| ileft = 4 |
| xv = knot(4)+eps |
| else if(mflag .eq. 1) then |
| ileft = nk |
| xv = knot(nk+1) - eps |
| endif |
| j=ileft-3 |
| C call bspvd(knot,4,1,xv,ileft,4,vnikx,work) |
| call bsplvd(knot,lenkno,4,xv,ileft,work,vnikx,1) |
| b0=vnikx(1,1) |
| b1=vnikx(2,1) |
| b2=vnikx(3,1) |
| b3=vnikx(4,1) |
| lev(i) = ( |
| & p1ip(4,j)*b0**2 + 2.d0*p1ip(3,j)*b0*b1 + |
| * 2.d0*p1ip(2,j)*b0*b2 + 2.d0*p1ip(1,j)*b0*b3 + |
| * p1ip(4,j+1)*b1**2 + 2.d0*p1ip(3,j+1)*b1*b2 + |
| * 2.d0*p1ip(2,j+1)*b1*b3 + p1ip(4,j+2)*b2**2 + |
| & 2.d0*p1ip(3,j+2)*b2*b3 + p1ip(4,j+3)*b3**2 |
| & )*w(i)**2 |
| end do |
| |
| |
| C Evaluate Criterion |
| |
| df = 0d0 |
| if(icrit .eq. 1) then ! Generalized CV -------------------- |
| rss = ssw |
| sumw = 0d0 |
| c w(i) are sqrt( wt[i] ) weights scaled in ../R/smspline.R such |
| c that sumw = number of observations with w(i) > 0 |
| do i=1,n |
| rss = rss + ((y(i)-sz(i))*w(i))**2 |
| df = df + lev(i) |
| sumw = sumw + w(i)**2 |
| end do |
| |
| crit = (rss/sumw)/((1d0-(dofoff + penalt*df)/sumw)**2) |
| c call dblepr("spar", 4, spar, 1) |
| c call dblepr("crit", 4, crit, 1) |
| |
| else if(icrit .eq. 2) then ! Ordinary CV ------------------ |
| crit = 0d0 |
| do i = 1,n |
| crit = crit + (((y(i)-sz(i))*w(i))/(1-lev(i)))**2 |
| end do |
| crit = crit/n |
| c call dblepr("spar", 4, spar, 1) |
| c call dblepr("crit", 4, crit, 1) |
| |
| else ! df := sum( lev[i] ) |
| do i=1,n |
| df = df + lev(i) |
| end do |
| if(icrit .eq. 3) then ! df matching -------------------- |
| crit = 3 + (dofoff-df)**2 |
| else ! if(icrit .eq. 4) then df - dofoff (=> zero finding) |
| crit = df - dofoff |
| endif |
| endif |
| endif |
| C Criterion evaluation |
| return |
| end |
