R/mcmcFit0.R

In thednainus/pangeaZA: Analyses and models for PANGEA data from South Africa (ZA).

library( phydynR)
library(BayesianTools )

#' Runs parallel mcmc chains in order to fit model0
#'
#' Priors are defined in the function
#'
#' @param BDTR A phydynR::DatedTree object
#' @param mode0 A list with model objects generated by generate_model
#' @param nc Integer number of CPUs and number of MCMC chains
#' @param iterations Length of MCMC chains
#' @param likelihood Switch for different coalescent likelihood approximations, passed to phydynR::colik
#' @value A BayesianTools::mcmcSamplerList
#' @export
mcmcfit0.0 <- function( BDTR , mod0 , nc =1 , iterations = 10, likelihood = 'PL2' , continueFits = NULL )
{
attach(mod0)
	ESTNAMES <- c( 'srcGrowth', 'src1980', 'i1980', 'importRate', 'rho2017', betaNames)
	ELOWER <-  c(0.01, 1.00, 0.01, 1/100, 1/24 , rep(0.001, length(betaNames) ) )
	EUPPER <-  c(0.10,  1e5, 99.9,   1/5, 1/1. , rep(2, length(betaNames) ) )

	priordensity <- function(theta)
	{
		names(theta) <- ESTNAMES
		o = dlnorm(theta['srcGrowth'], log(.035), 1 , log = TRUE ) +
		dlnorm(theta['src1980'], log(2e4), 1/4 , log = TRUE ) +
		dexp(theta['i1980'], 1/10 , log = TRUE ) +
		dlnorm(theta['importRate'], log(1/20), 1/4 , log = TRUE ) +
		sum( dlnorm(theta[betaNames], log(2/12), 1 , log = TRUE )  )
		unname( o )
	}

	sampler <- function(n = 1)
	{
		srcGrowth <- rlnorm ( n, log(.035), 1/8 )
		src1980 <- rlnorm( n, log(2e4), 1/2)
		i1980 <- rexp ( n, 1/10 )
		importRate <- rlnorm( n , log( 1/20 ), 1/4 )
		rho2017 <- runif( n, 1/8, 1/2)
		o = cbind( srcGrowth, src1980, i1980 , importRate, rho2017, matrix( rlnorm( n*length(betaNames), log(2/12), 1/4), nrow = n) )
		for ( i in 1:n){
			o[i,] <- pmax( ELOWER, o[i, ] )
			o[i, ] <- pmin( EUPPER, o[i,] )
		}
		o
	}

	prior <- createPrior(density = priordensity,
						 sampler = sampler,
						 lower = ELOWER ,
						 upper = EUPPER
	)

	co0 <- function(theta )
	{
		names(theta) <- ESTNAMES
	print(theta)
		p <- mod0$parms
		p[ names(theta)] <- theta
		x0 <- c(  src = p$src1980, i = p$i1980 , z =0)
		s <-  mod0$dm(p,  x0 , t0 = 1980, t1 = BDTR$maxSampleTime, res = 100 )
		y1 <- tail( s[[5]], 1 )[1,]
		zprior <- dnorm( y1['z'] / (y1['z'] + y1['i'] ) , .5, sd = .05, log = TRUE )
#~ browser()
		o = suppressWarnings( colik( BDTR, p , mod0$dm, x0 , t0 = 1980, res = 100, maxHeight = BDTR$maxSampleTime - 1980  , likelihood = likelihood) )
	print(o + zprior )
		o + zprior
	}

	bayesianSetup3 <- createBayesianSetup(likelihood = co0 , prior = prior, names = ESTNAMES , parallel = FALSE )
	settings3 = list(iterations = iterations, nrChains = 1, thin = 1)
	if(!is.null( continueFits))  {
		if ( inherits(continueFits, 'mcmcSamplerList' ) & ( length( continueFits ) != nc )){
			stop('nc must match length of continueFits')
		}
#~ browser()
		if ( inherits(continueFits, 'mcmcSamplerList' ) & nc > 1 ) {
			chs3 <- parallel::mclapply( 1:nc, function(i) {
			  runMCMC(continueFits[[i]], sampler = "DEzs", settings = settings3)
			}, mc.cores = nc )
			ch <- createMcmcSamplerList( chs3 )
		} else{
			ch <- runMCMC( continueFits, sampler = "DEzs", settings = settings3)
		}
	}else if ( nc > 1 ){
		chs3 <- parallel::mclapply( 1:nc, function(i) {
		 runMCMC(bayesianSetup = bayesianSetup3, sampler = "DEzs", settings = settings3)
		}, mc.cores = nc )
		ch <- createMcmcSamplerList( chs3 )
	} else {
		ch <- runMCMC(bayesianSetup = bayesianSetup3, sampler = "DEzs", settings = settings3)
	}
	ch
}