BeastJar {BeastJar} | R Documentation |
BeastJar
Description
Convenient packaging of the Bayesian Evolutionary Analysis Sampling Trees (BEAST) software package to facilitate Markov chain Monte Carlo sampling techniques including Hamiltonian Monte Carlo, bouncy particle sampling and zig-zag sampling.
Examples
# Example MCMC simulation using BEAST
#
# This function generates a Markov chain to sample from a simple normal distribution.
# It uses a random walk Metropolis kernel that is auto-tuning.
if (supportsJava8()) {
# Set seed
seed <- 123
rJava::J("dr.math.MathUtils")$setSeed(rJava::.jlong(seed));
# Set up simple model - Normal(mean = 1, sd = 2)
mean <- 1; sd <- 2
distribution <- rJava::.jnew("dr.math.distributions.NormalDistribution",
as.numeric(mean), as.numeric(sd))
model <- rJava::.jnew("dr.inference.distribution.DistributionLikelihood",
rJava::.jcast(distribution, "dr.math.distributions.Distribution"))
parameter <- rJava::.jnew("dr.inference.model.Parameter$Default", "p", 1.0,
as.numeric(-1.0 / 0.0), as.numeric(1.0 / 0.0))
model$addData(parameter)
# Construct posterior
dummy <- rJava::.jnew("dr.inference.model.DefaultModel",
rJava::.jcast(parameter, "dr.inference.model.Parameter"))
joint <- rJava::.jnew("java.util.ArrayList")
joint$add(rJava::.jcast(model, "dr.inference.model.Likelihood"))
joint$add(rJava::.jcast(dummy, "dr.inference.model.Likelihood"))
joint <- rJava::new(rJava::J("dr.inference.model.CompoundLikelihood"), joint)
# Specify auto-adapting random-walk Metropolis-Hastings transition kernel
operator <- rJava::.jnew("dr.inference.operators.RandomWalkOperator",
rJava::.jcast(parameter, "dr.inference.model.Parameter"),
0.75,
rJava::J(
"dr.inference.operators.RandomWalkOperator"
)$BoundaryCondition$reflecting,
1.0,
rJava::J("dr.inference.operators.AdaptationMode")$DEFAULT
)
schedule <- rJava::.jnew("dr.inference.operators.SimpleOperatorSchedule",
as.integer(1000), as.numeric(0.0))
schedule$addOperator(operator)
# Set up what features of posterior to log
subSampleFrequency <- 100
memoryFormatter <- rJava::.jnew("dr.inference.loggers.ArrayLogFormatter", FALSE)
memoryLogger <-
rJava::.jnew("dr.inference.loggers.MCLogger",
rJava::.jcast(memoryFormatter, "dr.inference.loggers.LogFormatter"),
rJava::.jlong(subSampleFrequency), FALSE)
memoryLogger$add(parameter)
# Execute MCMC
mcmc <- rJava::.jnew("dr.inference.mcmc.MCMC", "mcmc1")
mcmc$setShowOperatorAnalysis(FALSE)
chainLength <- 100000
mcmcOptions <- rJava::.jnew("dr.inference.mcmc.MCMCOptions",
rJava::.jlong(chainLength),
rJava::.jlong(10),
as.integer(1),
as.numeric(0.1),
TRUE,
rJava::.jlong(chainLength/100),
as.numeric(0.234),
FALSE,
as.numeric(1.0))
mcmc$init(mcmcOptions,
joint,
schedule,
rJava::.jarray(memoryLogger, contents.class = "dr.inference.loggers.Logger"))
mcmc$run()
# Summarize logged posterior quantities
traces <- memoryFormatter$getTraces()
trace <- traces$get(as.integer(1))
obj <- trace$getValues(as.integer(0),
as.integer(trace$getValueCount()))
sample <- rJava::J("dr.inference.trace.Trace")$toArray(obj)
outputStream <- rJava::.jnew("java.io.ByteArrayOutputStream")
printStream <- rJava::.jnew("java.io.PrintStream",
rJava::.jcast(outputStream, "java.io.OutputStream"))
rJava::J("dr.inference.operators.OperatorAnalysisPrinter")$showOperatorAnalysis(
printStream, schedule, TRUE)
operatorAnalysisString <- outputStream$toString("UTF8")
# Report auto-optimization information
cat(operatorAnalysisString)
# Report posterior quantities
c(mean(sample), sd(sample))
}
[Package BeastJar version 1.10.6 Index]