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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 | # 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))
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.