# inst/Devel/Colloc.MCMC1.R In CollocInfer: Collocation Inference for Dynamic Systems

```Colloc.MCMC = function(times,data,pars,coefs,lik,proc,prior,walk.var,cscale,nstep,in.meth='SplineEst',control.in=NULL)
{
cdims = dim(coefs)

parhist = matrix(0,nstep,length(pars))                     # Store parameter history
coefhist = matrix(0,nstep,length(coefs))
posthist = rep(0,nstep)
bestcoefhist = matrix(0,nstep,length(coefs))
acc = rep(0,nstep)

ei = eigen(walk.var)                                       # Random walk covariance
S = ei\$vectors%*%diag(sqrt(ei\$values))

pr = prior(pars)

f = SplineCoefsErr(coefs,times,data,lik,proc,pars,sgn=-1)  + pr  # Complete data log posterior

#    H = SplineCoefsDC2(coefs,times,data,lik,proc,pars,sgn=-1)
#    G = SplineCoefsDCDP(coefs,times,data,lik,proc,pars,sgn=-1)

#    dcdp = ginv(H)%*%G

tcoefs = coefs
#
#     ei = eigen(tH)
#     eta = rnorm(length(coefs))
#     logq = sum(eta^2)/2
#
#     coefs = coefs + ei\$vectors%*%diag(1/sqrt(ei\$values*(ei\$values>0)))%*%eta
#

parhist[1,] = pars                              # Initialize
bestcoefhist[1,]=tcoefs
coefhist[1,] = tcoefs

#  eta = 0        # cscale = 0 defaults to Collocation MCMC
logq = 0
for(i in 2:nstep){
print(c(i,pars))

delta = S %*% rnorm(length(pars))   # Random walk step

# Find the nearest parameters in the history

dists = apply( (parhist[1:(i-1),] - matrix(pars+delta,i-1,length(pars),byrow=TRUE))^2,1,sum)
whichrow = sort( dists,decreasing=FALSE,index.return=TRUE)\$ix[1]
tcoefs = matrix(bestcoefhist[whichrow,],cdims[1],cdims[2])

# Minimization

<<<<<<< .mine
res = inneropt(times=times,data=data,coefs=tcoefs,lik=lik,proc=proc,
=======
res = inneropt(times=times,data=data,coefs=matrix(as.vector(tcoefs), nrow=ncol(lik\$bvals),length(coefs)/ncol(lik\$bvals)),lik=lik,proc=proc,
>>>>>>> .r262
pars=pars+delta,in.meth=in.meth,control.in=control.in)

tcoefs = res\$coefs                                                  # Now a Gaussian error about the minimum
tH = SplineCoefsDC2(tcoefs,times,data,lik,proc,pars+delta,sgn=-1)

tpr = prior(pars+delta)

if(cscale>0){

for (j in 1:length(coefs)){

eta = rnorm(1)*cscale

logq = logq+dnorm(eta,log=TRUE)

ttcoefs=as.vector(tcoefs)
ttcoefs[j]=ttcoefs[j]+eta
ttcoefs=matrix(ttcoefs,dim(coefs))

posthist[i] =  SplineCoefsErr(ttcoefs,times,data,lik,proc,pars+delta,sgn=-1)

tf =  SplineCoefsErr(ttcoefs,times,data,lik,proc,pars+delta,sgn=-1)+ tpr - logq

if( runif(1) < exp( tf-f)){             # Acceptance probability
coefs=as.vector(coefs)
coefs[j]= as.vector(ttcoefs)[j]
coefs=matrix(coefs,dim(ttcoefs))
f = tf
acc[i] = 1

pars = pars + delta
pr = tpr
}
}
}

#   if(cscale>0){
#    ei = eigen(tH)
#    eta = rnorm(length(coefs))
#       logq = -sum(eta^2)/2 - sum(log(ei\$values[ei\$values>0])/2)
#       logq = sum(dnorm(eta,log=TRUE))

#      eta = cscale*ei\$vectors%*%diag(sqrt(ei\$values*(ei\$values>0)))%*%eta
#   }

#  tf =  SplineCoefsErr(as.vector(tcoefs)+eta,times,data,lik,proc,pars+delta,sgn=-1) + tpr - logq

print(c(tf,f,logq,pars+delta))
#    if( runif(1) < exp( tf-f)){             # Acceptance probability
#      coefs = as.vector(tcoefs)+eta
#      pars = pars + delta
#     pr = tpr

#    f = tf
#          H = SplineCoefsDC2(tcoefs,times,data,lik,proc,pars,sgn=-1)
#          G = SplineCoefsDCDP(tcoefs,times,data,lik,proc,pars,sgn=-1)
#          dcdp = ginv(H)%*%G

#   acc[i] = 1
#}

parhist[i,] = pars
coefhist[i,] = coefs
bestcoefhist[i,] = res\$coefs
}

return(list(parhist=parhist, coefhist=coefhist, acc = acc, posthist=posthist))
}
```

## Try the CollocInfer package in your browser

Any scripts or data that you put into this service are public.

CollocInfer documentation built on May 2, 2019, 4:03 a.m.