https://github.com/cran/fields
Tip revision: 8f8fb01c96e0bd7cbf33a3c3066eb0cd7c8f1274 authored by Doug Nychka on 09 February 2006, 12:55:37 UTC
version 2.3
version 2.3
Tip revision: 8f8fb01
cvrcss.f
subroutine cvrcss(
c arguments passed through to rcss
+ n,x,y,wt,sy,diag,din,dout,
c
c arguments needed for cv search
+ nstep,maxit,trmin, trmax, hmin,hmax,hopt,vopt,
+ tropt,mxstep,tabout,
+ ierr)
implicit double precision (a-h,o-z)
real*8 tabout(mxstep,4)
real*8 x(n),y(n),wt(n),sy(n)
real*8 diag(n), din(10), dout(10)
real*8 hopt,vopt,tropt,trmin,trmax,hmin, hmax
data tau,tausq/.6180339d0,.3819660d0/
**** coarse search in bandwidth h: this search feeds (hl,hr) to golden search
c compute range for coarse search
if (mxstep .lt. nstep) then
ierr=10
return
endif
c********* begin block for finding hmin and hmax
if( hmin.ge.hmax) then
C*** for any spline problem it is assumed that lambda=1 is
c*** computable and can be used as a start value
hmax=0.0
do 1020 k=1,25
call cvrf(hmax,n,x,y,wt,sy,trace,diag,din,dout,cv,ierr)
if( trace.le.trmin) goto 1030
hmax= hmax + 1
1020 continue
ierr=31
1030 continue
hmin=0.0
do 1000 k=1,25
call cvrf(hmin,n,x,y,wt,sy,trace,diag,din,dout,cv,ierr)
if( trace.ge.trmax) goto 1010
hmin= hmin - 1
1000 continue
ierr=21
1010 continue
endif
c********* end block for finding hmin and hmax
hstep=(hmax-hmin)/(nstep-1)
do 3 j=0,nstep-1
h=hmax - j*hstep
call cvrf(h,n,x,y,wt,sy,trace,diag,din,dout,cvh,ierr)
tabout(j+1,1)=h
tabout(j+1,2)=dout(3)
tabout(j+1,3)=cvh
tabout(j+1,4)=dout(1)
if ((cvh .lt. cvmin) .or. (j .eq. 0)) then
hopt=h
best=h
cvmin=cvh
trbest=trace
endif
5000 format(5e12.4)
c
3 continue
c
c**** fast return if crude search minimum cv is hmin or hmax
hopt= best
vopt=cvmin
tropt=trbest
if( (best.le.hmin) .or. (best.ge.hmax) ) then
ierr=-1
return
endif
c**** start values for golden search
hl=best-hstep
hr=best+hstep
*** Golden section search for min cv on interval (hl,hr)--maxit iterations
*** cv(h) must be quasiconvex on initial interval (hl,hr).
*** On return, h=abscissa of minimum, v=cv(h)=actual minimum achieved.
*** Interval of uncertainty is (hl,hr), hl <= hlm <= hrm <= hr.
call
* cvrf(hl,n,x,y,wt,sy,trace,diag,din,dout,cvhl,ierr)
call
* cvrf(hr,n,x,y,wt,sy,trace,diag,din,dout,cvhr,ierr)
hlm=hl*tau+hr*tausq
hrm=hl+hr-hlm
call
* cvrf(hlm,n,x,y,wt,sy,trchlm,diag,din,dout,cvhlm,ierr)
call
* cvrf(hrm,n,x,y,wt,sy,trchrm,diag,din,dout,cvhrm,ierr)
do 5 it=1,maxit
if( cvhlm .ge. cvhrm ) then
if (cvhl .lt. cvhlm) then
err= cvhl/cvhlm
ierr=2
return
endif
c
hl=hlm
cvhl=cvhlm
hlm=hrm
hrm=hrm+(hrm-hl)*tau
cvhlm=cvhrm
call
* cvrf(hrm,n,x,y,wt,sy,trchrm,diag,din,dout,cvhrm,ierr)
else
if (cvhr .lt. cvhrm) then
err= cvhrm/cvhr
ierr=2
return
endif
c
hr=hrm
cvhr=cvhrm
hrm=hlm
hlm=hlm+(hlm-hr)*tau
cvhrm=cvhlm
call
* cvrf(hlm,n,x,y,wt,sy,trchlm,diag,din,dout,cvhlm,ierr)
endif
5 continue
**** finished -- take the best h so far
if( cvhlm .ge. cvhrm) then
hopt=hrm
vopt=cvhrm
tropt=trchrm
else
hopt=hlm
vopt=cvhlm
tropt=trchlm
endif
c
10 format(' cv(h) is NOT quasiconvex ',/,
- ' iter,error,hl,hlm: ',I4,e15.5,2e15.5)
11 format(' cv(h) is NOT quasiconvex ',/,
- ' iter,error,hr,hrm: ',I4,e15.5,2e15.5)
c
c
return
end