978 lines
31 KiB
FortranFixed
978 lines
31 KiB
FortranFixed
! Multivariate Universal Maximum Likelihood Estimation (MUMLE)
|
|
! Developed by
|
|
! Lianhong Gu
|
|
! Environmental Sciences Division
|
|
! Oak Ridge National Lab
|
|
! lianhong-gu@ornl.gov
|
|
!
|
|
! Version 1.0, July, 2006
|
|
|
|
subroutine coapproxyfunc(nkrigpoints,nxdim,xsampkrig,
|
|
& ysampkrig,nvarscp,indexvarscp,iupdatelikelihood,
|
|
& iwhattodo,ithetap0,ixnewpos,xnew,
|
|
& ypred,grad_xnew,rootmeanerror)
|
|
implicit none
|
|
include '../maxlikelihood/colikelihood.h'
|
|
|
|
!!!!!!!!!!!!!!!!!!!! Arguments !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
|
!-------------------Inputs--------------------------------------
|
|
!nkrigpoints: the number of available points for kriging.
|
|
!nxdim: the dimension in the design parameter vector x
|
|
!xsampkrig: the coordinates of the sampled points for kriging. Note that these
|
|
! points may or may not be the same as those points used for estimating
|
|
! the parameters in the likelihood function
|
|
!
|
|
! ****************************************************************************
|
|
! *Note different variables are stacked in the vector (xsampkrig,ysampkrig) *
|
|
! *in an orderly fashion. The number of different variables is given by *
|
|
! *nvarscp. The number of points in each variable is indicated by indexvarscp*
|
|
! ****************************************************************************
|
|
!
|
|
!ysampkrig: the value at xsampkrig
|
|
!nvarscp: the number of variables for cokriging
|
|
!indexvarscp: the number of samples in each variable for cokriging
|
|
!iupdatelikelihood: whether to update the likelihood function
|
|
! =0: the first entering, so initilize the parameters to be optimized
|
|
! and do likelihood function optimization
|
|
! =1: update the likelihood parameters using the saved values from
|
|
! the last optimization as the initial guesses
|
|
! =2: do not do likelihood optimization
|
|
!iwhattodo: =1, only make a predition at xnew
|
|
! =2, make a prediction and calculate the derivative at xnew
|
|
! =3, only calculate the derivative at xnew
|
|
!ithetap0: =ithetap: =1 set pdist to 2.0, =2: pdist to be optimized
|
|
!ixnewpos: which variable does xnew belong to.
|
|
!xnew: the coordinates of the point to be kriged
|
|
!
|
|
!---------------------- Outputs --------------------------------
|
|
!ypred: the predicted value at xnew
|
|
!rootmeanerror: the root mean square error of ypred
|
|
!grad_xnew: the derivative at xnew
|
|
|
|
integer nkrigpoints,nxdim,iupdatelikelihood,iwhattodo,
|
|
& nvarscp,indexvarscp(nvarscp),ithetap0,ixnewpos
|
|
double precision xsampkrig(nkrigpoints,nxdim),
|
|
& ysampkrig(nkrigpoints),xnew(nxdim),ypred,
|
|
& rootmeanerror,grad_xnew(nxdim)
|
|
!
|
|
!!!!!!!!!!!!!!!!!!! Local variables !!!!!!!!!!!!!!!!!!!!!!!!!
|
|
!
|
|
double precision ftol
|
|
parameter(ftol=1.0d-8)
|
|
external df1dimd_like,gradlikefunc,
|
|
& likelihoodfunc,f1dim_like,colikelihoodgrad
|
|
integer i,j,nplike,info,nii,index,irowinc,istart,
|
|
& icompactpostri,icompactposoff
|
|
integer nbd(maxnvars*(maxnvars+1)*
|
|
& (2*maxdesignparam+1)/2)
|
|
double precision plike(maxnvars*(maxnvars+1)*
|
|
& (2*maxdesignparam+1)/2),
|
|
& x(nxdim),flikemin,xsamp(nxdim),
|
|
& theta((maxnvars*(maxnvars+1))/2,maxdesignparam),
|
|
& pdist((maxnvars*(maxnvars+1))/2,maxdesignparam),
|
|
& theta1(nxdim),pdist1(nxdim),
|
|
& covars(maxnvars),correls(maxnvars*(maxnvars-1)/2),
|
|
& rx(nkrigpoints),fb(maxkbasvar),corr_like,
|
|
& xmean(maxdesignparam),
|
|
& gradient(maxkbasvar,nxdim)
|
|
|
|
double precision gradlike(maxnvars*(maxnvars+1)*
|
|
& (2*maxdesignparam+1)/2)
|
|
save nplike
|
|
save plike,theta,pdist,xmean,covars,correls
|
|
|
|
if(iupdatelikelihood.le.1)then
|
|
! update parameters and correlation functions in the likelihood function
|
|
!
|
|
! nbd is an INTEGER array of dimension n that must be set by the
|
|
! user to the type of bounds imposed on the variables:
|
|
! nbd(i)=0 if x(i) is unbounded,
|
|
! 1 if x(i) has only a lower bound,
|
|
! 2 if x(i) has both lower and upper bounds,
|
|
! 3 if x(i) has only an upper bound.
|
|
|
|
nlikepoints=nkrigpoints
|
|
nvars=nvarscp
|
|
ithetap=ithetap0
|
|
do i=1,nvars
|
|
indexvars(i)=indexvarscp(i)
|
|
enddo
|
|
ndesignparam=nxdim
|
|
do i=1,nlikepoints
|
|
ysamplike(i)=ysampkrig(i)
|
|
enddo
|
|
|
|
! Centralize the coordinates to stablize the linear system
|
|
do j=1,nxdim
|
|
xmean(j)=0.0d0
|
|
do i=1,nlikepoints
|
|
xmean(j)=xmean(j)+xsampkrig(i,j)
|
|
enddo
|
|
xmean(j)=xmean(j)/dble(nlikepoints)
|
|
do i=1,nlikepoints
|
|
xsamplike(i,j)=xsampkrig(i,j)-xmean(j)
|
|
enddo
|
|
enddo
|
|
|
|
call setinit_bounds(nvars,ndesignparam,ithetap,
|
|
& iupdatelikelihood,plikemin,plike,plikemax,nbd,nplike)
|
|
|
|
isitlastcall=.false.
|
|
|
|
! to use Lbfgsb_2_4, no need to call the cost function first
|
|
call Lbfgsb_2_4(nplike,plike,flikemin,plikemin,plikemax,
|
|
& nbd,colikelihoodgrad,ftol,info)
|
|
|
|
! to use frprmn, the cost function needs to be called first
|
|
! call likelihoodfunc(nplike,plike,flikemin)
|
|
! call frprmn(plike,nplike,ftol,flikemin,plikemin,
|
|
! & plikemax,gradlikefunc,f1dim_like,df1dimd_like)
|
|
|
|
! call dfpmin(plike,nplike,ftol,flikemin,likelihoodfunc,
|
|
! & plikemin,plikemax,gradlikefunc,f1dim_like,df1dimd_like)
|
|
|
|
write(*,*)flikemin,plike(1),plike(2),plike(3),info
|
|
call gradlikefunc(nplike,plike,flikemin,gradlike)
|
|
write(*,*)gradlike(1),gradlike(2),gradlike(3)
|
|
|
|
call nongradopt(nplike,likelihoodfunc,f1dim_like,plike,
|
|
& plikemin,plikemax,ftol,flikemin)
|
|
|
|
write(*,*)flikemin,plike(1),plike(2),plike(3)
|
|
call gradlikefunc(nplike,plike,flikemin,gradlike)
|
|
write(*,*)gradlike(1),gradlike(2),gradlike(3)
|
|
|
|
isitlastcall=.true.
|
|
call likelihoodfunc(nplike,plike,flikemin)
|
|
call colikeliparampart(nplike,plike,ithetap,
|
|
& ndesignparam,nvars,
|
|
& theta(1:(nvars*(nvars+1))/2,1:ndesignparam),
|
|
& pdist(1:(nvars*(nvars+1))/2,1:ndesignparam),
|
|
& covars,correls)
|
|
|
|
endif
|
|
!
|
|
do j=1,ndesignparam
|
|
x(j)=xnew(j)-xmean(j)
|
|
enddo
|
|
istart=0
|
|
do i=1,ixnewpos-1
|
|
istart=istart+kbasvar(i)
|
|
enddo
|
|
!
|
|
if(iwhattodo.eq.1.or.iwhattodo.eq.2)then
|
|
! make a prediction at x
|
|
call cofuncbasis(ixnewpos,index,ndesignparam,
|
|
& x,fb)
|
|
irowinc=0
|
|
do i=1,nvars
|
|
index=icompactpostri(ixnewpos,i)
|
|
do j=1,ndesignparam
|
|
theta1(j)=theta(index,j)
|
|
pdist1(j)=pdist(index,j)
|
|
enddo
|
|
if(ixnewpos.ne.i)then
|
|
index=icompactposoff(ixnewpos,i)
|
|
endif
|
|
do nii=1,indexvars(i)
|
|
irowinc=irowinc+1
|
|
do j=1,ndesignparam
|
|
xsamp(j)=xsamplike(irowinc,j)
|
|
enddo
|
|
if(ixnewpos.eq.i)then
|
|
rx(irowinc)=covars(i)*corr_like(x,xsamp,
|
|
& theta1,pdist1,ndesignparam)
|
|
else
|
|
rx(irowinc)=dsqrt(covars(i)*covars(ixnewpos))*
|
|
& correls(index)*corr_like(x,xsamp,
|
|
& theta1,pdist1,ndesignparam)
|
|
endif
|
|
enddo
|
|
enddo
|
|
|
|
ypred=0.0d0
|
|
do i=1,kbasvar(ixnewpos)
|
|
ypred=ypred+fb(i)*coeffbas(i+istart)
|
|
enddo
|
|
do i=1,nlikepoints
|
|
ypred=ypred+delvector(i)*rx(i)
|
|
enddo
|
|
|
|
rootmeanerror=covars(ixnewpos)
|
|
do i=1,kbasis
|
|
call symmatindex(kbasis,i,kpvt)
|
|
flikemin=0.0d0
|
|
do j=1,kbasis
|
|
if(j.gt.istart.and.j.le.(istart+kbasvar(ixnewpos)))then
|
|
flikemin=flikemin+FtRinvF_inv(kpvt(j))*fb(j-istart)
|
|
endif
|
|
enddo
|
|
if(i.gt.istart.and.i.le.(istart+kbasvar(ixnewpos)))then
|
|
rootmeanerror=rootmeanerror+fb(i-istart)*flikemin
|
|
endif
|
|
enddo
|
|
|
|
do i=1,kbasis
|
|
flikemin=0.0d0
|
|
do j=1,nlikepoints
|
|
flikemin=flikemin+fbassamp(j,i)*rx(j)
|
|
enddo
|
|
if(i.gt.istart.and.i.le.(istart+kbasvar(ixnewpos)))then
|
|
rootmeanerror=rootmeanerror-2.0d0*fb(i)*flikemin
|
|
endif
|
|
enddo
|
|
|
|
do i=1,nlikepoints
|
|
call symmatindex(nlikepoints,i,kpvt)
|
|
flikemin=0.0d0
|
|
do j=1,nlikepoints
|
|
flikemin=flikemin+work(kpvt(j))*rx(j)
|
|
enddo
|
|
rootmeanerror=rootmeanerror+rx(i)*flikemin
|
|
enddo
|
|
if(rootmeanerror.lt.0.0d0)then
|
|
! computational error
|
|
rootmeanerror=0.0d0
|
|
endif
|
|
rootmeanerror=dsqrt(rootmeanerror)
|
|
endif
|
|
!
|
|
if(iwhattodo.eq.2.or.iwhattodo.eq.3)then
|
|
! calculate the derivative of the approximated function with respect to x
|
|
call cogradbasx(ixnewpos,kbasvar(ixnewpos),ndesignparam,
|
|
& x,gradient(1:kbasvar(ixnewpos),1:ndesignparam))
|
|
do i=1,ndesignparam
|
|
grad_xnew(i)=0.0d0
|
|
do j=1,kbasvar(ixnewpos)
|
|
grad_xnew(i)=grad_xnew(i)+gradient(j,i)*coeffbas(j+istart)
|
|
enddo
|
|
enddo
|
|
|
|
irowinc=0
|
|
do i=1,nvars
|
|
if(i.eq.ixnewpos)then
|
|
flikemin=covars(ixnewpos)
|
|
else
|
|
index=icompactposoff(i,ixnewpos)
|
|
flikemin=dsqrt(covars(i)*covars(ixnewpos))*
|
|
& correls(index)
|
|
endif
|
|
index=icompactpostri(i,ixnewpos)
|
|
do j=1,ndesignparam
|
|
theta1(j)=theta(index,j)
|
|
pdist1(j)=pdist(index,j)
|
|
enddo
|
|
do j=1,indexvars(i)
|
|
irowinc=irowinc+1
|
|
do index=1,ndesignparam
|
|
xsamp(index)=xsamplike(irowinc,index)
|
|
enddo
|
|
call gradcorratx(x,xsamp,ndesignparam,theta1,
|
|
& pdist1,rx)
|
|
do index=1,ndesignparam
|
|
grad_xnew(index)=grad_xnew(index)+flikemin*
|
|
& rx(index)*delvector(irowinc)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
endif
|
|
return
|
|
end
|
|
|
|
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
|
|
subroutine colikelihoodgrad(nplike,plike,dogradient,
|
|
& fvalue,gradlike,ierr)
|
|
implicit none
|
|
include '../maxlikelihood/colikelihood.h'
|
|
!
|
|
! calculate the transformed likelihood function value and derivatives
|
|
! with respect to input parameters if specified
|
|
!
|
|
!----------- Arguments ------------------------------------
|
|
integer nplike,ierr
|
|
logical dogradient
|
|
double precision plike(nplike),fvalue,gradlike(nplike)
|
|
!
|
|
!=> nplike: the number of likelihood parameters to be estimated
|
|
!=> plike: the likelihood parameters to be estimated
|
|
!=> dogradient: whether to do derivative calculations
|
|
!<= fvalue: the transformed likelihood function value
|
|
!<= gradlike: derivatives at plike
|
|
!<= ierr: error messages
|
|
! =0: every thing ok
|
|
! <0: parameter out of bounds. |ierr| denotes the out-of-bound parameter
|
|
! =11: the correlation matrix is not postitive definite (the determinant is negative)
|
|
|
|
!----------- End of list of arguments ---------------------
|
|
!------------- Local variables ----------------------------
|
|
integer i,j,n,t,inert(3),index,nii,
|
|
& irowinc,getwhichvar,ivar,jvar,icount,k,
|
|
& icompactposoff,icompactpostri,iloop,istart,jstart
|
|
double precision lnabsdet,signunit,delta,corr_like,
|
|
& covars(nvars),theta(nvars*(nvars+1)/2,ndesignparam),
|
|
& pdist(nvars*(nvars+1)/2,ndesignparam),
|
|
& correls(nvars*(nvars-1)/2),xi(ndesignparam),
|
|
& xj(ndesignparam),theta1(ndesignparam),
|
|
& pdist1(ndesignparam),derR(nlikepoints*(nlikepoints+1)/2),
|
|
& sumi,sumt,gradtheta(ndesignparam),gradp(ndesignparam)
|
|
|
|
do i=1,nplike
|
|
if(plike(i).lt.plikemin(i).or.
|
|
& plike(i).gt.plikemax(i))then
|
|
! penalize out-of-bounds parameters
|
|
fvalue=1.0d+100
|
|
ierr=-i
|
|
return
|
|
endif
|
|
enddo
|
|
|
|
call colikeliparampart(nplike,plike,ithetap,
|
|
& ndesignparam,nvars,
|
|
& theta(1:(nvars*(nvars+1))/2,1:ndesignparam),
|
|
& pdist(1:(nvars*(nvars+1))/2,1:ndesignparam),
|
|
& covars,correls)
|
|
|
|
t=0
|
|
do j=1,nlikepoints
|
|
jvar=getwhichvar(j,nvars,indexvars)
|
|
do n=1,ndesignparam
|
|
xi(n)=xsamplike(j,n)
|
|
enddo
|
|
do i=1,j
|
|
t=t+1
|
|
if(i.eq.j)then
|
|
work(t)=covars(jvar)
|
|
else
|
|
ivar=getwhichvar(i,nvars,indexvars)
|
|
k=icompactpostri(ivar,jvar)
|
|
do n=1,ndesignparam
|
|
xj(n)=xsamplike(i,n)
|
|
theta1(n)=theta(k,n)
|
|
pdist1(n)=pdist(k,n)
|
|
enddo
|
|
if(ivar.eq.jvar)then
|
|
work(t)=covars(jvar)*corr_like(xj,xi,
|
|
& theta1,pdist1,ndesignparam)
|
|
else
|
|
k=icompactposoff(ivar,jvar)
|
|
work(t)=dsqrt(covars(jvar)*covars(ivar))*
|
|
& correls(k)*corr_like(xj,xi,
|
|
& theta1,pdist1,ndesignparam)
|
|
endif
|
|
endif
|
|
enddo
|
|
enddo
|
|
|
|
call dspfa(work,nlikepoints,kpvt,i)
|
|
j=11
|
|
call dspdi(work,nlikepoints,kpvt,derR,inert,delvector,j)
|
|
! work now stores the upper triangle of the inverse covariance matrix
|
|
! determinant=derR(1)*10.0d0**derR(2)
|
|
signunit=dsign(1.0d0,derR(1))
|
|
lnabsdet=dlog(dabs(derR(1)))+derR(2)*dlog(10.0d0)
|
|
|
|
if(signunit.lt.0.0d0)then
|
|
! write(*,*)'Wrong! the covariance matrix has a negative'
|
|
! write(*,*)'determinant. Program stops in likelihoodfunc'
|
|
write(*,*)signunit,lnabsdet
|
|
ierr=11
|
|
return
|
|
endif
|
|
|
|
t=0
|
|
kbasis=0
|
|
do i=1,nvars
|
|
do j=1,indexvars(i)
|
|
t=t+1
|
|
do n=1,ndesignparam
|
|
xi(n)=xsamplike(t,n)
|
|
enddo
|
|
call cofuncbasis(i,kbasvar(i),ndesignparam,
|
|
& xi,delvector)
|
|
do n=1,kbasvar(i)
|
|
fbassamp(t,kbasis+n)=delvector(n)
|
|
enddo
|
|
enddo
|
|
if(i.gt.1)then
|
|
do j=t-indexvars(i)+1,t
|
|
do n=1,kbasis
|
|
fbassamp(j,n)=0.0d0
|
|
enddo
|
|
enddo
|
|
do j=1,t
|
|
do n=kbasis+1,kbasis+kbasvar(i)
|
|
fbassamp(j,n)=0.0d0
|
|
enddo
|
|
enddo
|
|
endif
|
|
kbasis=kbasis+kbasvar(i)
|
|
enddo
|
|
call solvebassystem(delta)
|
|
|
|
fvalue=delta+lnabsdet
|
|
!
|
|
! The actual likelihood function =exp(-(fvalue+nlikepoints*ln(2*pi))/2)
|
|
|
|
if(dogradient.eqv..false..or.dogradient.eqv..FALSE.)then
|
|
ierr=0
|
|
return
|
|
endif
|
|
if(isitlastcall.eqv..true..or.isitlastcall.eqv..TRUE.)then
|
|
ierr=0
|
|
return
|
|
endif
|
|
!
|
|
! section for derivatives starts
|
|
!
|
|
! delvector=R^(-1)(Y-FB) was computed in solvebassystem and will not be altered from now
|
|
!----Section for calculating derivatives with respect to theta and pdist-------
|
|
icount=0
|
|
do iloop=1,ithetap
|
|
do k=1,nvars*(nvars+1)/2
|
|
call ivarjvartri(k,nvars,indexvars,
|
|
& ivar,jvar,istart,jstart)
|
|
if(ivar.ne.jvar)then
|
|
nii=icompactposoff(ivar,jvar)
|
|
endif
|
|
do n=1,ndesignparam
|
|
theta1(n)=theta(k,n)
|
|
pdist1(n)=pdist(k,n)
|
|
enddo
|
|
ierr=0
|
|
do j=1,nlikepoints
|
|
do i=1,j
|
|
ierr=ierr+1
|
|
derR(ierr)=0.0d0
|
|
enddo
|
|
enddo
|
|
do t=1,ndesignparam
|
|
icount=icount+1
|
|
do j=jstart,jstart+indexvars(jvar)-1
|
|
do n=1,ndesignparam
|
|
xj(n)=xsamplike(j,n)
|
|
enddo
|
|
if(ivar.eq.jvar)then
|
|
kpvt(1)=j-1
|
|
else
|
|
kpvt(1)=istart+indexvars(ivar)-1
|
|
endif
|
|
do i=istart,kpvt(1)
|
|
do n=1,ndesignparam
|
|
xi(n)=xsamplike(i,n)
|
|
enddo
|
|
ierr=icompactpostri(i,j)
|
|
if(ivar.eq.jvar)then
|
|
if(iloop.eq.1)then
|
|
call gradcorr_liketheta(xi,xj,
|
|
& theta1,pdist1,ndesignparam,gradtheta)
|
|
derR(ierr)=covars(ivar)*gradtheta(t)
|
|
else
|
|
call gradcorr_likep(xi,xj,
|
|
& theta1,pdist1,ndesignparam,gradp)
|
|
derR(ierr)=covars(ivar)*gradp(t)
|
|
endif
|
|
else
|
|
if(iloop.eq.1)then
|
|
call gradcorr_liketheta(xi,xj,
|
|
& theta1,pdist1,ndesignparam,gradtheta)
|
|
derR(ierr)=dsqrt(covars(ivar)*covars(jvar))*
|
|
& correls(nii)*gradtheta(t)
|
|
else
|
|
call gradcorr_likep(xi,xj,
|
|
& theta1,pdist1,ndesignparam,gradp)
|
|
derR(ierr)=dsqrt(covars(ivar)*covars(jvar))*
|
|
& correls(nii)*gradp(t)
|
|
endif
|
|
endif
|
|
enddo
|
|
enddo
|
|
gradlike(icount)=0.0d0
|
|
do i=1,nlikepoints
|
|
call symmatindex(nlikepoints,i,kpvt)
|
|
sumi=0.0d0
|
|
do j=1,nlikepoints
|
|
sumi=sumi+derR(kpvt(j))*delvector(j)
|
|
enddo
|
|
gradlike(icount)=
|
|
& gradlike(icount)+delvector(i)*sumi
|
|
enddo
|
|
gradlike(icount)=-gradlike(icount)
|
|
!
|
|
!Now, the trace part
|
|
sumi=0.0d0
|
|
do i=1,nlikepoints
|
|
call symmatindex(nlikepoints,i,kpvt)
|
|
sumt=0.0d0
|
|
do j=1,nlikepoints
|
|
sumt=sumt+work(kpvt(j))*derR(kpvt(j))
|
|
enddo
|
|
sumi=sumi+sumt
|
|
enddo
|
|
! sumi is the trace
|
|
gradlike(icount)=gradlike(icount)+sumi
|
|
enddo
|
|
enddo
|
|
enddo
|
|
!-------------------------------------------------------------------
|
|
! ------- Section for calculating derivatives with respect to variance
|
|
do k=1,nvars
|
|
icount=icount+1
|
|
ierr=0
|
|
do j=1,nlikepoints
|
|
do i=1,j
|
|
ierr=ierr+1
|
|
derR(ierr)=0.0d0
|
|
enddo
|
|
enddo
|
|
! first, vertical elements
|
|
do i=1,k-1
|
|
ierr=icompactpostri(i,k)
|
|
do t=1,ndesignparam
|
|
theta1(t)=theta(ierr,t)
|
|
pdist1(t)=pdist(ierr,t)
|
|
enddo
|
|
ierr=icompactposoff(i,k)
|
|
call istartjstart(i,k,nvars,indexvars,istart,jstart)
|
|
do n=istart,istart+indexvars(i)-1
|
|
do t=1,ndesignparam
|
|
xi(t)=xsamplike(n,t)
|
|
enddo
|
|
do j=jstart,jstart+indexvars(k)-1
|
|
do t=1,ndesignparam
|
|
xj(t)=xsamplike(j,t)
|
|
enddo
|
|
ivar=icompactpostri(n,j)
|
|
derR(ivar)=0.5d0*dsqrt(covars(i)/covars(k))*
|
|
& correls(ierr)*corr_like(xi,xj,
|
|
& theta1,pdist1,ndesignparam)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
ierr=icompactpostri(k,k)
|
|
do t=1,ndesignparam
|
|
theta1(t)=theta(ierr,t)
|
|
pdist1(t)=pdist(ierr,t)
|
|
enddo
|
|
call istartjstart(k,k,nvars,indexvars,istart,jstart)
|
|
do j=jstart,jstart+indexvars(k)-1
|
|
do t=1,ndesignparam
|
|
xi(t)=xsamplike(j,t)
|
|
enddo
|
|
do n=istart,j
|
|
ivar=icompactpostri(n,j)
|
|
if(n.eq.j)then
|
|
derR(ivar)=1.0d0
|
|
else
|
|
do t=1,ndesignparam
|
|
xj(t)=xsamplike(n,t)
|
|
enddo
|
|
derR(ivar)=corr_like(xi,xj,
|
|
& theta1,pdist1,ndesignparam)
|
|
endif
|
|
enddo
|
|
enddo
|
|
! then, horizontal elements
|
|
do i=k+1,nvars
|
|
ierr=icompactpostri(k,i)
|
|
do t=1,ndesignparam
|
|
theta1(t)=theta(ierr,t)
|
|
pdist1(t)=pdist(ierr,t)
|
|
enddo
|
|
ierr=icompactposoff(k,i)
|
|
call istartjstart(k,i,nvars,indexvars,istart,jstart)
|
|
do n=istart,istart+indexvars(i)-1
|
|
do t=1,ndesignparam
|
|
xi(t)=xsamplike(n,t)
|
|
enddo
|
|
do j=jstart,jstart+indexvars(k)-1
|
|
do t=1,ndesignparam
|
|
xj(t)=xsamplike(j,t)
|
|
enddo
|
|
ivar=icompactpostri(n,j)
|
|
derR(ivar)=0.5d0*dsqrt(covars(i)/covars(k))*
|
|
& correls(ierr)*corr_like(xi,xj,
|
|
& theta1,pdist1,ndesignparam)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
gradlike(icount)=0.0d0
|
|
do i=1,nlikepoints
|
|
call symmatindex(nlikepoints,i,kpvt)
|
|
sumi=0.0d0
|
|
do j=1,nlikepoints
|
|
sumi=sumi+derR(kpvt(j))*delvector(j)
|
|
enddo
|
|
gradlike(icount)=
|
|
& gradlike(icount)+delvector(i)*sumi
|
|
enddo
|
|
gradlike(icount)=-gradlike(icount)
|
|
!
|
|
!Now, the trace part
|
|
sumi=0.0d0
|
|
do i=1,nlikepoints
|
|
call symmatindex(nlikepoints,i,kpvt)
|
|
sumt=0.0d0
|
|
do j=1,nlikepoints
|
|
sumt=sumt+work(kpvt(j))*derR(kpvt(j))
|
|
enddo
|
|
sumi=sumi+sumt
|
|
enddo
|
|
|
|
! sumi is the trace
|
|
gradlike(icount)=gradlike(icount)+sumi
|
|
enddo
|
|
!--------------------------------------------------------------------
|
|
! ------- Section for calculating derivatives with respect to correls
|
|
do k=1,nvars*(nvars-1)/2
|
|
icount=icount+1
|
|
ierr=0
|
|
do j=1,nlikepoints
|
|
do i=1,j
|
|
ierr=ierr+1
|
|
derR(ierr)=0.0d0
|
|
enddo
|
|
enddo
|
|
call ivarjvaroff(k,nvars,indexvars,
|
|
& ivar,jvar,istart,jstart)
|
|
ierr=icompactpostri(ivar,jvar)
|
|
do t=1,ndesignparam
|
|
theta1(t)=theta(ierr,t)
|
|
pdist1(t)=pdist(ierr,t)
|
|
enddo
|
|
do i=istart,istart+indexvars(ivar)-1
|
|
do t=1,ndesignparam
|
|
xi(t)=xsamplike(i,t)
|
|
enddo
|
|
do j=jstart,jstart+indexvars(jvar)-1
|
|
do t=1,ndesignparam
|
|
xj(t)=xsamplike(j,t)
|
|
enddo
|
|
ierr=icompactpostri(i,j)
|
|
derR(ierr)=dsqrt(covars(ivar)*covars(jvar))*
|
|
& corr_like(xi,xj,theta1,pdist1,ndesignparam)
|
|
enddo
|
|
enddo
|
|
|
|
gradlike(icount)=0.0d0
|
|
do i=1,nlikepoints
|
|
call symmatindex(nlikepoints,i,kpvt)
|
|
sumi=0.0d0
|
|
do j=1,nlikepoints
|
|
sumi=sumi+derR(kpvt(j))*delvector(j)
|
|
enddo
|
|
gradlike(icount)=
|
|
& gradlike(icount)+delvector(i)*sumi
|
|
enddo
|
|
gradlike(icount)=-gradlike(icount)
|
|
|
|
!Now, the trace part
|
|
sumi=0.0d0
|
|
do i=1,nlikepoints
|
|
call symmatindex(nlikepoints,i,kpvt)
|
|
sumt=0.0d0
|
|
do j=1,nlikepoints
|
|
sumt=sumt+work(kpvt(j))*derR(kpvt(j))
|
|
enddo
|
|
sumi=sumi+sumt
|
|
enddo
|
|
! sumi is the trace
|
|
gradlike(icount)=gradlike(icount)+sumi
|
|
enddo
|
|
!--------------------------------------------------------------
|
|
ierr=0
|
|
return
|
|
end
|
|
|
|
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
|
|
subroutine solvebassystem(delta)
|
|
implicit none
|
|
include '../maxlikelihood/colikelihood.h'
|
|
!
|
|
! Calculating the best linear estimates of the coefficients of the basis functions,
|
|
! delta (the residual square error), and matrices needed for prediction and estimation
|
|
! of prediction error
|
|
|
|
integer kcopy(kbasis),inert(3)
|
|
double precision FR(kbasis,nlikepoints),
|
|
& a(kbasis*(kbasis+1)/2),rhs(kbasis),sum,
|
|
& acopy(max0(kbasis*(kbasis+1)/2,nlikepoints)),
|
|
& amat(kbasis*(kbasis+1)/2),
|
|
& rhscopy(kbasis),delta
|
|
|
|
integer n,j,t,i,mark,irowinc,index,nii
|
|
|
|
do i=1,nlikepoints
|
|
call symmatindex(nlikepoints,i,kpvt)
|
|
do j=1,kbasis
|
|
FR(j,i)=0.0d0
|
|
do t=1,nlikepoints
|
|
FR(j,i)=FR(j,i)+work(kpvt(t))*fbassamp(t,j)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
do i=1,kbasis
|
|
! right hand side of the linear system
|
|
rhs(i)=0.0d0
|
|
do j=1,nlikepoints
|
|
rhs(i)=rhs(i)+FR(i,j)*ysamplike(j)
|
|
enddo
|
|
! make a copy
|
|
rhscopy(i)=rhs(i)
|
|
enddo
|
|
|
|
! coefficient matrix, symmetric
|
|
t=0
|
|
do j=1,kbasis
|
|
do i=1,j
|
|
t=t+1
|
|
a(t)=0.0d0
|
|
do n=1,nlikepoints
|
|
a(t)=a(t)+FR(i,n)*fbassamp(n,j)
|
|
enddo
|
|
! make a copy
|
|
amat(t)=a(t)
|
|
enddo
|
|
enddo
|
|
|
|
! solve the basis function linear system
|
|
call dspfa(a,kbasis,kpvt,mark)
|
|
do i=1,kbasis*(kbasis+1)/2
|
|
acopy(i)=a(i)
|
|
enddo
|
|
do i=1,kbasis
|
|
kcopy(i)=kpvt(i)
|
|
enddo
|
|
if(isitlastcall.eqv..true..or.isitlastcall.eqv..TRUE.)then
|
|
! calculate the inverse only
|
|
do j=1,kbasis*(kbasis+1)/2
|
|
FtRinvF_inv(j)=a(j)
|
|
enddo
|
|
j=1
|
|
call dspdi(FtRinvF_inv,kbasis,kpvt,coeffbas,
|
|
& inert,delvector,j)
|
|
do i=1,kbasis
|
|
kpvt(i)=kcopy(i)
|
|
enddo
|
|
endif
|
|
call dspsl(a,kbasis,kpvt,rhs)
|
|
do i=1,kbasis
|
|
coeffbas(i)=rhs(i)
|
|
enddo
|
|
|
|
! one step improvement
|
|
do i=1,kbasis
|
|
call symmatindex(kbasis,i,kpvt)
|
|
sum=0.0d0
|
|
do t=1,kbasis
|
|
sum=sum+amat(kpvt(t))*coeffbas(t)
|
|
enddo
|
|
rhs(i)=(sum-rhscopy(i))*1.0d+9
|
|
enddo
|
|
call dspsl(acopy,kbasis,kcopy,rhs)
|
|
do i=1,kbasis
|
|
coeffbas(i)=coeffbas(i)-rhs(i)*1.0d-9
|
|
enddo
|
|
|
|
do i=1,nlikepoints
|
|
acopy(i)=0.0d0
|
|
do j=1,kbasis
|
|
acopy(i)=acopy(i)+coeffbas(j)*fbassamp(i,j)
|
|
enddo
|
|
acopy(i)=ysamplike(i)-acopy(i)
|
|
enddo
|
|
|
|
do t=1,nlikepoints
|
|
call symmatindex(nlikepoints,t,kpvt)
|
|
delvector(t)=0.0d0
|
|
do i=1,nlikepoints
|
|
delvector(t)=delvector(t)+work(kpvt(i))*acopy(i)
|
|
enddo
|
|
enddo
|
|
delta=0.0d0
|
|
do t=1,nlikepoints
|
|
delta=delta+delvector(t)*acopy(t)
|
|
enddo
|
|
|
|
if(isitlastcall.eqv..true..or.isitlastcall.eqv..TRUE.)then
|
|
! All parameters in the likelihood function have been optimized. Prepare
|
|
! matrices for MSE calculations
|
|
do i=1,kbasis
|
|
call symmatindex(kbasis,i,kpvt)
|
|
do j=1,nlikepoints
|
|
fbassamp(j,i)=0.0d0
|
|
do t=1,kbasis
|
|
fbassamp(j,i)=fbassamp(j,i)+FtRinvF_inv(kpvt(t))*
|
|
& FR(t,j)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
t=0
|
|
do j=1,nlikepoints
|
|
do i=1,j
|
|
t=t+1
|
|
work(t)=-work(t)
|
|
do n=1,kbasis
|
|
work(t)=work(t)+fbassamp(i,n)*FR(n,j)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
endif
|
|
return
|
|
end
|
|
!######################################################
|
|
|
|
DOUBLE PRECISION FUNCTION df1dimd_like(x)
|
|
INTEGER NMAX
|
|
double precision x
|
|
PARAMETER (NMAX=50)
|
|
CU USES gradlikefunc
|
|
INTEGER j,ncom
|
|
double precision df(NMAX),pcom(NMAX),xicom(NMAX),
|
|
& xt(NMAX)
|
|
,dummy
|
|
COMMON /f1com/ pcom,xicom,ncom
|
|
do 11 j=1,ncom
|
|
xt(j)=pcom(j)+x*xicom(j)
|
|
11 continue
|
|
call gradlikefunc(ncom,xt,dummy,df)
|
|
df1dimd_like=0.0d0
|
|
do 12 j=1,ncom
|
|
df1dimd_like=df1dimd_like+df(j)*xicom(j)
|
|
12 continue
|
|
return
|
|
END
|
|
C (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.
|
|
|
|
subroutine gradlikefunc(nplike,plike,fvalue,gradlike)
|
|
implicit none
|
|
integer nplike,ierr
|
|
double precision plike(nplike),fvalue,gradlike(nplike)
|
|
logical dogradient
|
|
|
|
dogradient=.true.
|
|
call colikelihoodgrad(nplike,plike,dogradient,
|
|
& fvalue,gradlike,ierr)
|
|
return
|
|
end
|
|
|
|
subroutine likelihoodfunc(nplike,plike,fvalue)
|
|
implicit none
|
|
integer nplike,ierr
|
|
double precision plike(nplike),fvalue,gradlike(nplike)
|
|
logical dogradient
|
|
|
|
dogradient=.false.
|
|
call colikelihoodgrad(nplike,plike,dogradient,
|
|
& fvalue,gradlike,ierr)
|
|
return
|
|
end
|
|
|
|
double precision function f1dim_like(x)
|
|
INTEGER NMAX
|
|
double precision x
|
|
PARAMETER (NMAX=50)
|
|
CU USES likelihood
|
|
INTEGER j,ncom
|
|
double precision pcom(NMAX),xicom(NMAX),xt(NMAX)
|
|
COMMON /f1com/ pcom,xicom,ncom
|
|
|
|
do 11 j=1,ncom
|
|
xt(j)=pcom(j)+x*xicom(j)
|
|
11 continue
|
|
call likelihoodfunc(ncom,xt,f1dim_like)
|
|
return
|
|
END
|
|
C (C) Copr. 1986-92 Numerical Recipes Software v%1jw#<0(9p#3.
|
|
|
|
!#################################################################
|
|
subroutine gradfunkmin(ndim,x,fatx,fjac)
|
|
implicit none
|
|
|
|
integer ndim,j
|
|
double precision x(ndim),fatx,fjac(ndim),h,eps,
|
|
& zero,temp,fxh
|
|
|
|
call likelihoodfunc(ndim,x,fatx)
|
|
|
|
eps = 1.0d-5
|
|
zero=0.0d0
|
|
|
|
do j = 1, ndim
|
|
temp = x(j)
|
|
h = eps*dabs(temp)
|
|
if (h .eq. zero) h = eps
|
|
x(j) = temp + h
|
|
call likelihoodfunc(ndim,x,fxh)
|
|
x(j) = temp
|
|
fjac(j) = (fxh - fatx)/h
|
|
enddo
|
|
|
|
return
|
|
end subroutine gradfunkmin
|
|
|
|
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
|
|
subroutine randpermut(npoints,ndim,x)
|
|
implicit none
|
|
!
|
|
! conduct random permutation
|
|
integer npoints,ndim,i,j,index
|
|
double precision x(npoints,ndim),xtemp(npoints),
|
|
& ran2
|
|
|
|
do i=1,ndim
|
|
do j=1,npoints
|
|
xtemp(j)=x(j,i)
|
|
enddo
|
|
do j=1,npoints
|
|
index=int(dble(npoints-j+1)*ran2()+1.0d0)
|
|
x(j,i)=xtemp(index)
|
|
xtemp(index)=xtemp(npoints-j+1)
|
|
enddo
|
|
enddo
|
|
return
|
|
end
|
|
|
|
c&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
|
|
c
|
|
c####################################################################
|
|
c random number generator
|
|
c
|
|
double precision function ran2()
|
|
implicit none
|
|
integer im1,im2,imm1,ia1,ia2,iq1,iq2,ir1,ir2,ntab,ndiv
|
|
double precision am,eps,rnmx
|
|
parameter(im1=2147483563,im2=2147483399,am=1.0d0/dble(im1),
|
|
&imm1=im1-1,ia1=40014,ia2=40692,iq1=53668,iq2=52774,ir1=
|
|
&12211,ir2=3791,ntab=32,ndiv=1+imm1/ntab,eps=1.2d-7,
|
|
&rnmx=1.0d0-eps)
|
|
integer idum2,j,k,iv(ntab),iy,idum
|
|
save iv,iy,idum2,idum
|
|
data idum2/123456789/,iv/ntab*0/,iy/0/,idum/-1/
|
|
|
|
if(idum.le.0) then
|
|
idum=max0(-idum,1)
|
|
idum2=idum
|
|
do 11 j=ntab+8,1,-1
|
|
k=idum/iq1
|
|
idum=ia1*(idum-k*iq1)-k*ir1
|
|
if(idum.lt.0) idum=idum+im1
|
|
if(j.le.ntab) iv(j)=idum
|
|
11 continue
|
|
iy=iv(1)
|
|
end if
|
|
k=idum/iq1
|
|
idum=ia1*(idum-k*iq1)-k*ir1
|
|
if(idum.lt.0)idum=idum+im1
|
|
k=idum2/iq2
|
|
idum2=ia2*(idum2-k*iq2)-k*ir2
|
|
if(idum2.lt.0) idum2=idum2+im2
|
|
j=1+iy/ndiv
|
|
iy=iv(j)-idum2
|
|
iv(j)=idum
|
|
if(iy.lt.1)iy=iy+imm1
|
|
ran2=dmin1(am*dble(iy),rnmx)
|
|
return
|
|
end
|
|
!######################################################
|