R/delayBUGS.r

In LobsterScience/bio.lobster: Lobster Stock Assessment Tools for DFO Maritimes

#' @export
delayBUGS <- function(model, input, priors, years, n = 60000, burn = 20000, thin = 10, nchains = 2, debug = F, wd=file.path( ecomod.directory, "lobster", "src" ), parameters=c(names(priors),'K','P','B','R','mu'),add.parameters=NULL,add.inits=NULL,pe=F,sw='bugs'){

	require(R2WinBUGS)
	require(rjags)
		start<-Sys.time()
		
		data=input	
		if(!missing(priors)){
			
			# Prepare priors for WinBUGS
			for(i in 1:length(priors)){
				if(priors[[i]]$d%in%c("dlnorm","dnorm"))priors[[i]]$b<-1/priors[[i]]$b^2
			#	if(priors[[i]]$d=="dgamma")priors[[i]]$b<-1/priors[[i]]$b
			}
			prior.dat<- data.frame(par=names(priors),do.call("rbind",lapply(priors,rbind)))
			prior.lst<-list()
			for(i in seq(1,nrow(prior.dat)*2,2)){
				prior.lst[[i]]<-prior.dat$a[[ceiling(i/2)]]
				prior.lst[[i+1]]<-prior.dat$b[[ceiling(i/2)]]
			}
			names(prior.lst)<-paste(rep(prior.dat$par,2)[order(rep(1:nrow(prior.dat),2))],rep(c('a','b'),nrow(prior.dat)),sep='.')
			 	
			data=c(prior.lst,input)
		}
		# parameters to save
		parameters<-c(parameters, add.parameters)
		
		# initial values
		inits<-list()
		prior.dat<-data.frame(do.call("rbind",priors))
		for(i in 1:nchains){
			inits[[i]]<-sapply(1:nrow(prior.dat),function(x){rep(unlist(prior.dat[paste('i',i,sep='')][x,]),prior.dat$l[[x]])})
			#browser()
			if(pe){
				prior.dat$l[prior.dat$l>1]<-input$NY
				inits[[i]]<-sapply(1:nrow(prior.dat),function(x){rep(unlist(prior.dat[paste('i',i,sep='')][x,]),prior.dat$l[[x]])})
			}
			names(inits[[i]])<-names(priors)
			inits[[i]]<-c(inits[[i]],list(P=rep(1,length(years))),add.inits[[i]])
		}
		
#		browser()
		if(missing(years))years=1:input$NY
		input$iyr<-years

		## Call to WinBUGS ##
		if(sw=='bugs'){
			tmp <- bugs(data = data, inits, parameters.to.save = parameters, model.file = file.path(wd,"bugs",paste(model,".bug",sep="")), n.chains = nchains, n.iter = n, n.burnin = burn, n.thin = thin, bugs.directory = file.path("C:","WinBUGS14"), debug = debug)
			out.lst<-list(data=input, sims.list=tmp$sims.list,median=tmp$median,mean=tmp$mean,summary=tmp$summary)
			
		}
		## Call to JAGS ##
		if(sw=='jags'){
			m = jags.model( file=file.path(wd,"bugs",paste(model,".bug",sep="")), data=data, n.chains=nchains, n.adapt=burn ) 
			tmp = jags.samples(m, variable.names=parameters, n.iter=n, thin=thin) # sample from posterior
			sims.list<-list()
			for(i in 1:length(tmp)){sims.list[[i]]<-sapply(1:dim(tmp[[i]])[1],function(y){rowMeans(tmp[[i]][y,,])})}
			names(sims.list)<-names(tmp)
			out.lst<-list(data=input, sims.list=sims.list, median=lapply(tmp,apply,1,median),mean=lapply(tmp,apply,1,mean))
		}
		#browser()


		print(Sys.time()-start)

		rm(tmp)

		return(out.lst)	
}