Mamba: Markov chain Monte Carlo (MCMC) for Bayesian analysis in julia¶
Version:  0.10.1 

Requires:  julia releases 0.5.0 or later 
Date:  Mar 26, 2018 
Maintainer:  Brian J Smith (brianjsmith@uiowa.edu) 
Contributors:  Benjamin Deonovic (benjamindeonovic@uiowa.edu), Brian J Smith (brianjsmith@uiowa.edu), and others 
Web site:  https://github.com/brianjsmith/Mamba.jl 
License:  MIT 
Overview¶
Purpose¶
Mamba is an open platform for the implementation and application of MCMC methods to perform Bayesian analysis in julia. The package provides a framework for (1) specification of hierarchical models through stated relationships between data, parameters, and statistical distributions; (2) blockupdating of parameters with samplers provided, defined by the user, or available from other packages; (3) execution of sampling schemes; and (4) posterior inference. It is intended to give users access to all levels of the design and implementation of MCMC simulators to particularly aid in the development of new methods.
Several software options are available for MCMC sampling of Bayesian models. Individuals who are primarily interested in data analysis, unconcerned with the details of MCMC, and have models that can be fit in JAGS, Stan, or OpenBUGS are encouraged to use those programs. Mamba is intended for individuals who wish to have access to lowerlevel MCMC tools, are knowledgeable of MCMC methodologies, and have experience, or wish to gain experience, with their application. The package also provides standalone convergence diagnostics and posterior inference tools, which are essential for the analysis of MCMC output regardless of the software used to generate it.
Features¶
An interactive and extensible interface.
Support for a wide range of model and distributional specifications.
An environment in which all interactions with the software are made through a single, interpreted programming language.
 Any julia operator, function, type, or package can be used for model specification.
 Custom distributions and samplers can be written in julia to extend the package.
Directed acyclic graph representations of models.
Arbitrary blocking of model parameters and designation of blockspecific samplers.
Samplers that can be used with the included simulation engine or apart from it, including
 adaptive Metropolis within Gibbs and multivariate Metropolis,
 approximate Bayesian computation,
 binary,
 Hamiltonian Monte Carlo (simple and NoUTurn),
 simplex, and
 slice samplers.
Automatic parallel execution of parallel MCMC chains on multiprocessor systems.
Restarting of chains.
Commandline access to all package functionality, including its simulation API.
Convergence diagnostics: Gelman, Rubin, and Brooks; Geweke; Heidelberger and Welch; Raftery and Lewis.
Posterior summaries: moments, quantiles, HPD, crosscovariance, autocorrelation, MCSE, ESS.
Gadfly plotting: trace, density, running mean, autocorrelation.
Importing of sampler output saved in the CODA file format.
Runtime performance on par with compiled MCMC software.
Contents¶
Introduction¶
MCMC Software¶
Markov chain Monte Carlo (MCMC) methods are a class of algorithms for simulating autocorrelated draws from probability distributions [13][28][39][78]. They are widely used to obtain empirical estimates for and make inference on multidimensional distributions that often arise in Bayesian statistical modelling, computational physics, and computational biology. Because MCMC provides estimates of distributions of interest, and is not limited to point estimates and asymptotic standard errors, it facilitates wide ranges of inferences and provides for more realistic prediction errors. An MCMC algorithm can be devised for any probability model. Implementations of algorithms are computational in nature, with the resources needed to execute algorithms directly related to the dimensionality of their associated problems. Rapid increases in computing power and emergence of MCMC software have enabled models of increasing complexity to be fit. For all its advantages, MCMC is regarded as one of the most important developments and powerful tools in modern statistical computing.
Several software programs provide Bayesian modelling via MCMC. Programs range from those designed for general model fitting to those for specific models. WinBUGS, its opensource incarnation OpenBUGS, and the ‘BUGS’ clone Just Another Gibbs Sampler (JAGS) are among the most widely used programs for general model fitting [57][71]. These three provide similar programming syntaxes with which users can specify statistical models by simply stating relationships between data, parameters, and statistical distributions. Once a model is specified, the programs automatically formulate an MCMC sampling scheme to simulate parameter values from their posterior distribution. All aforementioned tasks can be accomplished with minimal programming and without specific knowledge of MCMC methodology. Users who are adept at both and so inclined can write software modules to add new distributions and samplers to OpenBUGS and JAGS [94][96]. General model fitting is also available with the MCMC procedure found in the SAS/STAT ® software [49]. Stan is a newer and similarinscope program worth noting for its accessible syntax and automatically tuned Hamiltonian Monte Carlo sampling scheme [98]. PyMC is a Pythonbased program that allows all modelling tasks to be accomplished in its native language, and gives users more handson access to model and sampling scheme specifications [70]. Programs like GRIMS [66] and LaplacesDemon [99] represent another class of programs that fit general models. In their approaches, users work with the functional forms of (unnormalized) probability densities directly, rather a domain specific modelling language (DSL), for model specification. Examples of programs for specific models can be found in the R catalogue of packages. For instance, the arm package provides Bayesian inference for generalized linear, ordered logistic or probit, and mixedeffects regression models [34], MCMCpack fits a wide range of models commonly encountered in the social and behavioral sciences [60], and many others that are more focused on specific classes of models can be found in the “Bayesian Inference” task view on the Comprehensive R Archive Network [69].
The Mamba Package¶
Mamba [88] is a julia [4] package designed for general Bayesian model fitting via MCMC. Like OpenBUGS and JAGS, it supports a wide range of model and distributional specifications, and provides a syntax for model specification. Unlike those two, and like PyMC, Mamba provides a unified environment in which all interactions with the software are made through a single, interpreted language. Any julia operator, function, type, or package can be used for model specification; and custom distributions and samplers can be written in julia to extend the package. Conversely, interactions with and extensions to OpenBUGS and JAGS can involve three different programming environments — R wrappers used to call the programs, their DSLs, and the underlying implementations in Component Pascal and C++. Advantages of a unified environment include more flexible model specification; tighter integration with supplied functions for convergence diagnostics and posterior inference; and faster development, testing, and debugging of extensions. Advantages of the BUGS DSLs include more concise model specification and facilitation of automated sampling scheme formulation. Indeed, sampling schemes must be selected manually in the initial release of Mamba. Nevertheless, Mamba holds other distinct advantages over existing offerings. In particular, it provides arbitrary blocking of model parameters and designation of blockspecific samplers; samplers that can be used with the included simulation engine or apart from it; and commandline access to all package functionality, including its simulation API. Likewise, advantages of the julia language include its familiar syntax, focus on technical computing, and benchmarks showing it to be one or more orders of magnitude faster than R and Python [3]. Finally, the intended audience for Mamba includes individuals interested in programming in julia; who wish to have lowlevel access to model design and implementation; and, in some cases, are able to derive full conditional distributions of model parameters (up to normalizing constants).
Mamba allows for the implementation of an MCMC sampling scheme to simulate draws for a set of Bayesian model parameters from their joint posterior distribution. The package supports the general Gibbs [30][35] scheme outlined in the algorithm below. In its implementation with the package, the user may specify blocks of parameters and corresponding functions to sample each from its full conditional distribution . Simulation performance (efficiency and runtime) can be affected greatly by the choice of blocking scheme and sampling functions. For some models, an optimal choice may not be obvious, and different choices may need to be tried to find one that gives a desired level of performance. This can be a timeconsuming process. The Mamba package provides a set of julia types and method functions to facilitate the specification of different schemes and functions. Supported sampling functions include those provided by the package, userdefined functions, and functions from other packages; thus providing great flexibility with respect to sampling methods. Furthermore, a sampling engine is provided to save the user from having to implement tasks common to all MCMC simulators. Therefore, time and energy can be focused on implementation aspects that most directly affect performance.
A summary of the steps involved in using the package to perform MCMC simulation for a Bayesian model is given below.
Decide on names to use for julia objects that will represent the model data structures and parameters (). For instance, the Tutorial section describes a linear regression example in which predictor and response are represented by objects
x
andy
, and regression parameters , , and by objectsb0
,b1
, ands2
.Create a dictionary to store all structures considered to be fixed in the simulation; e.g., the
line
dictionary in the regression example.Specify the model using the constructors described in the MCMC Types section, to create the following:
 A
Stochastic
object for each model term that has a distributional specification. This includes parameters and data, such as the regression parametersb0
,b1
, ands2
that have prior distributions andy
that has a likelihood specification. A vector of
Sampler
objects containing supplied, userdefined, or external functions for sampling each parameter block . A
Model
object from the resulting stochastic nodes and sampler vector.Simulate parameter values with the
mcmc()
function.Use the MCMC output to check convergence and perform posterior inference.
Tutorial¶
The complete source code for the examples contained in this tutorial can be obtained here
.
Bayesian Linear Regression Model¶
In the sections that follow, the Bayesian simple linear regression example from the BUGS 0.5 manual [89] is used to illustrate features of the package. The example describes a regression relationship between observations and that can be expressed as
with prior distribution specifications
where , and is a design matrix with an intercept vector of ones in the first column and in the second. Primary interest lies in making inference about the , , and parameters, based on their posterior distribution. A computational consideration in this example is that inference cannot be obtain from the joint posterior directly because of its nonstandard form, derived below up to a constant of proportionality.
A common alternative is to make approximate inference based on parameter values simulated from the posterior with MCMC methods.
Model Specification¶
Nodes¶
In the Mamba package, terms that appear in the Bayesian model specification are referred to as nodes. Nodes are classified as one of three types:
 Stochastic nodes are any model terms that have likelihood or prior distributional specifications. In the regression example, , , and are stochastic nodes.
 Logical nodes are terms, like , that are deterministic functions of other nodes.
 Input nodes are any remaining model terms () and are considered to be fixed quantities in the analysis.
Note that the node has both a distributional specification and is a fixed quantity. It is designated a stochastic node in Mamba because of its distributional specification. This allows for the possibility of model terms with distributional specifications that are a mix of observed and unobserved elements, as in the case of missing values in response vectors.
Model implementation begins by instantiating stochastic and logical nodes using the Mamba–supplied constructors Stochastic
and Logical
. Stochastic and logical nodes for a model are defined with a call to the Model
constructor. The model constructor formally defines and assigns names to the nodes. Userspecified names are given on the lefthand sides of the arguments to which Logical
and Stochastic
nodes are passed.
using Mamba
## Model Specification
model = Model(
y = Stochastic(1,
(mu, s2) > MvNormal(mu, sqrt(s2)),
false
),
mu = Logical(1,
(xmat, beta) > xmat * beta,
false
),
beta = Stochastic(1,
() > MvNormal(2, sqrt(1000))
),
s2 = Stochastic(
() > InverseGamma(0.001, 0.001)
)
)
A single integer value for the first Stochastic
constructor argument indicates that the node is an array of the specified dimension. Absence of an integer value implies a scalar node. The next argument is a function that may contain any valid julia code. Functions should be defined to take, as their arguments, the inputs upon which their nodes depend and, for stochastic nodes, return distribution objects or arrays of objects compatible with the Distributions package [2]. Such objects represent the nodes’ distributional specifications. An optional boolean argument after the function can be specified to indicate whether values of the node should be monitored (saved) during MCMC simulations (default: true
).
Stochastic functions must return a single distribution object that can accommodate the dimensionality of the node, or return an array of (univariate or multivariate) distribution objects of the same dimension as the node. Examples of alternative, but equivalent, prior distributional specifications for the beta
node are shown below.
# Case 1: Multivariate Normal with independence covariance matrix
beta = Stochastic(1,
() > MvNormal(2, sqrt(1000))
)
# Case 2: One common univariate Normal
beta = Stochastic(1,
() > Normal(0, sqrt(1000))
)
# Case 3: Array of univariate Normals
beta = Stochastic(1,
() > UnivariateDistribution[Normal(0, sqrt(1000)), Normal(0, sqrt(1000))]
)
# Case 4: Array of univariate Normals
beta = Stochastic(1,
() > UnivariateDistribution[Normal(0, sqrt(1000)) for i in 1:2]
)
Case 1 is one of the multivariate normal distributions available in the Distributions package, and the specification used in the example model implementation. In Case 2, a single univariate normal distribution is specified to imply independent priors of the same type for all elements of beta
. Cases 3 and 4 explicitly specify a univariate prior for each element of beta
and allow for the possibility of differences among the priors. Both return arrays of Distribution objects, with the last case automating the specification of array elements.
In summary, y
and beta
are stochastic vectors, s2
is a stochastic scalar, and mu
is a logical vector. We note that the model could have been implemented without mu
. It is included here primarily to illustrate use of a logical node. Finally, note that nodes y
and mu
are not being monitored.
Sampling Schemes¶
The package provides a flexible system for the specification of schemes to sample stochastic nodes. Arbitrary blocking of nodes and designation of blockspecific samplers is supported. Furthermore, blockupdating of nodes can be performed with samplers provided, defined by the user, or available from other packages. Schemes are specified as vectors of Sampler
objects. Constructors are provided for several popular sampling algorithms, including adaptive Metropolis, NoUTurn (NUTS), and slice sampling. Example schemes are shown below. In the first one, NUTS is used for the sampling of beta
and slice for s2
. The two nodes are block together in the second scheme and sampled jointly with NUTS.
## Hybrid NoUTurn and Slice Sampling Scheme
scheme1 = [NUTS(:beta),
Slice(:s2, 3.0)]
## NoUTurn Sampling Scheme
scheme2 = [NUTS([:beta, :s2])]
Additionally, users are free to create their own samplers with the generic Sampler
constructor. This is particularly useful in settings were full conditional distributions are of standard forms for some nodes and can be sampled from directly. Such is the case for the full conditional of which can be written as
where and , and is recognizable as multivariate normal. Likewise,
whose form is inverse gamma with shape and scale parameters. A userdefined sampling scheme to generate draws from these full conditionals is constructed below. Another example is given in the Pollution example for the sampling of multiple nodes.
## UserDefined Samplers
Gibbs_beta = Sampler([:beta],
(beta, s2, xmat, y) >
begin
beta_mean = mean(beta.distr)
beta_invcov = invcov(beta.distr)
Sigma = inv(Symmetric(xmat' * xmat / s2 + beta_invcov))
mu = Sigma * (xmat' * y / s2 + beta_invcov * beta_mean)
rand(MvNormal(mu, Sigma))
end
)
Gibbs_s2 = Sampler([:s2],
(mu, s2, y) >
begin
a = length(y) / 2.0 + shape(s2.distr)
b = sumabs2(y  mu) / 2.0 + scale(s2.distr)
rand(InverseGamma(a, b))
end
)
## UserDefined Sampling Scheme
scheme3 = [Gibbs_beta, Gibbs_s2]
In these samplers, the respective MvNormal(2, sqrt(1000))
and InverseGamma(0.001, 0.001)
priors on stochastic nodes beta
and s2
are accessed directly through the distr
fields. Features of the Distributions objects returned by beta.distr
and s2.distr
can, in turn, be extracted with method functions defined in that package or through their own fields. For instance, mean(beta.distr)
and invcov(beta.distr)
apply method functions to extract the mean vector and inversecovariance matrix of the beta
prior. Whereas, shape(s2.distr)
and scale(s2.distr)
extract the shape and scale parameters from fields of the inversegamma prior. Distributions method functions can be found in that package’s documentation; whereas, fields are found in the source code.
When possible to do so, direct sampling from full conditions is often preferred in practice because it tends to be more efficient than generalpurpose algorithms. Schemes that mix the two approaches can be used if full conditionals are available for some model parameters but not for others. Once a sampling scheme is formulated in Mamba, it can be assigned to an existing model with a call to the setsamplers!
function.
## Sampling Scheme Assignment
setsamplers!(model, scheme1)
Alternatively, a predefined scheme can be passed in to the Model
constructor at the time of model implementation as the value to its samplers
argument.
Directed Acyclic Graphs¶
One of the internal structures created by Model
is a graph representation of the model nodes and their associations. Graphs are managed internally with the LightGraphs package [11]. Like OpenBUGS, JAGS, and other BUGS clones, Mamba fits models whose nodes form a directed acyclic graph (DAG). A DAG is a graph in which nodes are connected by directed edges and no node has a path that loops back to itself. With respect to statistical models, directed edges point from parent nodes to the child nodes that depend on them. Thus, a child node is independent of all others, given its parents.
The DAG representation of a Model
can be printed out at the commandline or saved to an external file in a format that can be displayed with the Graphviz software.
## Graph Representation of Nodes
>>> draw(model)
digraph MambaModel {
"mu" [shape="diamond", fillcolor="gray85", style="filled"];
"mu" > "y";
"xmat" [shape="box", fillcolor="gray85", style="filled"];
"xmat" > "mu";
"beta" [shape="ellipse"];
"beta" > "mu";
"s2" [shape="ellipse"];
"s2" > "y";
"y" [shape="ellipse", fillcolor="gray85", style="filled"];
}
>>> draw(model, filename="lineDAG.dot")
Either the printed or saved output can be passed manually to the Graphviz software to plot a visual representation of the model. If julia is being used with a frontend that supports inline graphics, like IJulia [50], and the GraphViz julia package [26] is installed, then the following code will plot the graph automatically.
using GraphViz
>>> display(Graph(graph2dot(model)))
A generated plot of the regression model graph is show in the figure below.
Stochastic, logical, and input nodes are represented by ellipses, diamonds, and rectangles, respectively. Graycolored nodes are ones designated as unmonitored in MCMC simulations. The DAG not only allows the user to visually check that the model specification is the intended one, but is also used internally to check that nodal relationships are acyclic.
MCMC Simulation¶
Data¶
For the example, observations are stored in a julia dictionary defined in the code block below. Included are predictor and response vectors :x
and :y
as well as a design matrix :xmat
corresponding to the model matrix .
## Data
line = Dict{Symbol, Any}(
:x => [1, 2, 3, 4, 5],
:y => [1, 3, 3, 3, 5]
)
line[:xmat] = [ones(5) line[:x]]
Initial Values¶
A julia vector of dictionaries containing initial values for all stochastic nodes must be created. The dictionary keys should match the node names, and their values should be vectors whose elements are the same type of structures as the nodes. Three sets of initial values for the regression example are created in with the following code.
## Initial Values
inits = [
Dict{Symbol, Any}(
:y => line[:y],
:beta => rand(Normal(0, 1), 2),
:s2 => rand(Gamma(1, 1))
)
for i in 1:3
]
Initial values for y
are those in the observed response vector. Likewise, the node is not updated in the sampling schemes defined earlier and thus retains its initial values throughout MCMC simulations. Initial values are generated for beta
from a normal distribution and for s2
from a gamma distribution.
MCMC Engine¶
MCMC simulation of draws from the posterior distribution of a declared set of model nodes and sampling scheme is performed with the mcmc()
function. As shown below, the first three arguments are a Model
object, a dictionary of values for input nodes, and a dictionary vector of initial values. The number of draws to generate in each simulation run is given as the fourth argument. The remaining arguments are named such that burnin
is the number of initial values to discard to allow for convergence; thin
defines the interval between draws to be retained in the output; and chains
specifies the number of times to run the simulator. Results are retuned as Chains
objects on which methods for posterior inference are defined.
## MCMC Simulations
setsamplers!(model, scheme1)
sim1 = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3)
setsamplers!(model, scheme2)
sim2 = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3)
setsamplers!(model, scheme3)
sim3 = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3)
Parallel Computing¶
The simulation of multiple chains will be executed in parallel automatically if julia is running in multiprocessor mode on a multiprocessor system. Multiprocessor mode can be started with the command line argument julia p n
, where n
is the number of available processors. See the julia documentation on parallel computing for details.
Posterior Inference¶
Convergence Diagnostics¶
Checks of MCMC output should be performed to assess convergence of simulated draws to the posterior distribution. Checks can be performed with a variety of methods. The diagnostic of Gelman, Rubin, and Brooks [33][12] is one method for assessing convergence of posterior mean estimates. Values of the diagnostic’s potential scale reduction factor (PSRF) that are close to one suggest convergence. As a ruleofthumb, convergence is rejected if the 97.5 percentile of a PSRF is greater than 1.2.
>>> gelmandiag(sim1, mpsrf=true, transform=true) > showall
Gelman, Rubin, and Brooks Diagnostic:
PSRF 97.5%
beta[1] 1.009 1.010
beta[2] 1.009 1.010
s2 1.008 1.016
Multivariate 1.006 NaN
The diagnostic of Geweke [37] tests for nonconvergence of posterior mean estimates. It provides chainspecific test pvalues. Convergence is rejected for significant pvalues, like those obtained for s2
.
>>> gewekediag(sim1) > showall
Geweke Diagnostic:
First Window Fraction = 0.1
Second Window Fraction = 0.5
Zscore pvalue
beta[1] 1.237 0.2162
beta[2] 1.568 0.1168
s2 1.710 0.0872
Zscore pvalue
beta[1] 1.457 0.1452
beta[2] 1.752 0.0797
s2 1.428 0.1534
Zscore pvalue
beta[1] 0.550 0.5824
beta[2] 0.440 0.6597
s2 0.583 0.5596
The diagnostic of Heidelberger and Welch [47] tests for nonconvergence (nonstationarity) and whether ratios of estimation interval halfwidths to means are within a target ratio. Stationarity is rejected (0) for significant test pvalues. Halfwidth tests are rejected (0) if observed ratios are greater than the target, as is the case for s2
and beta[1]
.
>>> heideldiag(sim1) > showall
Heidelberger and Welch Diagnostic:
Target Halfwidth Ratio = 0.1
Alpha = 0.05
Burnin Stationarity pvalue Mean Halfwidth Test
beta[1] 251 1 0.0680 0.57366275 0.053311283 1
beta[2] 738 1 0.0677 0.81285744 0.015404173 1
s2 738 1 0.0700 1.00825202 0.094300432 1
Burnin Stationarity pvalue Mean Halfwidth Test
beta[1] 251 1 0.1356 0.6293320 0.065092099 0
beta[2] 251 1 0.0711 0.7934633 0.019215278 1
s2 251 1 0.4435 1.4635400 0.588158612 0
Burnin Stationarity pvalue Mean Halfwidth Test
beta[1] 251 1 0.0515 0.5883602 0.058928034 0
beta[2] 1225 1 0.1479 0.8086080 0.018478999 1
s2 251 1 0.6664 0.9942853 0.127959523 0
The diagnostic of Raftery and Lewis [74][75] is used to determine the number of iterations required to estimate a specified quantile within a desired degree of accuracy. For example, below are required total numbers of iterations, numbers to discard as burnin sequences, and thinning intervals for estimating 0.025 quantiles so that their estimated cumulative probabilities are within 0.025±0.005 with probability 0.95.
>>> rafterydiag(sim1) > showall
Raftery and Lewis Diagnostic:
Quantile (q) = 0.025
Accuracy (r) = 0.005
Probability (s) = 0.95
Thinning Burnin Total Nmin Dependence Factor
beta[1] 2 267 17897 3746 4.7776295
beta[2] 2 267 17897 3746 4.7776295
s2 2 257 8689 3746 2.3195408
Thinning Burnin Total Nmin Dependence Factor
beta[1] 4 271 2.1759x104 3746 5.8085958
beta[2] 4 275 2.8795x104 3746 7.6868660
s2 2 257 8.3450x103 3746 2.2277096
Thinning Burnin Total Nmin Dependence Factor
beta[1] 2 269 2.0647x104 3746 5.5117459
beta[2] 2 263 1.4523x104 3746 3.8769354
s2 2 255 7.8770x103 3746 2.1027763
More information on the diagnostic functions can be found in the Convergence Diagnostics section.
Posterior Summaries¶
Once convergence has been assessed, sample statistics may be computed on the MCMC output to estimate features of the posterior distribution. The StatsBase package [53] is utilized in the calculation of many posterior estimates. Some of the available posterior summaries are illustrated in the code block below.
## Summary Statistics
>>> describe(sim1)
Iterations = 252:10000
Thinning interval = 2
Chains = 1,2,3
Samples per chain = 4875
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
beta[1] 0.5971183 1.14894446 0.0095006014 0.016925598 4607.9743
beta[2] 0.8017036 0.34632566 0.0028637608 0.004793345 4875.0000
s2 1.2203777 2.00876760 0.0166104638 0.101798287 389.3843
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
beta[1] 1.74343373 0.026573102 0.59122696 1.1878720 2.8308472
beta[2] 0.12168742 0.628297573 0.80357822 0.9719441 1.5051573
s2 0.17091385 0.383671702 0.65371989 1.2206381 6.0313970
## Highest Posterior Density Intervals
>>> hpd(sim1)
95% Lower 95% Upper
beta[1] 1.75436235 2.8109571
beta[2] 0.09721501 1.4733163
s2 0.08338409 3.8706865
## CrossCorrelations
>>> cor(sim1)
beta[1] beta[2] s2
beta[1] 1.000000000 0.905245029 0.027467317
beta[2] 0.905245029 1.000000000 0.024489462
s2 0.027467317 0.024489462 1.000000000
## LagAutocorrelations
>>> autocor(sim1)
Lag 2 Lag 10 Lag 20 Lag 100
beta[1] 0.24521566 0.021411797 0.0077424153 0.044989417
beta[2] 0.20402485 0.019107846 0.0033980453 0.053869216
s2 0.85931351 0.568056917 0.3248136852 0.024157524
Lag 2 Lag 10 Lag 20 Lag 100
beta[1] 0.28180489 0.031007672 0.03930888 0.0394895028
beta[2] 0.25905976 0.017946010 0.03613043 0.0227758214
s2 0.92905843 0.761339226 0.58455868 0.0050215824
Lag 2 Lag 10 Lag 20 Lag 100
beta[1] 0.38634357 0.0029361782 0.032310111 0.0028806786
beta[2] 0.32822879 0.0056670786 0.020100663 0.0062622517
s2 0.68812720 0.2420402859 0.080495078 0.0290205896
## State Space Change Rate (per Iteration)
>>> changerate(sim1)
Change Rate
beta[1] 0.844
beta[2] 0.844
s2 1.000
Multivariate 1.000
## Deviance Information Criterion
>>> dic(sim1)
DIC Effective Parameters
pD 13.828540 1.1661193
pV 22.624104 5.5639015
Output Subsetting¶
Additionally, sampler output can be subsetted to perform posterior inference on select iterations, parameters, and chains.
## Subset Sampler Output
>>> sim = sim1[1000:5000, ["beta[1]", "beta[2]"], :]
>>> describe(sim)
Iterations = 1000:5000
Thinning interval = 2
Chains = 1,2,3
Samples per chain = 2001
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
beta[1] 0.62445845 1.0285709 0.013275474 0.023818436 1864.8416
beta[2] 0.79392648 0.3096614 0.003996712 0.006516677 2001.0000
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
beta[1] 1.53050898 0.076745702 0.61120944 1.2174641 2.6906753
beta[2] 0.18846617 0.618849048 0.79323126 0.9619767 1.4502109
File I/O¶
For cases in which it is desirable to store sampler output in external files for processing in future julia sessions, read and write methods are provided.
## Write to and Read from an External File
write("sim1.jls", sim1)
sim1 = read("sim1.jls", ModelChains)
Restarting the Sampler¶
Convergence diagnostics or posterior summaries may indicate that additional draws from the posterior are needed for inference. In such cases, the mcmc()
function can be used to restart the sampler with previously generated output, as illustrated below.
## Restart the Sampler
>>> sim = mcmc(sim1, 5000)
>>> describe(sim)
Iterations = 252:15000
Thinning interval = 2
Chains = 1,2,3
Samples per chain = 7375
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
beta[1] 0.59655228 1.1225920 0.0075471034 0.014053505 6380.79199
beta[2] 0.80144540 0.3395731 0.0022829250 0.003954871 7372.28048
s2 1.18366563 1.8163096 0.0122109158 0.070481708 664.08995
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
beta[1] 1.70512374 0.031582137 0.58989089 1.1783924 2.8253668
beta[2] 0.12399073 0.630638800 0.80358526 0.9703569 1.4939817
s2 0.17075261 0.382963160 0.65372440 1.2210168 5.7449800
Plotting¶
Plotting of sampler output in Mamba is based on the Gadfly package [51]. Summary plots can be created and written to files using the plot
and draw
functions.
## Default summary plot (trace and density plots)
p = plot(sim1)
## Write plot to file
draw(p, filename="summaryplot.svg")
The plot
function can also be used to make autocorrelation and running means plots. Legends can be added with the optional legend
argument. Arrays of plots can be created and passed to the draw
function. Use nrow
and ncol
to determine how many rows and columns of plots to include in each drawing.
## Autocorrelation and running mean plots
p = plot(sim1, [:autocor, :mean], legend=true)
draw(p, nrow=3, ncol=2, filename="autocormeanplot.svg")
## Pairwise contour plots
p = plot(sim1, :contour)
draw(p, nrow=2, ncol=2, filename="contourplot.svg")
Computational Performance¶
Computing runtimes were recorded for different sampling algorithms applied to the regression example. Runs wer performed on a desktop computer with an Intel i52500 CPU @ 3.30GHz. Results are summarized in the table below. Note that these are only intended to measure the raw computing performance of the package, and do not account for different efficiencies in output generated by the sampling algorithms.
Adaptive Metropolis  Slice  

Within Gibbs  Multivariate  Gibbs  NUTS  Within Gibbs  Multivariate 
16,700  11,100  27,300  2,600  13,600  17,600 
Development and Testing¶
Commandline access is provided for all package functionality to aid in the development and testing of models. Examples of available functions are shown in the code block below. Documentation for these and other related functions can be found in the MCMC Types section.
## Development and Testing
setinputs!(model, line) # Set input node values
setinits!(model, inits[1]) # Set initial values
setsamplers!(model, scheme1) # Set sampling scheme
showall(model) # Show detailed node information
logpdf(model, 1) # Logdensity sum for block 1
logpdf(model, 2) # Block 2
logpdf(model) # All blocks
sample!(model, 1) # Sample values for block 1
sample!(model, 2) # Block 2
sample!(model) # All blocks
In this example, functions setinputs!
, setinits!
, and setsampler!
allow the user to manually set the input node values, the initial values, and the sampling scheme form the model
object, and would need to be called prior to logpdf
and sample!
. Updated model objects should be returned when called; otherwise, a problem with the supplied values may exist. Method showall
prints a detailed summary of all model nodes, their values, and attributes; logpdf
sums the logdensities over nodes associated with a specified sampling block (second argument); and sample!
generates an MCMC sample of values for the nodes. Nonnumeric results may indicate problems with distributional specifications in the second case or with sampling functions in the last case. The block arguments are optional; and, if left unspecified, will cause the corresponding functions to be applied over all sampling blocks. This allows testing of some or all of the samplers.
MCMC Types¶
The MCMC types and their relationships are depicted below with a Unified Modelling Language (UML) diagram. In the diagram, types are represented with boxes that display their respective names in the topmost panels, and fields in the second panels. By convention, plus signs denote fields that are publicly accessible, which is always the case for these structures in julia. Hollow triangle arrows point to types that the originator extends. Solid diamond arrows indicate that a number of instances of the type being pointed to are contained in the originator. The undirected line between Sampler
and Stochastic
represents a bidirectional association. Numbers on the graph indicate that there is one (1), zero or more (0..*), or one or more (1..*) instances of a type at the corresponding end of a relationship.
The relationships are as follows. Type Model
contains a dictionary field (Dict{Symbol, Any}
) of model nodes and a field (Vector{Sampler}
) of one or more sampling functions. Nodes can be one of three types:
 Stochastic nodes (
ScalarStochastic
orArrayStochastic
) are any model terms that have likelihood or prior distributional specifications. Logical nodes (
ScalarLogical
orArrayLogical
) are terms that are deterministic functions of other nodes. Input nodes (not shown) are any other model terms and data types that are considered to be fixed quantities in the analysis.
Stochastic
and Logical
are inherited from the Variate types and can be used with operators and in functions defined for that type. The sampling functions in Model
each correspond to a block of one or more model parameters (stochastic nodes) to be sampled from a target distribution (e.g. full conditional) during the simulation. Finally, ModelChains
stores simulation output for a given model. Detailed information about each type is provided in the subsequent sections.
Variate¶
ScalarVariate
and ArrayVariate{N}
are abstract types that serve as the basis for several concrete types in the Mamba package. Conceptually, they represent data structures that store numeric values simulated from target distributions. Being abstract, these variate types cannot be instantiated and cannot have fields. They can, however, have method functions, which descendant subtypes will inherit. Such inheritance allows one to endow a core set of functionality to all subtypes by simply defining method functions once on the abstract types (see julia Types). Accordingly, a core set of functionality is defined for the variate types through the field and method functions discussed below. Although the (abstract) types do not have fields, their method functions assume that all subtypes will be declared with a value
field.
Declarations¶
abstract ScalarVariate <: Real
abstract ArrayVariate{N} <: DenseArray{Float64, N}
typealias AbstractVariate Union{ScalarVariate, ArrayVariate}
typealias VectorVariate ArrayVariate{1}
typealias MatrixVariate ArrayVariate{2}
Type Hierarchy¶
Subtypes of the variate types include the Dependent, Logical, Stochastic, and SamplerVariate types.
Field¶
value::T
: scalar or array ofFloat64
values that represent simulated values from a target distribution.
Methods¶
Methods for ScalarVariate
and ArrayVariate
include mathematical operators, mathematical functions, and statistics defined in the base julia language for parent types Real
and DenseArray
. In addition, the following functions are provided.
Function  Description 

logit(x) 
logodds 
invlogit(x) 
inverse logodds 
Dependent¶
AbstractDependent
is an abstract type designed to store values and attributes of model nodes, including parameters to be simulated via MCMC, functions of the parameters, and likelihood specifications on observed data. It extends the base Variate
types with method functions defined for the fields summarized below. Like the type it extends, values are stored in a value
field and can be used with method functions that accept Float64
or Array{Float64, N}
type objects.
Since parameter values in the AbstractDependent
structure are stored as a scalar or array, objects of this type can be created for model parameters of corresponding dimensions, with the choice between the two being user and applicationspecific. At one end of the spectrum, a model might be formulated in terms of parameters that are all scalars, with a separate instances of AbstractDependent
for each one. At the other end, a formulation might be made in terms of a single parameter array, with one corresponding instance of AbstractDependent
. Whether to formulate parameters as scalars or arrays will depend on the application at hand. Array formulations should be considered for parameters and data that have multivariate distributions, or are to be used as such in numeric operations and functions. In other cases, scalar parametrizations may be preferable. Situations in which parameter arrays are often used include the specification of regression coefficients and random effects.
Declaration¶
typealias AbstractDependent Union{AbstractLogical, AbstractStochastic}
Fields¶
value::T
: scalar or array ofFloat64
values that represent samples from a target distribution.symbol::Symbol
: identifying symbol for the node.monitor::Vector{Int}
: indices identifying elements of thevalue
field to include in monitored MCMC sampler output.eval::Function
: function for updating the state of the node.sources::Vector{Symbol}
: other nodes upon whom the values of this one depends.targets::Vector{Symbol}
:Dependent
nodes that depend on this one. Elements oftargets
are topologically sorted so that a given node in the vector is conditionally independent of subsequent nodes, given the previous ones.
Display¶

show
(d::AbstractDependent)¶ Write a text representation of nodal values and attributes to the current output stream.

showall
(d::AbstractDependent)¶ Write a verbose text representation of nodal values and attributes to the current output stream.
Initialization¶

setmonitor!
(d::AbstractDependent, monitor::Bool)¶ 
setmonitor!
(d::AbstractDependent, monitor::Vector{Int}) Specify node elements to be included in monitored MCMC sampler output.
Arguments
d
: node whose elements contain sampled MCMC values.monitor
: boolean indicating whether all elements are monitored, or vector of elementwise indices of elements to monitor.
Value
Returnsd
with itsmonitor
field updated to reflect the specified monitoring.
Node Operations¶

logpdf
(d::AbstractDependent, transform::Bool=false)¶ 
logpdf
(d::AbstractDependent, x, transform::Bool=false) Evaluate the logdensity function for a node. In this method, no density function is assumed for the node, and a constant value of 0 is returned. This method function may be redefined for subtypes of
AbstractDependent
that have distributional specifications.Arguments
d
: node for which to evaluate the logdensity.x
: value, of the same type and shape as the node value, at which to perform the evaluation. If not specified, the node value is used.transform
: whether the evaluation is on the linktransformed scale.
Value
The resulting numeric value of the logdensity.

unlist
(d::AbstractDependent, transform::Bool=false)¶ 
unlist
(d::AbstractDependent, x::Real, transform::Bool=false) 
unlist
(d::AbstractDependent, x::AbstractArray, transform::Bool=false) 
relist
(d::AbstractDependent, x::AbstractArray, transform::Bool=false)¶ Extract (unlist) node values to a vector, or reassemble (relist) values to be put into a node. In this generic method, all values are listed. The methods are used internally for the extraction of unique stochastic node values to sample, and can be redefined to implement different behaviors for
AbstractDependent
subtypes.Arguments
d
: node for which to unlist or relist values.x
: values to be listed. If not specified, the node values are used.transform
: whether to apply a link or inverselink transformation to the values. In this generic method, transformations are defined to be the identity function.
Value
Returns unmodifiedx
values as a vector (unlist) or in the same shape as the specified node (relist).
Logical¶
The Logical
types inherit fields and method functions from the AbstractDependent
type, and adds the constructors and methods listed below. It is designed for nodes that are deterministic functions of model parameters and data.
Declarations¶
type ScalarLogical <: ScalarVariate
type ArrayLogical{N} <: ArrayVariate{N}
typealias AbstractLogical Union{ScalarLogical, ArrayLogical}
Fields¶
value
: values of typeFloat64
forScalarLogical
nodes andArray{Float64}
forArrayLogical
nodes that represent samples from a target distribution.symbol::Symbol
: identifying symbol for the node.monitor::Vector{Int}
: indices identifying elements of thevalue
field to include in monitored MCMC sampler output.eval::Function
: function for updating values stored invalue
.sources::Vector{Symbol}
: other nodes upon whom the values of this one depends.targets::Vector{Symbol}
:Dependent
nodes that depend on this one. Elements oftargets
are topologically sorted so that a given node in the vector is conditionally independent of subsequent nodes, given the previous ones.
Constructors¶

Logical
(f::Function, monitor::Union{Bool, Vector{Int}}=true)¶ 
Logical
(d::Integer, f::Function, monitor::Union{Bool, Vector{Int}}=true) Construct a
Logical
object that defines a logical model node.Arguments
d
: number of dimensions for array nodes.f
: function whose untyped arguments are the other model nodes upon which this one depends. The function may contain any valid julia expression or code block. It will be saved in theeval
field of the constructed logical node and should return a value in the same type as and with which to update the node’svalue
field.monitor
: boolean indicating whether all elements are monitored, or vector of elementwise indices of elements to monitor.
Value
Returns anArrayLogical
if the dimension argumentd
is specified, and aScalarLogical
if not.Example
See the Model Specification section of the tutorial.
Stochastic¶
The Stochastic
types inherit fields and method functions from the AbstractDependent
type, and adds the additional ones listed below. It is designed for model parameters or data that have distributional or likelihood specifications, respectively. Its stochastic relationship to other nodes and data structures is represented by the structure stored in distr
field.
Declarations¶
type ScalarStochastic <: ScalarVariate
type ArrayStochastic{N} <: ArrayVariate{N}
typealias AbstractStochastic Union{ScalarStochastic, ArrayStochastic}
Fields¶
value
: values of typeFloat64
forScalarStochastic
nodes andArray{Float64}
forArrayStochastic
nodes that represent samples from a target distribution.symbol::Symbol
: identifying symbol for the node.monitor::Vector{Int}
: indices identifying elements of thevalue
field to include in monitored MCMC sampler output.eval::Function
: function for updating thedistr
field for the node.sources::Vector{Symbol}
: other nodes upon whom the distributional specification for this one depends.targets::Vector{Symbol}
:Dependent
nodes that depend on this one. Elements oftargets
are topologically sorted so that a given node in the vector is conditionally independent of subsequent nodes, given the previous ones.distr
: distributional specification of typeUnivariateDistribution
forScalarStochastic
nodes andDistributionStruct
forArrayStochastic
nodes.
Distribution Structures¶
The DistributionStruct
alias defines the types of distribution structures supported for AbstractStochastic
nodes. Single Distribution
types from the Distributions section, arrays of UnivariateDistribution
, and arrays of MultivariateDistribution
objects are supported. When a MultivariateDistribution
array is specified for a stochastic node, the node is assumed to be one dimension bigger than the array, with the last dimension containing values from the distributions stored in the previous dimensions. Such arrays may contain distributions of different lengths. Model specification syntax for all three types of distribution structures can be seen in the Birats Example.
typealias DistributionStruct Union{Distribution,
Array{UnivariateDistribution},
Array{MultivariateDistribution}}
Constructors¶

Stochastic
(f::Function, monitor::Union{Bool, Vector{Int}}=true)¶ 
Stochastic
(d::Integer, f::Function, monitor::Union{Bool, Vector{Int}}=true) Construct a
Stochastic
object that defines a stochastic model node.Arguments
d
: number of dimensions for array nodes.f
: function whose untyped arguments are the other model nodes upon which this one depends. The function may contain any valid julia expression or code block. It will be saved in theeval
field of the constructed stochastic node and should return aDistributionStruct
object to be stored in the node’sdistr
field.monitor
: boolean indicating whether all elements are monitored, or vector of elementwise indices of elements to monitor.
Value
Returns anArrayStochastic
if the dimension argumentd
is specified, and aScalarStochastic
if not.Example
See the Model Specification section of the tutorial.
Initialization¶

setinits!
(s::Stochastic, m::Model, x=nothing) Set initial values for a stochastic node.
Arguments
s
: stochastic node to which to assign initial values.m
: model containing the node.x
: values to assign to the node.
Value
Returns the node with its assigned initial values.
Node Operations¶

logpdf
(s::AbstractStochastic, transform::Bool=false) 
logpdf
(s::AbstractStochastic, x, transform::Bool=false) Evaluate the logdensity function for a stochastic node.
Arguments
s
: stochastic node for which to evaluate the logdensity.x
: value, of the same type and shape as the node value, at which to perform the evaluation. If not specified, the node value is used.transform
: whether the evaluation is on the linktransformed scale.
Value
The resulting numeric value of the logdensity.

rand
(s::AbstractStochastic)¶ Draw a sample from the distributional specification on a stochastic node.
Arguments
s
: stochastic node from which to generate a random sample.
Value
Returns the sampled value(s).

unlist
(s::AbstractStochastic, transform::Bool=false) 
unlist
(s::AbstractStochastic, x::Real, transform::Bool=false) 
unlist
(s::AbstractStochastic, x::AbstractArray, transform::Bool=false) 
relist
(s::AbstractStochastic, x::AbstractArray, transform::Bool=false) Extract (unlist) stochastic node values to a vector, or reassemble (relist) values into a format that can be put into a node. These methods are used internally to extract the unique and sampled values of stochastic nodes. They are used, for instance, to extract only the unique, uppertriangular portions of (symmetric) covariance matrices and only the sampled values of
Array{MultivariateDistribution}
specifications whose distributions may be of different lengths.Arguments
s
: stochastic node for which to unlist or relist values.x
: values to be listed. If not specified, the node values are used.transform
: whether to apply a link transformation, or its inverse, to map values in a constrained distributional support to an unconstrained space. Supports for continuous, univariate distributions and positivedefinite matrix distributions (Wishart or inverseWishart) are transformed as follows: Lower and upper bounded: scaled and shifted to the unit interval and logittransformed.
 Lower bounded: shifted to zero and logtransformed.
 Upper bounded: scaled by 1, shifted to zero, and logtransformed.
 Positivedefinite matrix: compute the (uppertriangular) Cholesky decomposition, and return it with the diagonal elements logtransformed.
Value
Returns the extractedx
values as a vector or the reassembled values in the same shape as the specified node.

update!
(s::AbstractStochastic, m::Model) Update the values of a stochastic node according to its relationship with others in a model.
Arguments
s
: stochastic node to update.m
: model containing the node.
Value
Returns the node with its values updated.
Distributions¶
Given in this section are distributions, as provided by the Distributions [2] and Mamba packages, supported for the specification of Stochastic nodes. Truncated versions of continuous univariate distributions are also supported.
Univariate Distributions¶
Distributions Package Univariate Types¶
The following univariate types from the Distributions package are supported.
Arcsine DiscreteUniform InverseGamma NoncentralChisq SymTriangularDist
Bernoulli Edgeworth InverseGaussian NoncentralF TDist
Beta Epanechnikov Kolmogorov NoncentralHypergeometric TriangularDist
BetaPrime Erlang KSDist NoncentralT Triweight
Binomial Exponential KSOneSided Normal Uniform
Biweight FDist Laplace NormalCanon VonMises
Categorical Frechet Levy Pareto Weibull
Cauchy Gamma Logistic PoissonBinomial
Chi Geometric LogNormal Poisson
Chisq Gumbel NegativeBinomial Rayleigh
Cosine Hypergeometric NoncentralBeta Skellam
Flat Distribution¶
A Flat distribution is supplied with the degenerate probability density function:
Flat() # Flat distribution
UserDefined Univariate Distributions¶
New known, unknown, or unnormalized univariate distributions can be created and added to Mamba as subtypes of the Distributions package ContinuousUnivariateDistribution
or DiscreteUnivariateDistribution
types. Mamba requires only a partial implementation of the method functions described in the full instructions for creating univariate distributions. The specific workflow is given below.
Create a
quote
block for the new distribution. Assign the block a variable name, sayextensions
, preceded by the@everywhere
macro to ensure compatibility when julia is run in multiprocessor mode.The Distributions package contains types and method definitions for new distributions. Load the package and import any of its methods (indicated below) that are extended.
Declare the new distribution subtype, say
D
, within the block. Any constructors explicitly defined for the subtype should accept untyped or abstracttype (Real
,AbstractArray
, orDenseArray
) arguments. Implementing constructors in this way ensures that they will be callable with the MambaStochastic
andLogical
types.Extend/define the following Distributions package methods for the new distribution
D
.Test the subtype.
Add the
quote
block (new distribution) to Mamba with the following calls.using Mamba @everywhere eval(extensions)
Below is a univariate example based on the linear regression model in the Tutorial.
## Define a new univariate Distribution type for Mamba.
## The definition must be placed within an unevaluated quote block.
@everywhere extensions = quote
## Load needed packages and import methods to be extended
using Distributions
import Distributions: minimum, maximum, logpdf
## Type declaration
type NewUnivarDist <: ContinuousUnivariateDistribution
mu::Float64
sigma::Float64
end
## The following method functions must be implemented
## Minimum and maximum support values
minimum(d::NewUnivarDist) = Inf
maximum(d::NewUnivarDist) = Inf
## Normalized or unnormalized logdensity value
function logpdf(d::NewUnivarDist, x::Real)
log(d.sigma)  0.5 * ((x  d.mu) / d.sigma)^2
end
end
## Test the extensions in a temporary module (optional)
module Testing end
eval(Testing, extensions)
d = Testing.NewUnivarDist(0.0, 1.0)
Testing.minimum(d)
Testing.maximum(d)
Testing.insupport(d, 2.0)
Testing.logpdf(d, 2.0)
## Add the extensions
using Mamba
@everywhere eval(extensions)
## Implement a Mamba model using the new distribution
model = Model(
y = Stochastic(1,
(mu, s2) >
begin
sigma = sqrt(s2)
UnivariateDistribution[
NewUnivarDist(mu[i], sigma) for i in 1:length(mu)
]
end,
false
),
mu = Logical(1,
(xmat, beta) >
xmat * beta,
false
),
beta = Stochastic(1,
() > MvNormal(2, sqrt(1000))
),
s2 = Stochastic(
() > InverseGamma(0.001, 0.001)
)
)
## Sampling Scheme
scheme = [NUTS(:beta),
Slice(:s2, 3.0)]
## Sampling Scheme Assignment
setsamplers!(model, scheme)
## Data
line = Dict{Symbol, Any}(
:x => [1, 2, 3, 4, 5],
:y => [1, 3, 3, 3, 5]
)
line[:xmat] = [ones(5) line[:x]]
## Initial Values
inits = [
Dict{Symbol, Any}(
:y => line[:y],
:beta => rand(Normal(0, 1), 2),
:s2 => rand(Gamma(1, 1))
)
for i in 1:3
]
## MCMC Simulation
sim = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3)
describe(sim)
Multivariate Distributions¶
Distributions Package Multivariate Types¶
The following multivariate types from the Distributions package are supported.
Dirichlet MvNormal MvTDist Multinomial MvNormalCanon VonMisesFisher
BlockDiagonal Multivariate Normal Distribution¶
A BlockDiagonal Multivariate Normal distribution is supplied with the probability density function:
where
BDiagNormal(mu, C) # multivariate normal with mean vector mu and block
# diagonal covariance matrix Sigma such that
# length(mu) = dim(Sigma), and Sigma_1 = ... = Sigma_m = C
# for a matrix C or Sigma_1 = C[1], ..., Sigma_m = C[m]
# for a vector of matrices C.
UserDefined Multivariate Distributions¶
New known, unknown, or unnormalized multivariate distributions can be created and added to Mamba as subtypes of the Distributions package ContinuousMultivariateDistribution
or DiscreteMultivariateDistribution
types. Mamba requires only a partial implementation of the method functions described in the full instructions for creating multivariate distributions. The specific workflow is given below.
Create a
quote
block for the new distribution. Assign the block a variable name, sayextensions
, preceded by the@everywhere
macro to ensure compatibility when julia is run in multiprocessor mode.The Distributions package contains types and method definitions for new distributions. Load the package and import any of its methods (indicated below) that are extended.
Declare the new distribution subtype, say
D
, within the block. Any constructors explicitly defined for the subtype should accept untyped or abstracttype (Real
,AbstractArray
, orDenseArray
) arguments. Implementing constructors in this way ensures that they will be callable with the MambaStochastic
andLogical
types.Extend/define the following Distributions package methods for the new distribution
D
.Test the subtype.
Add the
quote
block (new distribution) to Mamba with the following calls.using Mamba @everywhere eval(extensions)
Below is a multivariate example based on the linear regression model in the Tutorial.
## Define a new multivariate Distribution type for Mamba.
## The definition must be placed within an unevaluated quote block.
@everywhere extensions = quote
## Load needed packages and import methods to be extended
using Distributions
import Distributions: length, insupport, _logpdf
## Type declaration
type NewMultivarDist <: ContinuousMultivariateDistribution
mu::Vector{Float64}
sigma::Float64
end
## The following method functions must be implemented
## Dimension of the distribution
length(d::NewMultivarDist) = length(d.mu)
## Logical indicating whether x is in the support
function insupport{T<:Real}(d::NewMultivarDist, x::AbstractVector{T})
length(d) == length(x) && all(isfinite(x))
end
## Normalized or unnormalized logdensity value
function _logpdf{T<:Real}(d::NewMultivarDist, x::AbstractVector{T})
length(x) * log(d.sigma)  0.5 * sumabs2(x  d.mu) / d.sigma^2
end
end
## Test the extensions in a temporary module (optional)
module Testing end
eval(Testing, extensions)
d = Testing.NewMultivarDist([0.0, 0.0], 1.0)
Testing.insupport(d, [2.0, 3.0])
Testing.logpdf(d, [2.0, 3.0])
## Add the extensions
using Mamba
@everywhere eval(extensions)
## Implement a Mamba model using the new distribution
model = Model(
y = Stochastic(1,
(mu, s2) > NewMultivarDist(mu, sqrt(s2)),
false
),
mu = Logical(1,
(xmat, beta) > xmat * beta,
false
),
beta = Stochastic(1,
() > MvNormal(2, sqrt(1000))
),
s2 = Stochastic(
() > InverseGamma(0.001, 0.001)
)
)
## Sampling Scheme
scheme = [NUTS(:beta),
Slice(:s2, 3.0)]
## Sampling Scheme Assignment
setsamplers!(model, scheme)
## Data
line = Dict{Symbol, Any}(
:x => [1, 2, 3, 4, 5],
:y => [1, 3, 3, 3, 5]
)
line[:xmat] = [ones(5) line[:x]]
## Initial Values
inits = [
Dict{Symbol, Any}(
:y => line[:y],
:beta => rand(Normal(0, 1), 2),
:s2 => rand(Gamma(1, 1))
)
for i in 1:3
]
## MCMC Simulation
sim = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3)
describe(sim)
MatrixVariate Distributions¶
Distributions Package MatrixVariate Types¶
The following matrixvariate types from the Distributions package are supported.
InverseWishart Wishart
Sampler¶
The Sampler
type stores modelbased Sampling Functions for use in the Mamba Gibbs sampling scheme. Developers can use it as a wrapper for calling standalone samplers or as a structure for implementing selfcontained samplers.
Declaration¶
type Sampler{T}
Fields¶
params::Vector{Symbol}
: stochastic nodes in the block being updated by the sampler.eval::Function
: sampling function that updates values of theparams
nodes.tune::T
: tuning parameters needed by the sampling function.targets::Vector{Symbol}
:Dependent
nodes that depend on and whose states must be updated afterparams
. Elements oftargets
are topologically sorted so that a given node in the vector is conditionally independent of subsequent nodes, given the previous ones.
Constructors¶

Sampler
(param::Symbol, f::Function, tune::Any=Dict())¶ 
Sampler
(params::Vector{Symbol}, f::Function, tune::Any=Dict()) Construct a
Sampler
object that defines a sampling function for a block of stochastic nodes.Arguments
param/params
: node(s) being blockupdated by the sampler.f
: function for theeval
field of the constructed sampler and whose arguments are the other model nodes upon which the sampler depends, typed argumentmodel::Model
that contains all model nodes, and/or typed argumentblock::Integer
that is an index identifying the corresponding sampling function in a vector of all samplers for the associated model. Through the arguments, all model nodes and fields can be accessed in the body of the function. The function may return an updated sample for the nodes identified in itsparams
field. Such a return value can be a structure of the same type as the node if the block consists of only one node, or a dictionary of node structures with keys equal to the block node symbols if one or more. Alternatively, a value ofnothing
may be returned. Return values that are notnothing
will be used to automatically update the node values and propagate them to dependent nodes. No automatic updating will be done ifnothing
is returned.tune
: tuning parameters needed by the sampling function.
Value
Returns aSampler{typeof(tune)}
type object.Example
See the Sampling Schemes section of the tutorial.
SamplerVariate¶
The SamplerVariate
type is designed to store simulated values from and tuning parameters for standalone Sampling Functions. It is a parametric type whose parameter can be any subtype of the abstract SamplerTune
type and serves to identify the family of sampling functions to which the variate belongs.
Declaration¶
abstract SamplerTune
type SamplerVariate{T<:SamplerTune} <: VectorVariate
Fields¶
value::Vector{Float64}
: simulated values.tune::T
: tuning parameters.
Constructors¶

SamplerVariate
(x::AbstractVector{U<:Real}, tune::SamplerTune)¶ 
SamplerVariate{T<:SamplerTune}
(x::AbstractVector{U<:Real}, tune::T)¶ 
SamplerVariate{T<:SamplerTune}
(x::AbstractVector{U<:Real}, pargs...; kargs...) Construct a
SamplerVariate
object for storing simulated values and tuning parameters.Arguments
x
: simulated values.tune
: tuning parameters. If not specified, the tuning parameter constructor is called with thevalue
field of the variate to instantiate the parameters.T
: explicit tuning parameter type for the variate. If not specified, the type is inferred from thetune
argument.pargs...
,kargs...
: variable positional and keyword arguments that are passed, along with thevalue
field of the variate, to the tuning parameter constructor asT(value, pargs...; kargs...)
. Accordingly, the arguments that this version of theSamplerVariate{T}
constructor accepts are defined by theT
constructors implemented for it.
Value
Returns a
SamplerVariate{T}
type object with fields set to the suppliedx
and tuning parameter values.
Model¶
The Model
type is designed to store the set of all model nodes, including parameter set as denoted in the Mamba Gibbs sampling scheme. In particular, it stores Dependent
type objects in its nodes
dictionary field. Valid models are ones whose nodes form directed acyclic graphs (DAGs). Sampling functions are saved as Sampler
objects in the vector of field samplers
. Vector elements correspond to sampling blocks
Declaration¶
type Model
Fields¶
nodes::Dict{Symbol, Any}
: all input, logical, and stochastic model nodes.samplers::Vector{Sampler}
: sampling functions for updating blocks of stochastic nodes.states::Vector{ModelState}
: states of chains at the end of an MCMC run in a possible series of runs, whereModelState
has fieldsvalue::Vector{Float64}
andtune::Vector{Any}
to store the last values of sampled nodes and blocksampler tuning parameters, respectively.iter::Int
: current MCMC draw from the target distribution.burnin::Int
: number of initial draws to discard as a burnin sequence to allow for convergence.hasinputs::Bool
: whether values have been assigned to input nodes.hasinits::Bool
: whether initial values have been assigned to stochastic nodes.
Constructor¶

Model
(; iter::Integer=0, burnin::Integer=0, samplers::Vector{Sampler}=Sampler[], nodes...)¶ Construct a
Model
object that defines a model for MCMC simulation.Arguments
iter
: current iteration of the MCMC simulation.burnin
: number of initial draws to be discarded as a burnin sequence to allow for convergence.samplers
: blockspecific sampling functions.nodes...
: arbitrary number of userspecified arguments defining logical and stochastic nodes in the model. Argument values must beLogical
orStochastic
type objects. Their names in the model will be taken from the argument names.
Value
Returns aModel
type object.Example
See the Model Specification section of the tutorial.
MCMC Engine¶

mcmc
(m::Model, inputs::Dict{Symbol}, inits::Vector{Dict{Symbol, Any}}, iters::Integer; burnin::Integer=0, thin::Integer=1, chains::Integer=1, verbose::Bool=true)¶ 
mcmc
(mc::ModelChains, iters::Integer; verbose::Bool=true) Simulate MCMC draws for a specified model.
Arguments
m
: specified model.mc
: chains from a previous call tomcmc
for which to simulate additional draws.inputs
: values for input model nodes. Dictionary keys and values should be given for each input node.inits
: dictionaries that contain initial values for stochastic model nodes. Dictionary keys and values should be given for each stochastic node. Consecutive runs of the simulator will iterate through the vector’s dictionary elements.iters
: number of draws to generate for each simulation run.burnin
: numer of initial draws to discard as a burnin sequence to allow for convergence.thin
: stepsize between draws to output.chains
: number of simulation runs to perform.verbose
: whether to print sampler progress at the console.
Value
AModelChains
type object of simulated draws.Example
See the MCMC Simulation section of the tutorial.
Indexing¶

getindex
(m::Model, nodekey::Symbol)¶ Returns a model node identified by its symbol. The syntax
m[nodekey]
is converted togetindex(m, nodekey)
.Arguments
m
: model containing the node to get.nodekey
: node to get.
Value
The specified node.

keys
(m::Model)¶ 
keys
(m::Model, ntype::Symbol, at...) Extract the symbols (keys) for all existing nodes or for nodes of a specified type.
Arguments
m
: model containing the nodes of interest.ntype
: type of nodes to return. Options are:all
: all input, logical, and stochastic model nodes.:assigned
: nodes that have been assigned values.:block
: stochastic nodes being updated by the sampling block(s)at::Integer=0
(default: all blocks).:dependent
: logical and stochastic (dependent) nodes in topologically sorted order.:independent
or:input
: input (independent) nodes.:logical
: logical nodes.:monitor
: stochastic nodes being monitored in MCMC sampler output.:output
: stochastic nodes upon which no other stochastic nodes depend.:source
: nodes upon which the nodeat::Symbol
or vector of nodesat::Vector{Symbol}
depends.:stochastic
: stochastic nodes.:target
: topologically sorted nodes that depend on the sampling block(s)at::Integer=0
(default: all blocks), nodeat::Symbol
, or vector of nodesat::Vector{Symbol}
.
at...
: additional positional arguments to be passed to thentype
options, as described above.
Value
A vector of node symbols.
Display¶

draw
(m::Model; filename::AbstractString="")¶ Draw a GraphViz DOTformatted graph representation of model nodes and their relationships.
Arguments
m
: model for which to construct a graph.filename
: external file to which to save the resulting graph, or an empty string to draw to standard output (default). If a supplied external file name does not include a dot (.
), the file extension.dot
will be appended automatically.
Value
The model drawn to an external file or standard output. Stochastic, logical, and input nodes will be represented by ellipses, diamonds, and rectangles, respectively. Nodes that are unmonitored in MCMC simulations will be graycolored.Example
See the Directed Acyclic Graphs section of the tutorial.

graph
(m::Model)¶ Construct a graph representation of model nodes and their relationships.
Arguments
m
: model for which to construct a graph.
Value
Returns aModelGraph
type object with fieldgraph
containing aDiGraph
representation of indices, as defined in the LightGraphs package, to a vector of node symbols in fieldkeys
.

graph2dot
(m::Model)¶ Draw a GraphViz DOTformatted graph representation of model nodes and their relationships.
Arguments
m
: model for which to construct a graph.
Value
A character string representation of the graph suitable for inline processing. Stochastic, logical, and input nodes will be represented by ellipses, diamonds, and rectangles, respectively. Nodes that are unmonitored in MCMC simulations will be graycolored.Example
See the Directed Acyclic Graphs section of the tutorial.

show
(m::Model)¶ Write a text representation of the model, nodes, and attributes to the current output stream.

showall
(m::Model)¶ Write a verbose text representation of the model, nodes, and attributes to the current output stream.
Initialization¶

setinits!
(m::Model, inits::Dict{Symbol, Any})¶ Set the initial values of stochastic model nodes.
Arguments
m
: model with nodes to be initialized.inits
: initial values for stochastic model nodes. Dictionary keys and values should be given for each stochastic node.
Value
Returns the model with stochastic nodes initialized and theiter
field set equal to 0.Example
See the Development and Testing section of the tutorial.

setinputs!
(m::Model, inputs::Dict{Symbol, Any})¶ Set the values of input model nodes.
Arguments
m
: model with input nodes to be assigned.inputs
: values for input model nodes. Dictionary keys and values should be given for each input node.
Value
Returns the model with values assigned to input nodes.Example
See the Development and Testing section of the tutorial.

setsamplers!
(m::Model, samplers::Vector{T<:Sampler})¶ Set the blocksamplers for stochastic model nodes.
Arguments
m
: model with stochastic nodes to be sampled.samplers
: blockspecific samplers.
Values:
Returns the model updated with the blocksamplers.Example
See the Model Specification and MCMC Simulation sections of the tutorial.
Parameter Block Operations¶

gettune
(m::Model, block::Integer=0)¶ Get blocksampler tuning parameters.
Arguments
m
: model with blocksamplers.block
: block for which to get the tuning parameters (default: all blocks).
Value
AVector{Any}
of all blockspecific tuning parameters ifblock=0
, and turning parameters for the specified block otherwise.

gradlogpdf
(m::Model, block::Integer=0, transform::Bool=false; dtype::Symbol=:forward)¶ 
gradlogpdf
(m::Model, x::AbstractVector{T<:Real}, block::Integer=0, transform::Bool=false; dtype::Symbol=:forward) 
gradlogpdf!
(m::Model, x::AbstractVector{T<:Real}, block::Integer=0, transform::Bool=false; dtype::Symbol=:forward)¶ Compute the gradient of logdensities for stochastic nodes.
Arguments
m
: model containing the stochastic nodes for which to compute the gradient.block
: sampling block of stochastic nodes for which to compute the gradient (default: all stochastic nodes).x
: value (possibly different than the current one) at which to compute the gradient.transform
: whether to compute the gradient of block parameters on the link–transformed scale.dtype
: type of differentiation for gradient calculations. Options are:central
: central differencing.:forward
: forward differencing.
Value
The resulting gradient vector. Methodgradlogpdf!()
additionally updates modelm
with supplied valuesx
.Note
Numerical approximation of derivatives by central and forward differencing is performed with the Calculus package [97].

logpdf
(m::Model, block::Integer=0, transform::Bool=false)¶ 
logpdf
(m::Model, nodekeys::Vector{Symbol}, transform::Bool=false) 
logpdf
(m::Model, x::AbstractArray{T<:Real}, block::Integer=0, transform::Bool=false) 
logpdf!
(m::Model, x::AbstractArray{T<:Real}, block::Integer=0, transform::Bool=false)¶ Compute the sum of logdensities for stochastic nodes.
Arguments
m
: model containing the stochastic nodes for which to evaluate logdensities.block
: sampling block of stochastic nodes over which to sum densities (default: all stochastic nodes).nodekeys
: nodes over which to sum densities.x
: value (possibly different than the current one) at which to evaluate densities.transform
: whether to evaluate evaluate logdensities of block parameters on the link–transformed scale.
Value
The resulting numeric value of summed logdensities. Methodlogpdf!()
additionally updates modelm
with supplied valuesx
.

sample!
(m::Model, block::Integer=0)¶ Generate one MCMC sample of values for a specified model.
Argument:
m
: model specification.block
: block for which to sample values (default: all blocks).
Value
Returns the model updated with the MCMC sample and, in the case ofblock=0
, theiter
field incremented by 1.Example
See the Development and Testing section of the tutorial.

unlist
(m::Model, block::Integer=0, transform::Bool=false)¶ 
unlist
(m::Model, nodekeys::Vector{Symbol}, transform::Bool=false) 
relist
(m::Model, x::AbstractArray{T<:Real}, block::Integer=0, transform::Bool=false)¶ 
relist
(m::Model, x::AbstractArray{T<:Real}, nodekeys::Vector{Symbol}, transform::Bool=false) 
relist!
(m::Model, x::AbstractArray{T<:Real}, block::Integer=0, transform::Bool=false)¶ 
relist!
(m::Model, x::AbstractArray{T<:Real}, nodekey::Symbol, transform::Bool=false) Convert (unlist) sets of logical and/or stochastic node values to vectors, or reverse (relist) the process.
Arguments
m
: model containing nodes to be unlisted or relisted.block
: sampling block of nodes to be listed (default: all blocks).nodekey/nodekeys
: node(s) to be listed.x
: values to relist.transform
: whether to apply a link transformation in the conversion.
Value
Theunlist
methods return vectors of concatenated node values,relist
return dictionaries of symbol keys and values for the specified nodes, andrelist!
return their model argument with values copied to the nodes.

update!
(m::Model, block::Integer=0)¶ 
update!
(m::Model, nodekeys::Vector{Symbol}) Update values of logical and stochastic model node according to their relationship with others in a model.
Arguments
m
: mode with nodes to be updated.block
: sampling block of nodes to be updated (default: all blocks).nodekeys
: nodes to be updated in the given order.
Value
Returns the model with updated nodes.
Chains¶
AbstractChains
subtypes store output from one or more runs (chains) of an MCMC sampler. They serve as containers for output generated by the mcmc()
function, and supply methods for convergence diagnostics and posterior inference. Moreover, they can be used as standalone containers for any usergenerated MCMC output, and are thus a julia analogue to the boa [86][87] and coda [72][73] R packages.
Declarations¶
abstract AbstractChains
immutable Chains <: AbstractChains
immutable ModelChains <: AbstractChains
Fields¶
value::Array{Float64, 3}
: 3dimensional array of sampled values whose first, second, and third dimensions index the iterations, parameter elements, and runs of an MCMC sampler, respectively.range::Range{Int}
: range of iterations stored in the rows of thevalue
array.names::Vector{AbstractString}
: names assigned to the parameter elements.chains::Vector{Int}
: indices to the MCMC runs.model::Model
: model from which the sampled values were generated (ModelChains
only).
Constructors¶

Chains
(iters::Integer, params::Integer; start::Integer=1, thin::Integer=1, chains::Integer=1, names::Vector{T<:AbstractString}=AbstractString[])¶ 
Chains
(value::Array{T<:Real, 3}; start::Integer=1, thin::Integer=1, names::Vector{U<:AbstractString}=AbstractString[], chains::Vector{V<:Integer}=Int[]) 
Chains
(value::Matrix{T<:Real}; start::Integer=1, thin::Integer=1, names::Vector{U<:AbstractString}=AbstractString[], chains::Integer=1) 
Chains
(value::Vector{T<:Real}; start::Integer=1, thin::Integer=1, names::AbstractString="Param1", chains::Integer=1) 
ModelChains
(c::Chains, m::Model)¶ Construct a
Chains
orModelChains
object that stores MCMC sampler output.Arguments
iters
: total number of iterations in each sampler run, of whichlength(start:thin:iters)
outputted iterations will be stored in the object.params
: number of parameters to store.value
: array whose first, second (optional), and third (optional) dimensions index outputted iterations, parameter elements, and runs of an MCMC sampler, respectively.start
: number of the first iteration to be stored.thin
: number of steps between consecutive iterations to be stored.chains
: number of simulation runs for which to store output, or indices to the runs (default: 1, 2, …).names
: names to assign to the parameter elements (default:"Param1"
,"Param2"
, …).m
: model for which simulated values were generated.
Value
Returns an object of typeChains
orModelChains
according to the name of the constructor called.Example
Indexing and Concatenation¶

cat
(dim::Integer, chains::AbstractChains...)¶ 
vcat
(chains::AbstractChains...)¶ 
hcat
(chains::AbstractChains...)¶ Concatenate input MCMC chains along a specified dimension. For dimensions other than the specified one, all input chains must have the same sizes, which will also be the sizes of the output chain. The size of the output chain along the specified dimension will be the sum of the sizes of the input chains in that dimension.
vcat
concatenates vertically along dimension 1, and has the alternative syntax[chain1; chain2; ...]
.hcat
concatenates horizontally along dimension 2, and has the alternative syntax[chain1 chain2 ...]
.Arguments
dim
: dimension (1, 2, or 3) along which to concatenate the input chains.chains
: chains to concatenate.
Value
AChains
object containing the concatenated input.Example
See thereadcoda()
example.

first
(c::AbstractChains)¶ 
step
(c::AbstractChains)¶ 
last
(c::AbstractChains)¶ Get the first iteration, stepsize (thinning), or last iteration of MCMC sampler output.
Arguments
c
: sampler output for which to return results.
Value
Integer value of the requested iteration type.

getindex
(c::Chains, window, names, chains)¶ 
getindex
(mc::ModelChains, window, names, chains) Subset MCMC sampler output. The syntax
c[i, j, k]
is converted togetindex(c, i, j, k)
.Arguments
c
: sampler output to subset.window
: indices of the formstart:stop
orstart:thin:stop
can be used to subset iterations, wherestart
andstop
define a range for the subset andthin
will apply additional thinning to existing sampler output.names
: indices for subsetting of parameters that can be specified as strings, integers, or booleans identifying parameters to be kept.ModelChains
may additionally be indexed by model node symbols.chains
: indices for chains can be integers or booleans.
A value of
:
can be specified for any of the dimensions to indicate no subsetting.Value
Subsetted sampler output stored in the same type of object as that supplied in the call.Example
See the Output Subsetting section of the tutorial.

setindex!
(c::AbstractChains, value, iters, names, chains)¶ Store MCMC sampler output at a given index. The syntax
c[i, j, k] = value
is converted tosetindex!(c, value, i, j, k)
.Arguments
c
: object within which to store sampler output.value
: sampler output.iters
: iterations can be indexed as astart:stop
orstart:thin:stop
range, a single numeric index, or a vector of indices; and are taken to be relative to the index range store in thec.range
field.names
: indices for subsetting of parameters can be specified as strings, integers, or booleans.chains
: indices for chains can be integers or booleans.
A value of
:
can be specified for any of the dimensions to index all corresponding elements.Value
An object of the same type asc
with the sampler output stored in the specified indices.Example
File I/O¶

read
(name::AbstractString, ::Type{T<:AbstractChains})¶ 
write
(name::AbstractString, c::AbstractChains)¶ Read a chain from or write one to an external file.
Arguments
name
: file to read or write. Recommended convention is for the file name to be specified with a.jls
extension.T
: chain type to read.c
: chain to write.
Value
AnAbstractChains
subtype read from an external file, or a written external file containing a subtype.Example
See the File I/O section of the tutorial.

readcoda
(output::AbstractString, index::AbstractString)¶ Read MCMC sampler output generated in the CODA format by OpenBUGS [90]. The function only retains those sampler iterations at which all model parameters were monitored.
Arguments
output
: text file containing the iteration numbers and sampled values for the model parameters.index
: text file containing the names of the parameters, followed by the first and last rows in which their output can be found in theoutput
file.
Value
AChains
object containing the read sampler output.Example
The following example reads sampler output contained in the CODA files
line1.out
,line1.ind
,line2.out
, andline2.ind
.using Mamba ## Get the directory in which the example CODA files are saved dir = dirname(@__FILE__) ## Read the MCMC sampler output from the files c1 = readcoda(joinpath(dir, "line1.out"), joinpath(dir, "line1.ind")) c2 = readcoda(joinpath(dir, "line2.out"), joinpath(dir, "line2.ind")) ## Concatenate the resulting chains c = cat(3, c1, c2) ## Compute summary statistics describe(c)
Convergence Diagnostics¶
MCMC simulation provides autocorrelated samples from a target distribution. Because of computational complexities in implementing MCMC algorithms, the autocorrelated nature of samples, and the need to choose initial sampling values at different points in target distributions; it is important to evaluate the quality of resulting output. Specifically, one should check that MCMC samples have converged to the target (or, more commonly, are stationary) and that the number of convergent samples provides sufficiently accurate and precise estimates of posterior statistics.
Several established convergence diagnostics are supplied by Mamba. The diagnostics and their features are summarized in the table below and described in detail in the subsequent function descriptions. They differ with respect to the posterior statistic being assessed (mean vs. quantile), whether the application is to parameters univariately or multivariately, and the number of chains required for calculations. Diagnostics may assess stationarity, estimation accuracy and precision, or both. A more comprehensive comparative review can be found in [18]. Since diagnostics differ in their focus and design, it is often good practice to employ more than one to assess convergence. Note too that diagnostics generally test for nonconvergence and that nonsignificant test results do not prove convergence. Thus, nonsignificant results should be interpreted with care.
Convergence Assessments  

Diagnostic  Statistic  Parameters  Chains  Stationarity  Estimation 
Gelman, Rubin, and Brooks  Mean  Univariate  2+  Yes  No 
Multivariate  2+  Yes  No  
Geweke  Mean  Univariate  1  Yes  No 
Heidelberger and Welch  Mean  Univariate  1  Yes  Yes 
Raftery and Lewis  Quantile  Univariate  1  Yes  Yes 
Gelman, Rubin, and Brooks Diagnostics¶

gelmandiag
(c::AbstractChains; alpha::Real=0.05, mpsrf::Bool=false, transform::Bool=false)¶ Compute the convergence diagnostics of Gelman, Rubin, and Brooks [33][12] for MCMC sampler output. The diagnostics are designed to asses convergence of posterior means estimated with multiple autocorrelated samples (chains). They does so by comparing the between and withinchain variances with metrics called potential scale reduction factors (PSRF). Both univariate and multivariate factors are available to assess the convergence of parameters individually and jointly. Scale factors close to one are indicative of convergence. As a rule of thumb, convergence is concluded if the 0.975 quantile of an estimated factor is less than 1.2. Multiple chains are required for calculations. It is recommended that at least three chains be generated, each with different starting values chosen to be diffuse with respect to the anticipated posterior distribution. Use of multiple chains in the diagnostic provides for more robust assessment of convergence than is possible with single chain diagnostics.
Arguments
c
: sampler output on which to perform calculations.alpha
: quantile (1  alpha / 2
) at which to estimate the upper limits of scale reduction factors.mpsrf
: whether to compute the multivariate potential scale reduction factor. This factor will not be calculable if any one of the parameters in the output is a linear combination of others.transform
: whether to apply log or logit transformations, as appropriate, to parameters in the chain to potentially produce output that is more normally distributed, an assumption of the PSRF formulations.
Value
A
ChainSummary
type object of the form:immutable ChainSummary value::Array{Float64, 3} rownames::Vector{AbstractString} colnames::Vector{AbstractString} header::AbstractString end
with parameters contained in the rows of the
value
field, and scale reduction factors and upperlimit quantiles in the first and second columns.Example
See the Convergence Diagnostics section of the tutorial.
Geweke Diagnostic¶

gewekediag
(c::AbstractChains; first::Real=0.1, last::Real=0.5, etype=:imse, args...)¶ Compute the convergence diagnostic of Geweke [37] for MCMC sampler output. The diagnostic is designed to asses convergence of posterior means estimated with autocorrelated samples. It computes a normalbased test statistic comparing the sample means in two windows containing proportions of the first and last iterations. Users should ensure that there is sufficient separation between the two windows to assume that their samples are independent. A nonsignificant test pvalue indicates convergence. Significant pvalues indicate nonconvergence and the possible need to discard initial samples as a burnin sequence or to simulate additional samples.
Arguments
c
: sampler output on which to perform calculations.first
: proportion of iterations to include in the first window.last
: proportion of iterations to include in the last window.etype
: method for computing Monte Carlo standard errors. Seemcse()
for options.args...
: additional arguments to be passed to theetype
method.
Value
AChainSummary
type object with parameters contained in the rows of thevalue
field, and test Zscores and pvalues in the first and second columns. Results are chainspecific.Example
See the Convergence Diagnostics section of the tutorial.
Heidelberger and Welch Diagnostic¶

heideldiag
(c::AbstractChains; alpha::Real=0.05, eps::Real=0.1, etype=:imse, args...)¶ Compute the convergence diagnostic of Heidelberger and Welch [47] for MCMC sampler output. The diagnostic is designed to assess convergence of posterior means estimated with autocorrelated samples and to determine whether a target degree of accuracy is achieved. A stationarity test is performed for convergence assessment by iteratively discarding 10% of the initial samples until the test pvalue is nonsignificant and stationarity is concluded or until 50% have been discarded and stationarity is rejected, whichever occurs first. Then, a halfwidth test is performed by calculating the relative halfwidth of a posterior mean estimation interval as ; where is a standard normal quantile, is the Monte Carlo standard error, and is the estimated posterior mean. If the relative halfwidth is greater than a target ratio, the test is rejected. Rejection of the stationarity or halfwidth test suggests that additional samples are needed.
Arguments
c
: sampler output on which to perform calculations.alpha
: significance level for evaluations of stationarity tests and calculations of relative estimation interval halfwidths.eps
: target ratio for the relative halfwidths.etype
: method for computing Monte Carlo standard errors. Seemcse()
for options.args...
: additional arguments to be passed to theetype
method.
Value
AChainSummary
type object with parameters contained in the rows of thevalue
field, and numbers of burnin sequences to discard, whether the stationarity tests are passed (1 = yes, 0 = no), their pvalues ( implies stationarity), posterior means, halfwidths of their estimation intervals, and whether the halfwidth tests are passed (1 = yes, 0 = no) in the columns. Results are chainspecific.Example
See the Convergence Diagnostics section of the tutorial.
Raftery and Lewis Diagnostic¶

rafterydiag
(c::AbstractChains; q::Real=0.025, r::Real=0.005, s::Real=0.95, eps::Real=0.001)¶ Compute the convergence diagnostic of Raftery and Lewis [74][75] for MCMC sampler output. The diagnostic is designed to determine the number of autocorrelated samples required to estimate a specified quantile , such that , within a desired degree of accuracy. In particular, if is the estimand and the estimated cumulative probability, then accuracy is specified in terms of and , where . Thinning may be employed in the calculation of the diagnostic to satisfy its underlying assumptions. However, users may not want to apply the same (or any) thinning when estimating posterior summary statistics because doing so results in a loss of information. Accordingly, sample sizes estimated by the diagnostic tend to be conservative (too large).
Arguments
c
: sampler output on which to perform calculations.q
: posterior quantile of interest.r
: margin of error for estimated cumulative probabilities.s
: probability for the margin of error.eps
: tolerance within which the probabilities of transitioning from initial to retained iterations are within the equilibrium probabilities for the chain. This argument determines the number of samples to discard as a burnin sequence and is typically left at its default value.
Value
AChainSummary
type object with parameters contained in the rows of thevalue
field, and thinning intervals employed, numbers of samples to discard as burnin sequences, total numbers () to burnin and retain, numbers of independent samples that would be needed (), and dependence factors () in the columns. Results are chainspecific.Example
See the Convergence Diagnostics section of the tutorial.
Posterior Summary Statistics¶

autocor
(c::AbstractChains; lags::Vector=[1, 5, 10, 50], relative::Bool=true)¶ Compute lagk autocorrelations for MCMC sampler output.
Arguments
c
: sampler output on which to perform calculations.lags
: lags at which to compute autocorrelations.relative
: whether the lags are relative to the thinning interval of the output (true
) or relative to the absolute iteration numbers (false
).
Value
AChainSummary
type object with model parameters indexed by the first dimension ofvalue
, lagautocorrelations by the second, and chains by the third.Example
See the Posterior Summaries section of the tutorial.

changerate
(c::AbstractChains)¶ Estimate the probability, or rate per iteration, of a state space change for iterations in MCMC sampler output. Estimation is performed for each parameter univariately as well as for the full parameter vector multivariately. For continuous output generated from samplers, like MetropolisHastings, whose algorithms conditionally accept candidate draws, the probability can be viewed as the acceptance rate.
Arguments
c
: sampler output on which to perform calculations.
Value
AChainSummary
type object with parameters in the rows of thevalue
field, and the estimated rates in the column. Results are for all chains combined.Example
See the Posterior Summaries section of the tutorial.

cor
(c::AbstractChains)¶ Compute crosscorrelations for MCMC sampler output.
Arguments
c
: sampler output on which to perform calculations.
Value
AChainSummary
type object with the first and second dimensions of thevalue
field indexing the model parameters between which correlations. Results are for all chains combined.Example
See the Posterior Summaries section of the tutorial.

describe
(c::AbstractChains; q::Vector=[0.025, 0.25, 0.5, 0.75, 0.975], etype=:bm, args...)¶ Compute summary statistics for MCMC sampler output.
Arguments
c
: sampler output on which to perform calculations.q
: probabilities at which to calculate quantiles.etype
: method for computing Monte Carlo standard errors. Seemcse()
for options.args...
: additional arguments to be passed to theetype
method.
Value
Results from calls tosummarystats(c, etype, args...)
andquantile(c, q)
are printed for all chains combined, and a value ofnothing
is returned.Example
See the Posterior Summaries section of the tutorial.

hpd
(c::AbstractChains; alpha::Real=0.05)¶ Compute highest posterior density (HPD) intervals of Chen and Shao [16] for MCMC sampler output. HPD intervals have the desirable property of being the smallest intervals that contain a given probability. However, their calculation assumes unimodal marginal posterior distributions, and they are not invariant to transformations of parameters like central (quantilebased) posterior intervals.
Arguments
c
: sampler output on which to perform calculations.alpha
: the100 * (1  alpha)
% interval to compute.
Value
AChainSummary
type object with parameters contained in the rows of thevalue
field, and lower and upper intervals in the first and second columns. Results are for all chains combined.Example
See the Posterior Summaries section of the tutorial.

mcse
(x::Vector{T<:Real}, method::Symbol=:imse; args...)¶ Compute Monte Carlo standard errors.
Arguments
x
: time series of values on which to perform calculations.method
: method used for the calculations. Options are:bm
: batch means [41], with optional argumentsize::Integer=100
determining the number of sequential values to include in each batch. This method requires that the number of values inx
is at least 2 times the batch size.:imse
: initial monotone sequence estimator [38].:ipse
: initial positive sequence estimator [38].
args...
: additional arguments for the calculation method.
Value
The numeric standard error value.

quantile
(c::AbstractChains; q::Vector=[0.025, 0.25, 0.5, 0.75, 0.975])¶ Compute posterior quantiles for MCMC sampler output.
Arguments
c
: sampler output on which to perform calculations.q
: probabilities at which to compute quantiles.
Value
AChainSummary
type object with parameters contained in the rows of thevalue
field, and quantiles in the columns. Results are for all chains combined.

summarystats
(c::AbstractChains; etype=:bm, args...)¶ Compute posterior summary statistics for MCMC sampler output.
Arguments
c
: sampler output on which to perform calculations.etype
: method for computing Monte Carlo standard errors. Seemcse()
for options.args...
: additional arguments to be passed to theetype
method.
Value
AChainSummary
type object with parameters in the rows of thevalue
field; and the sample mean, standard deviation, standard error, Monte Carlo standard error, and effective sample size in the columns. Results are for all chains combined.
ModelBased Inference¶

dic
(mc::ModelChains)¶ Compute the Deviance Information Criterion (DIC) of Spiegelhalter et al. [91] and Gelman et al. [31] from MCMC sampler output.
Arguments
mc
: sampler output from a model fit with themcmc()
function.
Value
A
ChainSummary
type object with DIC results from the methods of Spiegelhalter and Gelman in the first and second rows of thevalue
field, and the DIC value and effective numbers of parameters in the first and second columns; wheresuch that is the loglikelihood of model outputs given the expected values of model parameters , and is the effective number of parameters. The latter is defined as for the method of Spiegelhalter and as for the method of Gelman. Results are for all chains combined.
Example
See the Posterior Summaries section of the tutorial.

logpdf
(mc::ModelChains, nodekey::Symbol)¶ 
logpdf
(mc::ModelChains, nodekeys::Vector{Symbol}=keys(mc.model, :stochastic)) Compute the sum of logdensities at each iteration of MCMC output for stochastic nodes.
Arguments
mc
: sampler output from a model fit with themcmc()
function.nodekey/nodekeys
: stochastic model node(s) over which to sum densities (default: all).
Value
AModelChains
object of resulting summed logdensities at each MCMC iteration of the supplied chain.

predict
(mc::ModelChains, nodekey::Symbol)¶ 
predict
(mc::ModelChains, nodekeys::Vector{Symbol}=keys(mc.model, :output)) Generate MCMC draws from a posterior predictive distribution.
Arguments
mc
: sampler output from a model fit with themcmc()
function.nodekey/nodekeys
: observed Stochastic model node(s) for which to generate draws from the predictive distribution (default: all observed data nodes).
Value
A
ModelChains
object of draws simulated at each MCMC iteration of the supplied chain. For observed data node , simulation is from the posterior predictive distributionwhere is an unknown observation on the node, is the data likelihood, and is the posterior distribution of unobserved parameters .
Example
See the Pumps example.
Plotting¶

plot
(c::AbstractChains, ptype::Vector{Symbol}=[:trace, :density]; legend::Bool=false, args...)¶ 
plot
(c::AbstractChains, ptype::Symbol; legend::Bool=false, args...) Various plots to summarize sampler output stored in
AbstractChains
subtypes. Separate plots are produced for each sampled parameter.Arguments
c
: sampler output to plot.ptype
: plot type(s). Options are:autocor
: autocorrelation plots, with optional argumentmaxlag::Integer=round(Int, 10 * log10(length(c.range)))
determining the maximum autocorrelation lag to plot. Lags are plotted relative to the thinning interval of the output.:bar
: bar plots. Optional argumentposition::Symbol=:stack
controls whether bars should be stacked on top of each other (default) or placed side by side (:dodge
).:contour
: pairwise posterior density contour plots. Optional argumentbins::Integer=100
controls the plot resolution.:density
: density plots. Optional argumenttrim::Tuple{Real, Real}=(0.025, 0.975)
trims off lower and upper quantiles of density.:mean
: running mean plots.:mixeddensity
: bar plots (:bar
) for parameters with integer values within bounds defined by optional argumentbarbounds::Tuple{Real, Real}=(0, Inf)
, and density plots (:density
) otherwise.:trace
: trace plots.
legend
: whether to include legends in the plots to identify chainspecific results.args...
: additional keyword arguments to be passed to theptype
options, as described above.
Value
Returns aVector{Plot}
whose elements are individual parameter plots of the specified type ifptype
is a symbol, and aMatrix{Plot}
with plot types in the rows and parameters in the columns ifptype
is a vector. The result can be displayed or saved to a file withdraw()
.Note
Plots are created using the Gadfly package [51].Example
See the Plotting section of the tutorial.

draw
(p::Array{Plot}; fmt::Symbol=:svg, filename::AbstractString="", width::MeasureOrNumber=8inch, height::MeasureOrNumber=8inch, nrow::Integer=3, ncol::Integer=2, byrow::Bool=true, ask::Bool=true)¶ Draw plots produced by
plot()
into display grids containing a default of 3 rows and 2 columns of plots.Arguments
p
: plots to be drawn. Elements ofp
are read in the order stored by julia (e.g. columnmajor order for matrices) and written to the display grid according to thebyrow
argument. Grids will be filled sequentially until all plots have been drawn.fmt
: output format. Options are:pdf
: Portable Document Format (.pdf).:pgf
: Portable Graphics Format (.pgf).:png
: Portable Network Graphics (.png).:ps
: Postscript (.ps).:svg
: Scalable Vector Graphics (.svg).
filename
: file to which to save the display grids as they are drawn, or an empty string to draw to the display device (default). If a supplied external file name does not include a dot (.
), then a hyphen followed by the grid sequence number and then the format extension will be appended automatically. In the case of multiple grids, the former file name behavior will write all grids to the single named file, but prompt users before advancing to the next grid and overwriting the file; the latter behavior will write each grid to a different file.width/height
: grid widths/heights incm
,mm
,inch
,pt
, orpx
units.nrow/ncol
: number of rows/columns in the display grids.byrow
: whether the display grids should be filled by row.ask
: whether to prompt users before displaying subsequent grids to a single named file or the display device.
Value
Grids drawn to an external file or the display device.Example
See the Plotting section of the tutorial.
Sampling Functions¶
Listed below are the sampling methods for which functions are provided to simulating draws from distributions that can be specified up to constants of proportionalities. Modelbased Sampler constructors are available for use with the mcmc()
engine as well as standalone functions that can be used independently.
Approximate Bayesian Computation (ABC)¶
Approximate Bayesian Computation in the framework of MCMC (also known as LikelihoodFree MCMC) as proposed by [59] for simulating autocorrelated draws from a posterior distribution without evaluating its likelihood. Also see [85] for a thorough review of LikelihoodFree MCMC.
ModelBased Constructor¶

ABC
(params::ElementOrVector{Symbol}, scale::ElementOrVector{T<:Real}, summary::Function, epsilon::Real; kernel::KernelDensityType=SymUniform, dist::Function=(Tsim, Tobs) → sqrt(sumabs2(Tsim  Tobs)), proposal::SymDistributionType=Normal, maxdraw::Integer=1, nsim::Integer=1, decay::Real=1.0, randeps::Bool=false, args...)¶ Construct a
Sampler
object for ABC sampling. Parameters are assumed to be continuous, but may be constrained or unconstrained.Arguments
params
: stochastic node(s) to be updated with the sampler. Constrained parameters are mapped to unconstrained space according to transformations defined by the Stochasticunlist()
function.scale
: scaling value or vector of the same length as the combined elements of nodesparams
for theproposal
distribution. Values are relative to the unconstrained parameter space, where candidate draws are generated.summary
: function that takes a vector of observed or simulated data and returns a summary statistic or vector of statistics.epsilon
: target tolerance for determining how similar observed and simulated data summary statistics need to be in order to accept a candidate draw.kernel
: weighting kernel density of typeBiweight
,Cosine
,Epanechnikov
,Normal
,SymTriangularDist
,SymUniform
, orTriweight
to use in measuring similarity between observed and simulated data summary statistics. Specifiedepsilon
determines the standard deviation of Normal kernels and widths of the others.dist
: positive function for the kernel density to compute distance between vectors of observed (Tobs
) and simulated (Tsim
) data summary statistics (default: Euclidean distance).proposal
: symmetric distribution of typeBiweight
,Cosine
,Epanechnikov
,Normal
,SymTriangularDist
,SymUniform
, orTriweight
to be centered around current parameter values and used to generate proposal draws. Specifiedscale
determines the standard deviations of Normal proposals and widths of the others.maxdraw
: maximum number of unaccepted candidates to draw in each call of the sampler. Draws are generated until one is accepted or the maximum is reached. Larger values increase acceptance rates at the expense of longer runtimes.nsim
: number of data sets to simulate in deciding whether to accept a candidate draw. Larger values lead to closer approximations of the target distribution at the expense of longer runtimes.decay
: if0 < decay <= 1
, the rate at which internal tolerances are monotonically decreased from the initial distance between observed and simulated summary statistics toward the maximum of each subsequent distance andepsilon
; ifdecay = 0
, internal tolerances are fixed atepsilon
.randeps
: whether to perturb internal tolerances by random exponential variates.args...
: additional keyword arguments to be passed to thedist
function.
Value
Returns aSampler{ABCTune}
type object.Example
ABCTune Type¶
Declaration¶
type ABCTune
Fields¶
datakeys::Vector{Symbol}
: stochastic “data” nodes in the full conditional distribution for parameters to be updated and nodes at which summary statistics are computed separately in the sampling algorithm.Tsim::Vector{Vector{Float64}}
: simulated data summary statistics for thensim
data sets.epsilon::Vector{Float64}
: internal tolerances for the data sets.epsilonprime::Vector{Float64}
: perturbed tolerances ifrandeps=true
orepsilon
otherwise.
Adaptive Mixture Metropolis (AMM)¶
Implementation of the Roberts and Rosenthal [79] adaptive (multivariate) mixture Metropolis [45][46][63] sampler for simulating autocorrelated draws from a distribution that can be specified up to a constant of proportionality.
ModelBased Constructor¶

AMM
(params::ElementOrVector{Symbol}, Sigma::Matrix{T<:Real}; adapt::Symbol=:all, args...)¶ Construct a
Sampler
object for AMM sampling. Parameters are assumed to be continuous, but may be constrained or unconstrained.Arguments
params
: stochastic node(s) to be updated with the sampler. Constrained parameters are mapped to unconstrained space according to transformations defined by the Stochasticunlist()
function.Sigma
: covariance matrix for the nonadaptive multivariate normal proposal distribution. The covariance matrix is relative to the unconstrained parameter space, where candidate draws are generated.adapt
: type of adaptation phase. Options are:all
: adapt proposal during all iterations.:burnin
: adapt proposal during burnin iterations.:none
: no adaptation (multivariate Metropolis sampling with fixed proposal).
args...
: additional keyword arguments to be passed to theAMMVariate
constructor.
Value
Returns aSampler{AMMTune}
type object.Example
StandAlone Function¶

sample!
(v::AMMVariate; adapt::Bool=true)¶ Draw one sample from a target distribution using the AMM sampler. Parameters are assumed to be continuous and unconstrained.
Arguments
v
: current state of parameters to be simulated. When running the sampler in adaptive mode, thev
argument in a successive call to the function will contain thetune
field returned by the previous call.adapt
: whether to adaptively update the proposal distribution.
Value
Returnsv
updated with simulated values and associated tuning parameters.Example
The following example samples parameters in a simple linear regression model. Details of the model specification and posterior distribution can be found in the Supplement.
################################################################################ ## Linear Regression ## y ~ N(b0 + b1 * x, s2) ## b0, b1 ~ N(0, 1000) ## s2 ~ invgamma(0.001, 0.001) ################################################################################ using Mamba ## Data data = Dict( :x => [1, 2, 3, 4, 5], :y => [1, 3, 3, 3, 5] ) ## Logtransformed Posterior(b0, b1, log(s2)) + Constant logf = function(x::DenseVector) b0 = x[1] b1 = x[2] logs2 = x[3] r = data[:y]  b0  b1 * data[:x] (0.5 * length(data[:y])  0.001) * logs2  (0.5 * dot(r, r) + 0.001) / exp(logs2)  0.5 * b0^2 / 1000  0.5 * b1^2 / 1000 end ## MCMC Simulation with Adaptive Multivariate Metopolis Sampling n = 5000 burnin = 1000 sim = Chains(n, 3, names = ["b0", "b1", "s2"]) theta = AMMVariate([0.0, 0.0, 0.0], eye(3), logf) for i in 1:n sample!(theta, adapt = (i <= burnin)) sim[i, :, 1] = [theta[1:2]; exp(theta[3])] end describe(sim)
AMMVariate Type¶
Declaration¶
typealias AMMVariate SamplerVariate{AMMTune}
Fields¶
value::Vector{Float64}
: simulated values.tune::AMMTune
: tuning parameters for the sampling algorithm.
Constructor¶

AMMVariate
(x::AbstractVector{T<:Real}, Sigma::Matrix{U<:Real}, logf::Function; beta::Real=0.05, scale::Real=2.38)¶ Construct an
AMMVariate
object that stores simulated values and tuning parameters for AMM sampling.Arguments
x
: initial values.Sigma
: covariance matrix for the nonadaptive multivariate normal proposal distribution. The covariance matrix is relative to the unconstrained parameter space, where candidate draws are generated.logf
: function that takes a singleDenseVector
argument of parameter values at which to compute the logtransformed density (up to a normalizing constant).beta
: proportion of weight given to draws from the nonadaptive proposal with covariance factorizationSigmaL
, relative to draws from the adaptively tuned proposal with covariance factorizationSigmaLm
, during adaptive updating.scale
: factor (scale^2 / length(x)
) by which the adaptively updated covariance matrix is scaled—default value adopted from Gelman, Roberts, and Gilks [32].
Value
Returns anAMMVariate
type object with fields set to the suppliedx
and tuning parameter values.
AMMTune Type¶
Declaration¶
type AMMTune <: SamplerTune
Fields¶
logf::Nullable{Function}
: function supplied to the constructor to compute the logtransformed density, or null if not supplied.adapt::Bool
: whether the proposal distribution is being adaptively tuned.beta::Float64
: proportion of weight given to draws from the nonadaptive proposal with covariance factorizationSigmaL
, relative to draws from the adaptively tuned proposal with covariance factorizationSigmaLm
, during adaptive updating.m::Int
: number of adaptive update iterations that have been performed.Mv::Vector{Float64}
: running mean of drawsv
during adaptive updating. Used in the calculation ofSigmaLm
.Mvv::Matrix{Float64}
: running mean ofv * v'
during adaptive updating. Used in the calculation ofSigmaLm
.scale::Float64
: factor (scale^2 / length(v)
) by which the adaptively updated covariance matrix is scaled.SigmaL::LowerTriangular{Float64}
: Cholesky factorization of the nonadaptive covariance matrix.SigmaLm::Matrix{Float64}
: pivoted factorization of the adaptively tuned covariance matrix.
Adaptive Metropolis within Gibbs (AMWG)¶
Implementation of a MetropoliswithinGibbs sampler [63][79][95] for iteratively simulating autocorrelated draws from a distribution that can be specified up to a constant of proportionality.
ModelBased Constructor¶

AMWG
(params::ElementOrVector{Symbol}, sigma::ElementOrVector{T<:Real}; adapt::Symbol=:all, args...)¶ Construct a
Sampler
object for AMWG sampling. Parameters are assumed to be continuous, but may be constrained or unconstrained.Arguments
params
: stochastic node(s) to be updated with the sampler. Constrained parameters are mapped to unconstrained space according to transformations defined by the Stochasticunlist()
function.sigma
: scaling value or vector of the same length as the combined elements of nodesparams
, defining initial standard deviations for univariate normal proposal distributions. Standard deviations are relative to the unconstrained parameter space, where candidate draws are generated.adapt
: type of adaptation phase. Options are:all
: adapt proposals during all iterations.:burnin
: adapt proposals during burnin iterations.:none
: no adaptation (MetropoliswithinGibbs sampling with fixed proposals).
args...
: additional keyword arguments to be passed to theAMWGVariate
constructor.
Value
Returns aSampler{AMWGTune}
type object.Example
StandAlone Function¶

sample!
(v::AMWGVariate; adapt::Bool=true)¶ Draw one sample from a target distribution using the AMWG sampler. Parameters are assumed to be continuous and unconstrained.
Arguments
v
: current state of parameters to be simulated. When running the sampler in adaptive mode, thev
argument in a successive call to the function will contain thetune
field returned by the previous call.adapt
: whether to adaptively update the proposal distribution.
Value
Returnsv
updated with simulated values and associated tuning parameters.Example
The following example samples parameters in a simple linear regression model. Details of the model specification and posterior distribution can be found in the Supplement. Also, see the Line: BlockSpecific Sampling with AMWG and Slice example.
################################################################################ ## Linear Regression ## y ~ N(b0 + b1 * x, s2) ## b0, b1 ~ N(0, 1000) ## s2 ~ invgamma(0.001, 0.001) ################################################################################ using Mamba ## Data data = Dict( :x => [1, 2, 3, 4, 5], :y => [1, 3, 3, 3, 5] ) ## Logtransformed Posterior(b0, b1, log(s2)) + Constant logf = function(x::DenseVector) b0 = x[1] b1 = x[2] logs2 = x[3] r = data[:y]  b0  b1 * data[:x] (0.5 * length(data[:y])  0.001) * logs2  (0.5 * dot(r, r) + 0.001) / exp(logs2)  0.5 * b0^2 / 1000  0.5 * b1^2 / 1000 end ## MCMC Simulation with Adaptive MetopoliswithinGibbs Sampling n = 5000 burnin = 1000 sim = Chains(n, 3, names = ["b0", "b1", "s2"]) theta = AMWGVariate([0.0, 0.0, 0.0], 1.0, logf) for i in 1:n sample!(theta, adapt = (i <= burnin)) sim[i, :, 1] = [theta[1:2]; exp(theta[3])] end describe(sim)
AMWGVariate Type¶
Declaration¶
typealias AMWGVariate SamplerVariate{AMWGTune}
Fields¶
value::Vector{Float64}
: simulated values.tune::AMWGTune
: tuning parameters for the sampling algorithm.
Constructor¶

AMWGVariate
(x::AbstractVector{T<:Real}, sigma::ElementOrVector{U<:Real}, logf::Function; batchsize::Integer=50, target::Real=0.44)¶ Construct an
AMWGVariate
object that stores simulated values and tuning parameters for AMWG sampling.Arguments
x
: initial values.sigma
: scaling value or vector of the same length as the combined elements of nodesparams
, defining initial standard deviations for univariate normal proposal distributions. Standard deviations are relative to the unconstrained parameter space, where candidate draws are generated.logf
: function that takes a singleDenseVector
argument of parameter values at which to compute the logtransformed density (up to a normalizing constant).batchsize
: number of samples that must be accumulated before applying an adaptive update to the proposal distributions.target
: target acceptance rate for the algorithm.
Value
Returns anAMWGVariate
type object with fields set to the suppliedx
and tuning parameter values.
AMWGTune Type¶
Declaration¶
type AMWGTune <: SamplerTune
Fields¶
logf::Nullable{Function}
: function supplied to the constructor to compute the logtransformed density, or null if not supplied.adapt::Bool
: whether the proposal distribution is being adaptively tuned.accept::Vector{Int}
: number of accepted candidate draws generated for each element of the parameter vector during adaptive updating.batchsize::Int
: number of samples that must be accumulated before applying an adaptive update to the proposal distributions.m::Int
: number of adaptive update iterations that have been performed.sigma::Vector{Float64}
: updated values of the proposal standard deviations ifm > 0
, and usersupplied values otherwise.target::Float64
: target acceptance rate for the adaptive algorithm.
Binary Hamiltonian Monte Carlo (BHMC)¶
Implementation of the binarystate Hamiltonian Monte Carlo sampler of Pakman [68]. The sampler simulates autocorrelated draws from a distribution that can be specified up to a constant of proportionality.
ModelBased Constructor¶

BHMC
(params::ElementOrVector{Symbol}, traveltime::Real)¶ Construct a
Sampler
object for BHMC sampling. Parameters are assumed to have binary numerical values (0 or 1).Arguments
params
: stochastic node(s) to be updated with the sampler.traveltime
: length of time over which particle paths are simulated. It is recommended that supplied values be of the form , where optimal choices of are expected to grow with the parameter space dimensionality.
Value
Returns aSampler{BHMCTune}
type object.Example
StandAlone Function¶

sample!
(v::BHMCVariate)¶ Draw one sample from a target distribution using the BHMC sampler. Parameters are assumed to have binary numerical values (0 or 1).
Arguments
v
: current state of parameters to be simulated.
Value
Returnsv
updated with simulated values and associated tuning parameters.Example
################################################################################ ## Linear Regression ## y ~ MvNormal(X * (beta0 .* gamma), 1) ## gamma ~ DiscreteUniform(0, 1) ################################################################################ using Mamba ## Data n, p = 25, 10 X = randn(n, p) beta0 = randn(p) gamma0 = rand(0:1, p) y = X * (beta0 .* gamma0) + randn(n) ## Logtransformed Posterior(gamma) + Constant logf = function(gamma::DenseVector) logpdf(MvNormal(X * (beta0 .* gamma), 1.0), y) end ## MCMC Simulation with Binary Hamiltonian Monte Carlo Sampling t = 10000 sim = Chains(t, p, names = map(i > "gamma[$i]", 1:p)) gamma = BHMCVariate(zeros(p), (2 * p + 0.5) * pi, logf) for i in 1:t sample!(gamma) sim[i, :, 1] = gamma end describe(sim)
BHMCVariate Type¶
Declaration¶
typealias BHMCVariate SamplerVariate{BHMCTune}
Fields¶
value::Vector{Float64}
: simulated values.tune::BHMCTune
: tuning parameters for the sampling algorithm.
Constructor¶

BHMCVariate
(x::AbstractVector{T<:Real}, traveltime::Real, logf::Function)¶ Construct a
BHMCVariate
object that stores simulated values and tuning parameters for BHMC sampling.Arguments
x
: initial values.traveltime
: length of time over which particle paths are simulated. It is recommended that supplied values be of the form , where optimal choices of are expected to grow with the parameter space dimensionality.logf
: function that takes a singleDenseVector
argument of parameter values at which to compute the logtransformed density (up to a normalizing constant).
Value
Returns aBHMCVariate
type object with fields set to the suppliedx
and tuning parameter values.
BHMCTune Type¶
Declaration¶
type BHMCTune <: SamplerTune
Fields¶
logf::Nullable{Function}
: function supplied to the constructor to compute the logtransformed density, or null if not supplied.traveltime::Float64
: length of time over which particle paths are simulated.position::Vector{Float64}
: initial particle positions.velocity::Vector{Float64}
: initial particle velocites.wallhits::Int
: number of times particles are reflected off the 0 threshold.wallcrosses::Int
: number of times particles travel through the threshold.
Binary Individual Adaptation (BIA)¶
Implementation of the binarystate Individual Adaptation sampler of Griffin, et al. [43] which adjusts a general proposal to the data. The sampler simulates autocorrelated draws from a distribution that can be specified up to a constant of proportionality.
ModelBased Constructor¶

BIA
(params::ElementOrVector{Symbol}; args...)¶ Construct a
Sampler
object for BIA sampling. Parameters are assumed to have binary numerical values (0 or 1).Arguments
params
: stochastic node(s) to be updated with the sampler.args...
: additional keyword arguments to be passed to theBIAVariate
constructor.
Value
Returns aSampler{BIATune}
type object.Example
StandAlone Function¶

sample!
(v::BIAVariate)¶ Draw one sample from a target distribution using the BIA sampler. Parameters are assumed to have binary numerical values (0 or 1).
Arguments
v
: current state of parameters to be simulated.
Value
Returnsv
updated with simulated values and associated tuning parameters.Example
################################################################################ ## Linear Regression ## y ~ MvNormal(X * (beta0 .* gamma), 1) ## gamma ~ DiscreteUniform(0, 1) ################################################################################ using Mamba ## Data n, p = 25, 10 X = randn(n, p) beta0 = randn(p) gamma0 = rand(0:1, p) y = X * (beta0 .* gamma0) + randn(n) ## Logtransformed Posterior(gamma) + Constant logf = function(gamma::DenseVector) logpdf(MvNormal(X * (beta0 .* gamma), 1.0), y) end ## MCMC Simulation with Binary Individual Adaptation Sampler t = 10000 sim = Chains(t, p, names = map(i > "gamma[$i]", 1:p)) gamma = BIAVariate(zeros(p), logf) for i in 1:t sample!(gamma) sim[i, :, 1] = gamma end describe(sim)
BIAVariate Type¶
Declaration¶
typealias BIAVariate SamplerVariate{BIATune}
Fields¶
value::Vector{Float64}
: simulated values.tune::BIATune
: tuning parameters for the sampling algorithm.
Constructor¶

BIAVariate
(x::AbstractVector{T<:Real}, logf::Function; A::Vector{Float64} = ones(x) / length(x), D::Vector{Float64} = ones(x) / length(x), epsilon::Real = 0.01 / length(x), decay::Real = 0.55, target::Real = 0.45)¶ Construct a
BIAVariate
object that stores simulated values and tuning parameters for BIA sampling.Arguments
x
: initial values.logf
: function that takes a singleDenseVector
argument of parameter values at which to compute the logtransformed density (up to a normalizing constant).A
: vector of probabilities to switch the elements ofx
from 0 to 1 (i.e. added).D
: vector of probabilities to switch elements from 1 to 0 (i.e. deleted).epsilon
: range(epsilon, 1  epsilon)
for the elements ofA
andD
, where0 < epsilon < 0.5
.decay
: rate of decay of the adaptation, where0.5 < decay <= 1.0
.target
: target mutation rate, where0.0 < target < 1.0
.
Value
Returns aBIAVariate
type object with fields set to the suppliedx
and tuning parameter values.
BIATune Type¶
Declaration¶
type BIATune <: SamplerTune
Fields¶
logf::Nullable{Function}
: function supplied to the constructor to compute the logtransformed density, or null if not supplied.A::Vector{Float64}
: vector of probabilities to switch from 0 to 1.D::Vector{Float64}
: vector of probabilities to switch from 1 to 0.epsilon::Float64
: range(epsilon, 1  epsilon)
for the elements ofA
andD
.decay::Float64
: rate of decay of the adaptation.target::Float64
: target mutation rate.iter::Int
: iteration number for adaptive updating.
Binary MCMC Model Composition (BMC3)¶
Implementation of the binarystate MCMC Model Composition of Madigan and York [58] in which proposed updates are always state changes. Liu [55] shows this sampler is more efficient than Gibbs sampling for a binary vector. Schafer [83][84] proposes a method for block updates of binary vectors using this sampler. The sampler simulates autocorrelated draws from a distribution that can be specified up to a constant of proportionality.
ModelBased Constructor¶

BMC3
(params::ElementOrVector{Symbol}; k::BMC3Form=1)¶ Construct a
Sampler
object for BMC3 sampling. Parameters are assumed to have binary numerical values (0 or 1).Arguments
params
: stochastic node(s) to be updated with the sampler.k
: number of parameters or vector of parameter indices to select at random for simultaneous updating in each call of the sampler.
Value
Returns aSampler{BMC3Tune{typeof(k)}}
type object.Example
StandAlone Function¶

sample!
(v::SamplerVariate{BMC3Tune{F<:BMC3Form}})¶ Draw one sample from a target distribution using the BMC3 sampler. Parameters are assumed to have binary numerical values (0 or 1).
Arguments
v
: current state of parameters to be simulated.
Value
Returnsv
updated with simulated values and associated tuning parameters.Example
################################################################################ ## Linear Regression ## y ~ MvNormal(X * (beta0 .* gamma), 1) ## gamma ~ DiscreteUniform(0, 1) ################################################################################ using Mamba ## Data n, p = 25, 10 X = randn(n, p) beta0 = randn(p) gamma0 = rand(0:1, p) y = X * (beta0 .* gamma0) + randn(n) ## Logtransformed Posterior(gamma) + Constant logf = function(gamma::DenseVector) logpdf(MvNormal(X * (beta0 .* gamma), 1.0), y) end ## MCMC Simulation with Binary MCMC Model Composition t = 10000 sim1 = Chains(t, p, names = map(i > "gamma[$i]", 1:p)) sim2 = Chains(t, p, names = map(i > "gamma[$i]", 1:p)) gamma1 = BMC3Variate(zeros(p), logf) gamma2 = BMC3Variate(zeros(p), logf, k=Vector{Int}[[i] for i in 1:p]) for i in 1:t sample!(gamma1) sample!(gamma2) sim1[i, :, 1] = gamma1 sim2[i, :, 1] = gamma2 end describe(sim1) describe(sim2) p = plot(sim1, [:trace, :mixeddensity]) draw(p, filename = "bmc3plot")
BMC3 Variate Type¶
Declaration¶
SamplerVariate{BMC3Tune{F<:BMC3Form}}
Fields¶
value::Vector{Float64}
: simulated values.tune::BMC3Tune{F}
: tuning parameters for the sampling algorithm.
Constructor¶

BMC3Variate
(x::AbstractVector{T<:Real}, logf::Function; k::BMC3Form=1)¶ Construct a
SamplerVariate
object that stores simulated values and tuning parameters for BMC3 sampling.Arguments
x
: initial values.logf
: function that takes a singleDenseVector
argument of parameter values at which to compute the logtransformed density (up to a normalizing constant).k
: number of parameters or vector of parameter indices to select at random for simultaneous updating in each call of the sampler.
Value
Returns aSamplerVariate{BMC3Tune{typeof(k)}}
type object with fields set to the suppliedx
and tuning parameter values.
BMC3Tune Type¶
Declaration¶
typealias BMC3Form Union{Int, Vector{Vector{Int}}}
type BMC3Tune{F<:BMC3Form} <: SamplerTune
Fields¶
logf::Nullable{Function}
: function supplied to the constructor to compute the logtransformed density, or null if not supplied.k::F
: number of parameters or vector of parameter indices to select at random for simultaneous updating in each call of the sampler.
Binary Metropolised Gibbs (BMG)¶
Implementation of the binarystate Metropolised Gibbs sampler described by Schafer [83][84] in which components are drawn sequentially from full conditional marginal distributions and accepted together in a single MetropolisHastings step. The sampler simulates autocorrelated draws from a distribution that can be specified up to a constant of proportionality.
ModelBased Constructor¶

BMG
(params::ElementOrVector{Symbol}; k::BMGForm=1)¶ Construct a
Sampler
object for BMG sampling. Parameters are assumed to have binary numerical values (0 or 1).Arguments
params
: stochastic node(s) to be updated with the sampler.k
: number of parameters or vector of parameter indices to select at random for simultaneous updating in each call of the sampler.
Value
Returns aSampler{BMGTune{typeof(k)}}
type object.Example
StandAlone Function¶

sample!
(v::SamplerVariate{BMGTune{F<:BMGForm}})¶ Draw one sample from a target distribution using the BMG sampler. Parameters are assumed to have binary numerical values (0 or 1).
Arguments
v
: current state of parameters to be simulated.
Value
Returnsv
updated with simulated values and associated tuning parameters.Example
################################################################################ ## Linear Regression ## y ~ MvNormal(X * (beta0 .* gamma), 1) ## gamma ~ DiscreteUniform(0, 1) ################################################################################ using Mamba ## Data n, p = 25, 10 X = randn(n, p) beta0 = randn(p) gamma0 = rand(0:1, p) y = X * (beta0 .* gamma0) + randn(n) ## Logtransformed Posterior(gamma) + Constant logf = function(gamma::DenseVector) logpdf(MvNormal(X * (beta0 .* gamma), 1.0), y) end ## MCMC Simulation with Binary Metropolised Gibbs t = 10000 sim1 = Chains(t, p, names = map(i > "gamma[$i]", 1:p)) sim2 = Chains(t, p, names = map(i > "gamma[$i]", 1:p)) gamma1 = BMGVariate(zeros(p), logf) gamma2 = BMGVariate(zeros(p), logf, k=Vector{Int}[[i] for i in 1:p]) for i in 1:t sample!(gamma1) sample!(gamma2) sim1[i, :, 1] = gamma1 sim2[i, :, 1] = gamma2 end describe(sim1) describe(sim2)
BMG Variate Type¶
Declaration¶
SamplerVariate{BMGTune{F<:BMGForm}}
Fields¶
value::Vector{Float64}
: simulated values.tune::BMGTune{F}
: tuning parameters for the sampling algorithm.
Constructor¶

BMGVariate
(x::AbstractVector{T<:Real}, logf::Function; k::BMGForm=1)¶ Construct a
SamplerVariate
object that stores simulated values and tuning parameters for BMG sampling.Arguments
x
: initial values.logf
: function that takes a singleDenseVector
argument of parameter values at which to compute the logtransformed density (up to a normalizing constant).k
: number of parameters or vector of parameter indices to select at random for simultaneous updating in each call of the sampler.
Value
Returns aSamplerVariate{BMGTune{typeof(k)}}
type object with fields set to the suppliedx
and tuning parameter values.
BMGTune Type¶
Declaration¶
typealias BMGForm Union{Int, Vector{Vector{Int}}}
type BMGTune{F<:BMGForm} <: SamplerTune
Fields¶
logf::Nullable{Function}
: function supplied to the constructor to compute the logtransformed density, or null if not supplied.k::F
: number of parameters or vector of parameter indices to select at random for simultaneous updating in each call of the sampler.
Discrete Gibbs Sampler (DGS)¶
Implementation of a sampler for the simulation of discrete or discretized model parameters with finite support. Draws are simulated directly from a probability mass function that can be specified up to a constant of proportionality. Note that versions of this sampler evaluate the probability function over all points in the parameter space; and, as a result, may be very computationally intensive for large spaces.
ModelBased Constructor¶

DGS
(params::ElementOrVector{Symbol})¶ Construct a
Sampler
object for which DGS sampling is to be applied separately to each of the supplied parameters. Parameters are assumed to have discrete univariate distributions with finite supports.Arguments
params
: stochastic node(s) to be updated with the sampler.
Value
Returns aSampler{DSTune{Function}}
type object.Example
StandAlone Functions¶

sample!
(v::DGSVariate)¶ 
sample!
(v::DiscreteVariate) Draw one sample directly from a target probability mass function. Parameters are assumed to have discrete and finite support.
Arguments
v
: current state of parameters to be simulated.
Value
Returnsv
updated with simulated values and associated tuning parameters.
Discrete Variate Types¶
Declaration¶
typealias DGSVariate SamplerVariate{DSTune{Function}}
typealias DiscreteVariate SamplerVariate{DSTune{Vector{Float64}}}
Fields¶
value::Vector{Float64}
: simulated values.tune::DSTune{F<:DSForm}
: tuning parameters for the sampling algorithm.
Constructors¶

DGSVariate
(x::AbstractVector{T<:Real}, support::Matrix{U<:Real}, mass::Function)¶ 
DiscreteVariate
(x::AbstractVector{T<:Real}, support::Matrix{U<:Real}, mass::Vector{Float64})¶ Construct an object that stores simulated values and tuning parameters for discrete sampling.
Arguments
x
: initial values.support
: matrix whose columns contain the vector coordinates in the parameter space from which to simulate values.mass
: function that takes a singleDenseVector
argument of parameter values at which to compute the density (up to a normalizing constant), or a vector of sampling probabilities for the parameter space.
Value
Returns aDGSVariate
orDiscreteVariate
type object with fields set to the suppliedx
and tuning parameter values.
DSTune Type¶
Declaration¶
typealias DSForm Union{Function, Vector{Float64}}
type DSTune{F<:DSForm} <: SamplerTune
Fields¶
mass::Nullable{F}
: density mass function or vector supplied to the constructor, or null if not supplied.support::Matrix{Real}
: matrix whose columns contain the vector coordinates in the parameter space from which to simulate values.
Hamiltonian Monte Carlo (HMC)¶
Implementation of the Hybrid Monte Carlo (also known as Hamiltonian Monte Carlo) of Duane [23]. The sampler simulates autocorrelated draws from a distribution that can be specified up to a constant of proportionality. Code is derived from Neal’s implementation [65].
ModelBased Constructors¶

HMC
(params::ElementOrVector{Symbol}, epsilon::Real, L::Integer; dtype::Symbol=:forward)¶ 
HMC
(params::ElementOrVector{Symbol}, epsilon::Real, L::Integer, Sigma::Matrix{T<:Real}; dtype::Symbol=:forward) Construct a
Sampler
object for HMC sampling. Parameters are assumed to be continuous, but may be constrained or unconstrained.Arguments
params
: stochastic node(s) to be updated with the sampler. Constrained parameters are mapped to unconstrained space according to transformations defined by the Stochasticunlist()
function.epsilon
: step size.L
: number of steps to take in the Leapfrog algorithm.Sigma
: covariance matrix for the multivariate normal proposal distribution. The covariance matrix is relative to the unconstrained parameter space, where candidate draws are generated. If omitted, the identity matrix is assumed.dtype
: type of differentiation for gradient calculations. Options are:central
: central differencing.:forward
: forward differencing.
Value
Returns aSampler{HMCTune}
type object.Example
StandAlone Function¶

sample!
(v::HMCVariate)¶ Draw one sample from a target distribution using the HMC sampler. Parameters are assumed to be continuous and unconstrained.
Arguments
v
: current state of parameters to be simulated.
Value
Returnsv
updated with simulated values and associated tuning parameters.Example
The following example samples parameters in a simple linear regression model. Details of the model specification and posterior distribution can be found in the Supplement.
################################################################################ ## Linear Regression ## y ~ N(b0 + b1 * x, s2) ## b0, b1 ~ N(0, 1000) ## s2 ~ invgamma(0.001, 0.001) ################################################################################ using Mamba ## Data data = Dict( :x => [1, 2, 3, 4, 5], :y => [1, 3, 3, 3, 5] ) ## Logtransformed Posterior(b0, b1, log(s2)) + Constant and Gradient Vector logfgrad = function(x::DenseVector) b0 = x[1] b1 = x[2] logs2 = x[3] r = data[:y]  b0  b1 * data[:x] logf = (0.5 * length(data[:y])  0.001) * logs2  (0.5 * dot(r, r) + 0.001) / exp(logs2)  0.5 * b0^2 / 1000  0.5 * b1^2 / 1000 grad = [ sum(r) / exp(logs2)  b0 / 1000, sum(data[:x] .* r) / exp(logs2)  b1 / 1000, 0.5 * length(data[:y])  0.001 + (0.5 * dot(r, r) + 0.001) / exp(logs2) ] logf, grad end ## MCMC Simulation with Hamiltonian Monte Carlo ## Without (1) and with (2) a userspecified proposal covariance matrix n = 5000 sim1 = Chains(n, 3, names = ["b0", "b1", "s2"]) sim2 = Chains(n, 3, names = ["b0", "b1", "s2"]) epsilon = 0.1 L = 50 Sigma = eye(3) theta1 = HMCVariate([0.0, 0.0, 0.0], epsilon, L, logfgrad) theta2 = HMCVariate([0.0, 0.0, 0.0], epsilon, L, Sigma, logfgrad) for i in 1:n sample!(theta1) sample!(theta2) sim1[i, :, 1] = [theta1[1:2]; exp(theta1[3])] sim2[i, :, 1] = [theta2[1:2]; exp(theta2[3])] end describe(sim1) describe(sim2)
HMCVariate Type¶
Declaration¶
typealias HMCVariate SamplerVariate{HMCTune}
Fields¶
value::Vector{Float64}
: simulated values.tune::HMCTune
: tuning parameters for the sampling algorithm.
Constructors¶

HMCVariate
(x::AbstractVector{T<:Real}, epsilon::Real, L::Integer, logfgrad::Function)¶ 
HMCVariate
(x::AbstractVector{T<:Real}, epsilon::Real, L::Integer, Sigma::Matrix{U<:Real}, logfgrad::Function) Construct an
HMCVariate
object that stores simulated values and tuning parameters for HMC sampling.Arguments
x
: initial values.epsilon
: step size.L
: number of steps to take in the Leapfrog algorithm.Sigma
: covariance matrix for the multivariate normal proposal distribution. The covariance matrix is relative to the unconstrained parameter space, where candidate draws are generated. If omitted, the identity matrix is assumed.logfgrad
: function that takes a singleDenseVector
argument at which to compute the logtransformed density (up to a normalizing constant) and gradient vector, and returns the respective results as a tuple.
Value
Returns anHMCVariate
type object with fields set to the suppliedx
and tuning parameter values.
HMCTune Type¶
Declaration¶
type HMCTune <: SamplerTune
Fields¶
logfgrad::Nullable{Function}
: function supplied to the constructor to compute the logtransformed density and gradient vector, or null if not supplied.epsilon::Float64
: step size.L::Int
: number of steps to take in the Leapfrog algorithm.SigmaL::Union{UniformScaling{Int}, LowerTriangular{Float64}}
: Cholesky factorization of the covariance matrix for the multivariate normal proposal distribution.
MetropolisAdjusted Langevin Algorithm (MALA)¶
Implementation of the MetropolisAdjusted Langevin Algorithm of Roberts and Tweedie [81] and Roberts and Stramer [80]. The sampler simulates autocorrelated draws from a distribution that can be specified up to a constant of proportionality. MALA is related to Hamiltonian Monte Carlo as described thoroughly by Girolami and Calderhead [40].
ModelBased Constructors¶

MALA
(params::ElementOrVector{Symbol}, epsilon::Real; dtype::Symbol=:forward)¶ 
MALA
(params::ElementOrVector{Symbol}, epsilon::Real, Sigma::Matrix{T<:Real}; dtype::Symbol=:forward) Construct a
Sampler
object for MALA sampling. Parameters are assumed to be continuous, but may be constrained or unconstrained.Arguments
params
: stochastic node(s) to be updated with the sampler. Constrained parameters are mapped to unconstrained space according to transformations defined by the Stochasticunlist()
function.epsilon
: factor by which the drift and covariance matrix of the proposal distribution are scaled.Sigma
: covariance matrix for the multivariate normal proposal distribution. The covariance matrix is relative to the unconstrained parameter space, where candidate draws are generated. If omitted, the identity matrix is assumed.dtype
: type of differentiation for gradient calculations. Options are:central
: central differencing.:forward
: forward differencing.
Value
Returns aSampler{MALATune}
type object.Example
StandAlone Function¶

sample!
(v::MALAVariate)¶ Draw one sample from a target distribution using the MALA sampler. Parameters are assumed to be continuous and unconstrained.
Arguments
v
: current state of parameters to be simulated.
Value
Returnsv
updated with simulated values and associated tuning parameters.Example
The following example samples parameters in a simple linear regression model. Details of the model specification and posterior distribution can be found in the Supplement.
################################################################################ ## Linear Regression ## y ~ N(b0 + b1 * x, s2) ## b0, b1 ~ N(0, 1000) ## s2 ~ invgamma(0.001, 0.001) ################################################################################ using Mamba ## Data data = Dict( :x => [1, 2, 3, 4, 5], :y => [1, 3, 3, 3, 5] ) ## Logtransformed Posterior(b0, b1, log(s2)) + Constant and Gradient Vector logfgrad = function(x::DenseVector) b0 = x[1] b1 = x[2] logs2 = x[3] r = data[:y]  b0  b1 * data[:x] logf = (0.5 * length(data[:y])  0.001) * logs2  (0.5 * dot(r, r) + 0.001) / exp(logs2)  0.5 * b0^2 / 1000  0.5 * b1^2 / 1000 grad = [ sum(r) / exp(logs2)  b0 / 1000, sum(data[:x] .* r) / exp(logs2)  b1 / 1000, 0.5 * length(data[:y])  0.001 + (0.5 * dot(r, r) + 0.001) / exp(logs2) ] logf, grad end ## MCMC Simulation with MetropolisAdjusted Langevin Algorithm ## Without (1) and with (2) a userspecified proposal covariance matrix n = 5000 sim1 = Chains(n, 3, names = ["b0", "b1", "s2"]) sim2 = Chains(n, 3, names = ["b0", "b1", "s2"]) epsilon = 0.1 Sigma = eye(3) theta1 = MALAVariate([0.0, 0.0, 0.0], epsilon, logfgrad) theta2 = MALAVariate([0.0, 0.0, 0.0], epsilon, Sigma, logfgrad) for i in 1:n sample!(theta1) sample!(theta2) sim1[i, :, 1] = [theta1[1:2]; exp(theta1[3])] sim2[i, :, 1] = [theta2[1:2]; exp(theta2[3])] end describe(sim1) describe(sim2)
MALAVariate Type¶
Declaration¶
typealias MALAVariate SamplerVariate{MALATune}
Fields¶
value::Vector{Float64}
: simulated values.tune::MALATune
: tuning parameters for the sampling algorithm.
Constructors¶

MALAVariate
(x::AbstractVector{T<:Real}, epsilon::Real, logfgrad::Function)¶ 
MALAVariate
(x::AbstractVector{T<:Real}, epsilon::Real, Sigma::Matrix{U<:Real}, logfgrad::Function) Construct a
MALAVariate
object that stores simulated values and tuning parameters for MALA sampling.Arguments
x
: initial values.epsilon
: factor by which the drift and covariance matrix of the proposal distribution are scaled.Sigma
: covariance matrix for the multivariate normal proposal distribution. The covariance matrix is relative to the unconstrained parameter space, where candidate draws are generated. If omitted, the identity matrix is assumed.logfgrad
: function that takes a singleDenseVector
argument of parameter values at which to compute the logtransformed density (up to a normalizing constant) and gradient vector, and returns the respective results as a tuple.
Value
Returns aMALAVariate
type object with fields set to the suppliedx
and tuning parameter values.
MALATune Type¶
Declaration¶
type MALATune <: SamplerTune
Fields¶
logfgrad::Nullable{Function}
: function supplied to the constructor to compute the logtransformed density and gradient vector, or null if not supplied.epsilon::Float64
: factor by which the drift and covariance matrix of the proposal distribution are scaled.SigmaL::Union{UniformScaling{Int}, LowerTriangular{Float64}}
: Cholesky factorization of the covariance matrix for the multivariate normal proposal distribution.
Missing Values Sampler (MISS)¶
A sampler to simulate missing output values from their likelihood distributions.
ModelBased Constructor¶

MISS
(params::ElementOrVector{Symbol})¶ Construct a
Sampler
object to sampling missing output values. The constructor should only be used to sample stochastic nodes upon which no other stochastic node depends. Socalled ‘output nodes’ can be identified with thekeys()
function. Moreover, when theMISS
constructor is included in a vector ofSampler
objects to define a sampling scheme, it should be positioned at the beginning of the vector. This ensures that missing output values are updated before any other samplers are executed.Arguments
params
: stochastic node(s) that contain missing values (NaN
) to be updated with the sampler.
Value
Returns aSampler{Dict{Symbol, MISSTune}}
type object.Example
NoUTurn Sampler (NUTS)¶
Implementation of the NoUTurn Sampler extension (algorithm 6) [48] to Hamiltonian Monte Carlo [65] for simulating autocorrelated draws from a distribution that can be specified up to a constant of proportionality.
ModelBased Constructor¶

NUTS
(params::ElementOrVector{Symbol}; dtype::Symbol=:forward, args...)¶ Construct a
Sampler
object for NUTS sampling, with the algorithm’s step size parameter adaptively tuned during burnin iterations. Parameters are assumed to be continuous, but may be constrained or unconstrained.Arguments
params
: stochastic node(s) to be updated with the sampler. Constrained parameters are mapped to unconstrained space according to transformations defined by the Stochasticunlist()
function.dtype
: type of differentiation for gradient calculations. Options are:central
: central differencing.:forward
: forward differencing.
args...
: additional keyword arguments to be passed to theNUTSVariate
constructor.
Value
Returns aSampler{NUTSTune}
type object.Example
StandAlone Function¶

sample!
(v::NUTSVariate; adapt::Bool=false)¶ Draw one sample from a target distribution using the NUTS sampler. Parameters are assumed to be continuous and unconstrained.
Arguments
v
: current state of parameters to be simulated. When running the sampler in adaptive mode, thev
argument in a successive call to the function will contain thetune
field returned by the previous call.adapt
: whether to adaptively update theepsilon
step size parameter.
Value
Returnsv
updated with simulated values and associated tuning parameters.Example
The following example samples parameters in a simple linear regression model. Details of the model specification and posterior distribution can be found in the Supplement.
################################################################################ ## Linear Regression ## y ~ N(b0 + b1 * x, s2) ## b0, b1 ~ N(0, 1000) ## s2 ~ invgamma(0.001, 0.001) ################################################################################ using Mamba ## Data data = Dict( :x => [1, 2, 3, 4, 5], :y => [1, 3, 3, 3, 5] ) ## Logtransformed Posterior(b0, b1, log(s2)) + Constant and Gradient Vector logfgrad = function(x::DenseVector) b0 = x[1] b1 = x[2] logs2 = x[3] r = data[:y]  b0  b1 * data[:x] logf = (0.5 * length(data[:y])  0.001) * logs2  (0.5 * dot(r, r) + 0.001) / exp(logs2)  0.5 * b0^2 / 1000  0.5 * b1^2 / 1000 grad = [ sum(r) / exp(logs2)  b0 / 1000, sum(data[:x] .* r) / exp(logs2)  b1 / 1000, 0.5 * length(data[:y])  0.001 + (0.5 * dot(r, r) + 0.001) / exp(logs2) ] logf, grad end ## MCMC Simulation with NoUTurn Sampling n = 5000 burnin = 1000 sim = Chains(n, 3, start = (burnin + 1), names = ["b0", "b1", "s2"]) theta = NUTSVariate([0.0, 0.0, 0.0], logfgrad) for i in 1:n sample!(theta, adapt = (i <= burnin)) if i > burnin sim[i, :, 1] = [theta[1:2]; exp(theta[3])] end end describe(sim)
NUTSVariate Type¶
Declaration¶
typealias NUTSVariate SamplerVariate{NUTSTune}
Fields¶
value::Vector{Float64}
: simulated values.tune::NUTSTune
: tuning parameters for the sampling algorithm.
Constructors¶

NUTSVariate
(x::AbstractVector{T<:Real}, epsilon::Real, logfgrad::Function; target::Real=0.6)¶ 
NUTSVariate
(x::AbstractVector{T<:Real}, logfgrad::Function; target::Real=0.6) Construct a
NUTSVariate
object that stores simulated values and tuning parameters for NUTS sampling.Arguments
x
: initial values.epsilon
: step size parameter.logfgrad
: function that takes a singleDenseVector
argument of parameter values at which to compute the logtransformed density (up to a normalizing constant) and gradient vector, and returns the respective results as a tuple. Ifepsilon
is not specified, the function is used by the constructor to generate an initial step size value.target
: target acceptance rate for the algorithm.
Value
Returns aNUTSVariate
type object with fields set to the suppliedx
and tuning parameter values.
NUTSTune Type¶
Declaration¶
type NUTSTune <: SamplerTune
Fields¶
logfgrad::Nullable{Function}
: function supplied to the constructor to compute the logtransformed density and gradient vector, or null if not supplied.adapt::Bool
: whether the proposal distribution is being adaptively tuned.alpha::Float64
: cumulative acceptance probabilities from leapfrog steps.epsilon::Float64
: updated value of the step size parameter ifm > 0
, and the usersupplied value otherwise.epsbar::Float64
: dual averaging parameter, defined as .gamma::Float64
: dual averaging parameter, fixed at .Hbar::Float64
: dual averaging parameter, defied as .kappa::Float64
: dual averaging parameter, fixed at .m::Int
: number of adaptive update iterations that have been performed.mu::Float64
: dual averaging parameter, defined as .nalpha::Int
: the total number of leapfrog steps performed.t0::Float64
: dual averaging parameter, fixed at .target::Float64
: target acceptance rate for the adaptive algorithm.
Random Walk Metropolis (RWM)¶
Random walk MetropolisHastings algorithm [46][63] in which parameters are sampled from symmetric distributions centered around the current values. The sampler simulates autocorrelated draws from a distribution that can be specified up to a constant of proportionality.
ModelBased Constructor¶

RWM
(params::ElementOrVector{Symbol}, scale::ElementOrVector{T<:Real}; args...)¶ Construct a
Sampler
object for RWM sampling. Parameters are assumed to be continuous, but may be constrained or unconstrained.Arguments
params
: stochastic node(s) to be updated with the sampler. Constrained parameters are mapped to unconstrained space according to transformations defined by the Stochasticunlist()
function.scale
: scaling value or vector of the same length as the combined elements of nodesparams
for theproposal
distribution. Values are relative to the unconstrained parameter space, where candidate draws are generated.args...
: additional keyword arguments to be passed to theRWMVariate
constructor.
Value
Returns aSampler{RWMTune}
type object.Example
StandAlone Function¶

sample!
(v::RWMVariate)¶ Draw one sample from a target distribution using the RWM sampler. Parameters are assumed to be continuous and unconstrained.
Arguments
v
: current state of parameters to be simulated.
Value
Returnsv
updated with simulated values and associated tuning parameters.Example
The following example samples parameters in a simple linear regression model. Details of the model specification and posterior distribution can be found in the Supplement.
################################################################################ ## Linear Regression ## y ~ N(b0 + b1 * x, s2) ## b0, b1 ~ N(0, 1000) ## s2 ~ invgamma(0.001, 0.001) ################################################################################ using Mamba ## Data data = Dict( :x => [1, 2, 3, 4, 5], :y => [1, 3, 3, 3, 5] ) ## Logtransformed Posterior(b0, b1, log(s2)) + Constant logf = function(x::DenseVector) b0 = x[1] b1 = x[2] logs2 = x[3] r = data[:y]  b0  b1 * data[:x] (0.5 * length(data[:y])  0.001) * logs2  (0.5 * dot(r, r) + 0.001) / exp(logs2)  0.5 * b0^2 / 1000  0.5 * b1^2 / 1000 end ## MCMC Simulation with Random Walk Metropolis n = 5000 burnin = 1000 sim = Chains(n, 3, names = ["b0", "b1", "s2"]) theta = RWMVariate([0.0, 0.0, 0.0], [0.5, 0.25, 1.0], logf, proposal = SymUniform) for i in 1:n sample!(theta) sim[i, :, 1] = [theta[1:2]; exp(theta[3])] end describe(sim)
RWMVariate Type¶
Declaration¶
typealias RWMVariate SamplerVariate{RWMTune}
Fields¶
value::Vector{Float64}
: simulated values.tune::RWMTune
: tuning parameters for the sampling algorithm.
Constructor¶

RWMVariate
(x::AbstractVector{T<:Real}, scale::ElementOrVector{U<:Real}, logf::Function; proposal::SymDistributionType=Normal)¶ Construct a
RWMVariate
object that stores simulated values and tuning parameters for RWM sampling.Arguments
x
: initial values.scale
: scalar or vector of the same length asx
for theproposal
distribution.logf
: function that takes a singleDenseVector
argument of parameter values at which to compute the logtransformed density (up to a normalizing constant).proposal
: symmetric distribution of typeBiweight
,Cosine
,Epanechnikov
,Normal
,SymTriangularDist
,SymUniform
, orTriweight
to be centered around current parameter values and used to generate proposal draws. Specifiedscale
determines the standard deviations of Normal proposals and widths of the others.
Value
Returns aRWMVariate
type object with fields set to the suppliedx
and tuning parameter values.
RWMTune Type¶
Declaration¶
type RWMTune <: SamplerTune
Fields¶
logf::Nullable{Function}
: function supplied to the constructor to compute the logtransformed density, or null if not supplied.scale::Union{Float64, Vector{Float64}}
: scaling for the proposal distribution.proposal::SymDistributionType
: proposal distribution.
Shrinkage Slice (Slice)¶
Implementation of the shrinkage slice sampler of Neal [64] for simulating autocorrelated draws from a distribution that can be specified up to a constant of proportionality.
ModelBased Constructor¶

Slice
(params::ElementOrVector{Symbol}, width::ElementOrVector{T<:Real}, ::Type{F<:SliceForm}=Multivariate; transform::Bool=false)¶ Construct a
Sampler
object for Slice sampling. Parameters are assumed to be continuous, but may be constrained or unconstrained.Arguments
params
: stochastic node(s) to be updated with the sampler.width
: scaling value or vector of the same length as the combined elements of nodesparams
, defining initial widths of a hyperrectangle from which to simulate values.F
: sampler type. Options areUnivariate
: sequential univariate sampling of parameters .Multivariate
: joint multivariate sampling.
transform
: whether to sample parameters on the linktransformed scale (unconstrained parameter space). Iftrue
, then constrained parameters are mapped to unconstrained space according to transformations defined by the Stochasticunlist()
function, andwidth
is interpreted as being relative to the unconstrained parameter space. Otherwise, sampling is relative to the untransformed space.
Value
Returns aSampler{SliceTune{Univariate}}
orSampler{SliceTune{Multivariate}}
type object if sampling univariately or multivariately, respectively.Example
StandAlone Functions¶

sample!
(v::SliceUnivariate)¶ 
sample!
(v::SliceMultivariate) Draw one sample from a target distribution using the Slice univariate or multivariate sampler. Parameters are assumed to be continuous, but may be constrained or unconstrained.
Arguments
v
: current state of parameters to be simulated.
Value
Returnsv
updated with simulated values and associated tuning parameters.Example
The following example samples parameters in a simple linear regression model. Details of the model specification and posterior distribution can be found in the Supplement. Also, see the Line: BlockSpecific Sampling with AMWG and Slice example.
################################################################################ ## Linear Regression ## y ~ N(b0 + b1 * x, s2) ## b0, b1 ~ N(0, 1000) ## s2 ~ invgamma(0.001, 0.001) ################################################################################ using Mamba ## Data data = Dict( :x => [1, 2, 3, 4, 5], :y => [1, 3, 3, 3, 5] ) ## Logtransformed Posterior(b0, b1, log(s2)) + Constant logf = function(x::DenseVector) b0 = x[1] b1 = x[2] logs2 = x[3] r = data[:y]  b0  b1 * data[:x] (0.5 * length(data[:y])  0.001) * logs2  (0.5 * dot(r, r) + 0.001) / exp(logs2)  0.5 * b0^2 / 1000  0.5 * b1^2 / 1000 end ## MCMC Simulation with Slice Sampling ## With multivariate (1) and univariate (2) updating n = 5000 sim1 = Chains(n, 3, names = ["b0", "b1", "s2"]) sim2 = Chains(n, 3, names = ["b0", "b1", "s2"]) width = [1.0, 1.0, 2.0] theta1 = SliceUnivariate([0.0, 0.0, 0.0], width, logf) theta2 = SliceMultivariate([0.0, 0.0, 0.0], width, logf) for i in 1:n sample!(theta1) sample!(theta2) sim1[i, :, 1] = [theta1[1:2]; exp(theta1[3])] sim2[i, :, 1] = [theta2[1:2]; exp(theta2[3])] end describe(sim1) describe(sim2)
Slice Variate Types¶
Declaration¶
typealias SliceUnivariate SamplerVariate{SliceTune{Univariate}}
typealias SliceMultivariate SamplerVariate{SliceTune{Multivariate}}
Fields¶
value::Vector{Float64}
: simulated values.tune::SliceTune{F<:SliceForm}
: tuning parameters for the sampling algorithm.
Constructors¶

SliceUnivariate
(x::AbstractVector{T<:Real}, width::ElementOrVector{U<:Real}, logf::Function)¶ 
SliceMultivariate
(x::AbstractVector{T<:Real}, width::ElementOrVector{U<:Real}, logf::Function)¶ Construct an object that stores simulated values and tuning parameters for Slice sampling.
Arguments
x
: initial values.width
: scaling value or vector of the same length as the combined elements of nodesparams
, defining initial widths of a hyperrectangle from which to simulate values.logf
: function that takes a singleDenseVector
argument of parameter values at which to compute the logtransformed density (up to a normalizing constant).
Value
Returns an object of the same type as the constructor name for univariate or multivariate sampling, respectively, with fields set to the suppliedx
and tuning parameter values.
SliceTune Type¶
Declaration¶
typealias SliceForm Union{Univariate, Multivariate}
type SliceTune{F<:SliceForm} <: SamplerTune
Fields¶
logf::Nullable{Function}
: function supplied to the constructor to compute the logtransformed density, or null if not supplied.width::Union{Float64, Vector{Float64}}
: initial widths of hyperrectangles from which to simulate values.
Slice Simplex (SliceSimplex)¶
Implementation of the slice simplex sampler as described by Cowles et al. [19] for simulating autocorrelated draws of parameters on the simplex and from a distribution that can be specified up to a constant of proportionality.
ModelBased Constructor¶

SliceSimplex
(params::ElementOrVector{Symbol}; args...)¶ Construct a
Sampler
object for which SliceSimplex sampling is to be applied separately to each of the supplied parameters. Parameters are assumed to be continuous and constrained to a simplex.Arguments
params
: stochastic node(s) to be updated with the sampler.args...
: additional keyword arguments to be passed to theSliceSimplexVariate
constructor.
Value
Returns aSampler{SliceSimplexTune}
type object.Example
StandAlone Function¶

sample!
(v::SliceSimplexVariate)¶ Draw one sample from a target distribution using the SliceSimplex sampler. Parameters are assumed to be continuous and constrained to a simplex.
Arguments
v
: current state of parameters to be simulated.
Value
Returnsv
updated with simulated values and associated tuning parameters.Example
################################################################################ ## Multinomial Model ## y ~ Multinomial(n, rho) ## rho ~ Dirichlet(1, ..., 1) ################################################################################ using Mamba ## Data n, k = 100, 5 rho0 = rand(Dirichlet(ones(k))) y = rand(Multinomial(n, rho0)) ## Logtransformed Posterior(rho) + Constant logf = function(rho::DenseVector) logpdf(Multinomial(n, rho), y) end ## MCMC Simulation with Slice Simplex Sampling t = 10000 sim = Chains(t, k, names = map(i > "rho[$i]", 1:k)) rho = SliceSimplexVariate(fill(1 / k, k), logf) for i in 1:t sample!(rho) sim[i, :, 1] = rho end describe(sim) p = plot(sim) draw(p, filename = "slicesimplexplot")
SliceSimplexVariate Type¶
Declaration¶
typealias SliceSimplexVariate SamplerVariate{SliceSimplexTune}
Fields¶
value::Vector{Float64}
: simulated values.tune::SliceSimplexTune
: tuning parameters for the sampling algorithm.
Constructor¶

SliceSimplexVariate
(x::AbstractVector{T<:Real}, logf::Function; scale::Real=1.0)¶ Construct a
SliceSimplexVariate
object that stores simulated values and tuning parameters for slice simplex sampling.Arguments
x
: initial values.scale
: value0 < scale <= 1
by which to scale the standard simplex to define an initial space from which to simulate values.logf
: function that takes a singleDenseVector
argument of parameter values at which to compute the logtransformed density (up to a normalizing constant).
Value
Returns aSliceSimplexVariate
type object with fields set to the suppliedx
and tuning parameter values.
SliceSimplexTune Type¶
Declaration¶
type SliceSimplexTune <: SamplerTune
Fields¶
logf::Nullable{Function}
: function supplied to the constructor to compute the logtransformed density, or null if not supplied.scale::Float64
: value0 < scale <= 1
by which to scale the standard simplex to define an initial space from which to simulate values.
The following table summarizes the (ddimensional) sample spaces over which each method simulates draws, whether draws are generated univariately or multivariately, and whether transformations are applied to map parameters to the sample spaces.
ModelBased Constructors  StandAlone Functions  

Method  Sample Space  Univariate  Multivariate  Transformations  Univariate  Multivariate 
ABC  No  Yes  Yes  No  No  
AMM  No  Yes  Yes  No  Yes  
AMWG  Yes  No  Yes  Yes  No  
BHMC  No  Yes  No  No  Yes  
BIA  No  Yes  No  No  Yes  
BMC3  Yes  Yes  No  Yes  Yes  
BMG  Yes  Yes  No  Yes  Yes  
DGS  Finite  Yes  No  No  No  Yes 
HMC  No  Yes  Yes  No  Yes  
MALA  No  Yes  Yes  No  Yes  
MISS  Parameterdefined  Yes  Yes  No  No  No 
NUTS  No  Yes  Yes  No  Yes  
RWM  No  Yes  Yes  No  Yes  
Slice  Yes  Yes  Optional  Yes  Yes  
SliceSimplex  dsimplex  No  Yes  No  No  Yes 
Examples¶
OpenBUGS¶
The following examples are taken from OpenBUGS [44], and were used in the development and testing of Mamba. They are provide to illustrate model specification and fitting with the package, and how its syntax compares to other Bayesian modelling software.
Volume I¶
Rats: A Normal Hierarchical Model¶
An example from OpenBUGS [44] and section 6 of Gelfand et al. [29] concerning 30 rats whose weights were measured at each of five consecutive weeks.
Weights are modeled as
where is repeated weight measurement on rat , and is the day on which the measurement was taken.
using Mamba
## Data
rats = Dict{Symbol, Any}(
:y =>
[151, 199, 246, 283, 320,
145, 199, 249, 293, 354,
147, 214, 263, 312, 328,
155, 200, 237, 272, 297,
135, 188, 230, 280, 323,
159, 210, 252, 298, 331,
141, 189, 231, 275, 305,
159, 201, 248, 297, 338,
177, 236, 285, 350, 376,
134, 182, 220, 260, 296,
160, 208, 261, 313, 352,
143, 188, 220, 273, 314,
154, 200, 244, 289, 325,
171, 221, 270, 326, 358,
163, 216, 242, 281, 312,
160, 207, 248, 288, 324,
142, 187, 234, 280, 316,
156, 203, 243, 283, 317,
157, 212, 259, 307, 336,
152, 203, 246, 286, 321,
154, 205, 253, 298, 334,
139, 190, 225, 267, 302,
146, 191, 229, 272, 302,
157, 211, 250, 285, 323,
132, 185, 237, 286, 331,
160, 207, 257, 303, 345,
169, 216, 261, 295, 333,
157, 205, 248, 289, 316,
137, 180, 219, 258, 291,
153, 200, 244, 286, 324],
:x => [8.0, 15.0, 22.0, 29.0, 36.0]
)
rats[:xbar] = mean(rats[:x])
rats[:N] = size(rats[:y], 1)
rats[:T] = size(rats[:y], 2)
rats[:rat] = Int[div(i  1, 5) + 1 for i in 1:150]
rats[:week] = Int[(i  1) % 5 + 1 for i in 1:150]
rats[:X] = rats[:x][rats[:week]]
rats[:Xm] = rats[:X]  rats[:xbar]
## Model Specification
model = Model(
y = Stochastic(1,
(alpha, beta, rat, Xm, s2_c) >
begin
mu = alpha[rat] + beta[rat] .* Xm
MvNormal(mu, sqrt(s2_c))
end,
false
),
alpha = Stochastic(1,
(mu_alpha, s2_alpha) > Normal(mu_alpha, sqrt(s2_alpha)),
false
),
alpha0 = Logical(
(mu_alpha, xbar, mu_beta) > mu_alpha  xbar * mu_beta
),
mu_alpha = Stochastic(
() > Normal(0.0, 1000),
false
),
s2_alpha = Stochastic(
() > InverseGamma(0.001, 0.001),
false
),
beta = Stochastic(1,
(mu_beta, s2_beta) > Normal(mu_beta, sqrt(s2_beta)),
false
),
mu_beta = Stochastic(
() > Normal(0.0, 1000)
),
s2_beta = Stochastic(
() > InverseGamma(0.001, 0.001),
false
),
s2_c = Stochastic(
() > InverseGamma(0.001, 0.001)
)
)
## Initial Values
inits = [
Dict(:y => rats[:y], :alpha => fill(250, 30), :beta => fill(6, 30),
:mu_alpha => 150, :mu_beta => 10, :s2_c => 1, :s2_alpha => 1,
:s2_beta => 1),
Dict(:y => rats[:y], :alpha => fill(20, 30), :beta => fill(0.6, 30),
:mu_alpha => 15, :mu_beta => 1, :s2_c => 10, :s2_alpha => 10,
:s2_beta => 10)
]
## Sampling Scheme
scheme = [Slice(:s2_c, 10.0),
AMWG(:alpha, 100.0),
Slice([:mu_alpha, :s2_alpha], [100.0, 10.0], Univariate),
AMWG(:beta, 1.0),
Slice([:mu_beta, :s2_beta], 1.0, Univariate)]
setsamplers!(model, scheme)
## MCMC Simulations
sim = mcmc(model, rats, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)
Iterations = 2502:10000
Thinning interval = 2
Chains = 1,2
Samples per chain = 3750
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
s2_c 37.2543133 6.026634572 0.0695895819 0.2337982327 664.4576
mu_beta 6.1830663 0.108042927 0.0012475723 0.0017921615 3634.4474
alpha0 106.6259925 3.459210115 0.0399435178 0.0526804390 3750.0000
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
s2_c 27.778388 33.0906026 36.4630047 40.5538472 51.5713716
mu_beta 5.969850 6.1110307 6.1836454 6.2538831 6.3964953
alpha0 99.815707 104.3369878 106.6105679 108.9124224 113.5045347
Pumps: GammaPoisson Hierarchical Model¶
An example from OpenBUGS [44] and George et al. [36] concerning the number of failures of 10 power plant pumps.
Pump failure are modelled as
where is the number of times that pump failed, and is the operation time of the pump (in 1000s of hours).
using Mamba
## Data
pumps = Dict{Symbol, Any}(
:y => [5, 1, 5, 14, 3, 19, 1, 1, 4, 22],
:t => [94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5]
)
pumps[:N] = length(pumps[:y])
## Model Specification
model = Model(
y = Stochastic(1,
(theta, t, N) >
UnivariateDistribution[
begin
lambda = theta[i] * t[i]
Poisson(lambda)
end
for i in 1:N
],
false
),
theta = Stochastic(1,
(alpha, beta) > Gamma(alpha, 1 / beta),
true
),
alpha = Stochastic(
() > Exponential(1.0)
),
beta = Stochastic(
() > Gamma(0.1, 1.0)
)
)
## Initial Values
inits = [
Dict(:y => pumps[:y], :alpha => 1.0, :beta => 1.0,
:theta => rand(Gamma(1.0, 1.0), pumps[:N])),
Dict(:y => pumps[:y], :alpha => 10.0, :beta => 10.0,
:theta => rand(Gamma(10.0, 10.0), pumps[:N]))
]
## Sampling Scheme
scheme = [Slice([:alpha, :beta], 1.0, Univariate),
Slice(:theta, 1.0, Univariate)]
setsamplers!(model, scheme)
## MCMC Simulations
sim = mcmc(model, pumps, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)
## Posterior Predictive Distribution
ppd = predict(sim, :y)
describe(ppd)
## MCMC Simulations
Iterations = 2502:10000
Thinning interval = 2
Chains = 1,2
Samples per chain = 3750
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
beta 0.93036099 0.517274134 0.00597296721 0.01824153419 804.11618
alpha 0.69679849 0.264420102 0.00305326034 0.00722593007 1339.06428
theta[1] 0.05991674 0.025050913 0.00028926303 0.00032725274 3750.00000
theta[2] 0.10125873 0.078887354 0.00091091270 0.00129985769 3683.18182
theta[3] 0.08910382 0.037307731 0.00043079257 0.00045660986 3750.00000
theta[4] 0.11529009 0.030274265 0.00034957710 0.00032583333 3750.00000
theta[5] 0.59971611 0.316811730 0.00365822675 0.00585119652 2931.65686
theta[6] 0.60967188 0.134761371 0.00155609028 0.00174949908 3750.00000
theta[7] 0.86767451 0.669943846 0.00773584520 0.02858200254 549.40360
theta[8] 0.85445727 0.668132036 0.00771492422 0.02485547216 722.57105
theta[9] 1.55721556 0.753564449 0.00870141275 0.03109274798 587.38464
theta[10] 1.98475207 0.405438843 0.00468160451 0.00912748779 1973.09587
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
beta 0.189558591 0.55034346 0.84013720 1.22683690 2.134324768
alpha 0.286434561 0.50214938 0.66295652 0.85427430 1.319308645
theta[1] 0.021291808 0.04164372 0.05696335 0.07435302 0.118041776
theta[2] 0.008789962 0.04279772 0.08116810 0.13771833 0.305675244
theta[3] 0.032503106 0.06202730 0.08293498 0.11024076 0.176883533
theta[4] 0.063587574 0.09389172 0.11222268 0.13443866 0.182555751
theta[5] 0.153474125 0.36947945 0.54009121 0.77234190 1.364481338
theta[6] 0.371096082 0.51395467 0.60089332 0.69531957 0.898866079
theta[7] 0.077416146 0.37391993 0.71230582 1.17659703 2.616100803
theta[8] 0.072479432 0.35973235 0.69319639 1.16347299 2.655857786
theta[9] 0.463821785 1.01169918 1.43745818 1.96871605 3.351798915
theta[10] 1.269842527 1.70167020 1.95748757 2.24000075 2.861197982
## Posterior Predictive Distribution
Iterations = 2502:10000
Thinning interval = 2
Chains = 1,2
Samples per chain = 3750
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
y[1] 5.5784000 3.3457654 0.038633571 0.042295418 3750.0000
y[2] 1.6009333 1.7532320 0.020244579 0.022399464 3750.0000
y[3] 5.5854667 3.2976980 0.038078537 0.040503717 3750.0000
y[4] 14.5666667 5.4547523 0.062986054 0.061718435 3750.0000
y[5] 3.1118667 2.4165260 0.027903638 0.035965796 3750.0000
y[6] 19.0780000 6.0320229 0.069651801 0.071700830 3750.0000
y[7] 0.9008000 1.1756177 0.013574864 0.031805581 1366.2355
y[8] 0.8922667 1.1550634 0.013337523 0.025610086 2034.1808
y[9] 3.2585333 2.4039269 0.027758156 0.069046712 1212.1504
y[10] 20.7933333 6.2068800 0.071670876 0.113975236 2965.6896
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
y[1] 1 3 5 8 13
y[2] 0 0 1 2 6
y[3] 1 3 5 7 13
y[4] 5 11 14 18 27
y[5] 0 1 3 4 9
y[6] 9 15 19 23 32
y[7] 0 0 1 1 4
y[8] 0 0 1 1 4
y[9] 0 1 3 5 9
y[10] 10 16 20 25 34
Dogs: Loglinear Model for Binary Data¶
An example from OpenBUGS [44], Lindley and Smith [54], and Kalbfleisch [52] concerning the SolomonWynne experiment on dogs. In the experiment, 30 dogs were subjected to 25 trials. On each trial, a barrier was raised, and an electric shock was administered 10 seconds later if the dog did not jump the barrier.
Failures to jump the barriers in time are modelled as
where if dog fails to jump the barrier before the shock on trial , and 0 otherwise; is the number of successful jumps prior to trial ; and is the probability of failure.
using Mamba
## Data
dogs = Dict{Symbol, Any}(
:Y =>
[0 0 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1
0 0 0 0 0 1 1 0 1 1 0 0 1 1 0 1 0 1 1 1 1 1 1 1 1
0 1 1 0 0 1 1 1 1 0 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1
0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 0 0 0 0 1 1 1 1 0 0 1 0 1 1 1 1 1 1 1 1 1 1 1
0 0 0 0 0 1 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 0 0 0 0 0 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 0 0 0 1 0 1 0 1 1 0 1 0 0 0 1 1 1 1 1 0 1 1 0
0 0 0 0 1 0 0 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 0 0 0 1 1 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 0 1 1 0 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 0 0 1 0 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 0 1 0 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 1 0 1 0 0 0 1 0 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1
0 0 0 0 1 0 1 0 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1
0 1 0 0 0 0 1 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 0 0 1 1 0 1 0 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1
0 0 0 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 1 0 1 0 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 1 1 1 1
0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 0 0 0 0 0 0 1 1 1 0 1 0 0 0 1 1 0 1 1 1 1 1 1
0 0 0 0 0 0 1 1 0 1 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1
0 0 1 0 1 1 1 0 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 0 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 0 1 0 1 0 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1
0 0 0 0 1 1 0 0 1 1 1 0 1 0 1 0 1 0 1 1 1 1 1 1 1
0 0 0 0 1 1 1 1 1 1 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1]
)
dogs[:Dogs] = size(dogs[:Y], 1)
dogs[:Trials] = size(dogs[:Y], 2)
dogs[:xa] = mapslices(cumsum, dogs[:Y], 2)
dogs[:xs] = mapslices(x > collect(1:25)  x, dogs[:xa], 2)
dogs[:y] = 1  dogs[:Y][:, 2:25]
## Model Specification
model = Model(
y = Stochastic(2,
(Dogs, Trials, alpha, xa, beta, xs) >
UnivariateDistribution[
begin
p = exp(alpha * xa[i, j] + beta * xs[i, j])
Bernoulli(p)
end
for i in 1:Dogs, j in 1:Trials1
],
false
),
alpha = Stochastic(
() > Truncated(Flat(), Inf, 1e5)
),
A = Logical(
alpha > exp(alpha)
),
beta = Stochastic(
() > Truncated(Flat(), Inf, 1e5)
),
B = Logical(
beta > exp(beta)
)
)
## Initial Values
inits = [
Dict(:y => dogs[:y], :alpha => 1, :beta => 1),
Dict(:y => dogs[:y], :alpha => 2, :beta => 2)
]
## Sampling Scheme
scheme = [Slice([:alpha, :beta], 1.0)]
setsamplers!(model, scheme)
## MCMC Simulations
sim = mcmc(model, dogs, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)
Iterations = 2502:10000
Thinning interval = 2
Chains = 1,2
Samples per chain = 3750
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
beta 0.0788601 0.011812880 0.00013640339 0.00018435162 3750.0000
B 0.9242336 0.010903201 0.00012589932 0.00017012187 3750.0000
alpha 0.2441654 0.024064384 0.00027787158 0.00045249868 2828.2310
A 0.7835846 0.018821132 0.00021732772 0.00035466850 2816.0882
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
beta 0.10326776 0.086633844 0.07848501 0.07046663 0.05703080
B 0.90188545 0.917012804 0.92451592 0.93195884 0.94456498
alpha 0.29269315 0.260167645 0.24381478 0.22784322 0.19872975
A 0.74625110 0.770922334 0.78363276 0.79624909 0.81977141
Seeds: Random Effect Logistic Regression¶
An example from OpenBUGS [44], Crowder [20], and Breslow and Clayton [10] concerning the proportion of seeds that germinated on each of 21 plates arranged according to a 2 by 2 factorial layout by seed and type of root extract.
Germinations are modelled as
where are the number of seeds, out of , that germinate on plate ; and and are the seed type and root extract.
using Mamba
## Data
seeds = Dict{Symbol, Any}(
:r => [10, 23, 23, 26, 17, 5, 53, 55, 32, 46, 10, 8, 10, 8, 23, 0, 3, 22, 15,
32, 3],
:n => [39, 62, 81, 51, 39, 6, 74, 72, 51, 79, 13, 16, 30, 28, 45, 4, 12, 41,
30, 51, 7],
:x1 => [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
:x2 => [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1]
)
seeds[:N] = length(seeds[:r])
## Model Specification
model = Model(
r = Stochastic(1,
(alpha0, alpha1, x1, alpha2, x2, alpha12, b, n, N) >
UnivariateDistribution[
begin
p = invlogit(alpha0 + alpha1 * x1[i] + alpha2 * x2[i] +
alpha12 * x1[i] * x2[i] + b[i])
Binomial(n[i], p)
end
for i in 1:N
],
false
),
b = Stochastic(1,
s2 > Normal(0, sqrt(s2)),
false
),
alpha0 = Stochastic(
() > Normal(0, 1000)
),
alpha1 = Stochastic(
() > Normal(0, 1000)
),
alpha2 = Stochastic(
() > Normal(0, 1000)
),
alpha12 = Stochastic(
() > Normal(0, 1000)
),
s2 = Stochastic(
() > InverseGamma(0.001, 0.001)
)
)
## Initial Values
inits = [
Dict(:r => seeds[:r], :alpha0 => 0, :alpha1 => 0, :alpha2 => 0,
:alpha12 => 0, :s2 => 0.01, :b => zeros(seeds[:N])),
Dict(:r => seeds[:r], :alpha0 => 0, :alpha1 => 0, :alpha2 => 0,
:alpha12 => 0, :s2 => 1, :b => zeros(seeds[:N]))
]
## Sampling Scheme
scheme = [AMM([:alpha0, :alpha1, :alpha2, :alpha12], 0.01 * eye(4)),
AMWG(:b, 0.01),
AMWG(:s2, 0.1)]
setsamplers!(model, scheme)
## MCMC Simulations
sim = mcmc(model, seeds, inits, 12500, burnin=2500, thin=2, chains=2)
describe(sim)
Iterations = 2502:12500
Thinning interval = 2
Chains = 1,2
Samples per chain = 5000
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
alpha2 1.310728093 0.26053104 0.0026053104 0.0153996801 286.21707
alpha1 0.088700176 0.26872879 0.0026872879 0.0128300598 438.70341
alpha0 0.556154341 0.17595432 0.0017595432 0.0101730837 299.15388
alpha12 0.746440855 0.43006756 0.0043006756 0.0251658152 292.04607
s2 0.085705306 0.09738014 0.0009738014 0.0080848189 145.07755
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
alpha2 0.8040593795 1.148881827 1.309947687 1.48076318 1.82815608
alpha1 0.4250164771 0.093637900 0.094390643 0.26007581 0.62353933
alpha0 0.9149197759 0.666632319 0.551292851 0.44262420 0.22244775
alpha12 1.5457041398 1.027576522 0.757250262 0.49149187 0.17029702
s2 0.0011739822 0.021117624 0.059376315 0.11140082 0.35234645
Surgical: Institutional Ranking¶
An example from OpenBUGS [44] concerning mortality rates in 12 hospitals performing cardiac surgery in infants.
Number of deaths are modelled as
where is the number of deaths, out of operations, at hospital .
using Mamba
## Data
surgical = Dict{Symbol, Any}(
:r => [0, 18, 8, 46, 8, 13, 9, 31, 14, 8, 29, 24],
:n => [47, 148, 119, 810, 211, 196, 148, 215, 207, 97, 256, 360]
)
surgical[:N] = length(surgical[:r])
## Model Specification
model = Model(
r = Stochastic(1,
(n, p, N) >
UnivariateDistribution[Binomial(n[i], p[i]) for i in 1:N],
false
),
p = Logical(1,
b > invlogit(b)
),
b = Stochastic(1,
(mu, s2) > Normal(mu, sqrt(s2)),
false
),
mu = Stochastic(
() > Normal(0, 1000)
),
pop_mean = Logical(
mu > invlogit(mu)
),
s2 = Stochastic(
() > InverseGamma(0.001, 0.001)
)
)
## Initial Values
inits = [
Dict(:r => surgical[:r], :b => fill(0.1, surgical[:N]), :s2 => 1, :mu => 0),
Dict(:r => surgical[:r], :b => fill(0.5, surgical[:N]), :s2 => 10, :mu => 1)
]
## Sampling Scheme
scheme = [NUTS(:b),
Slice([:mu, :s2], 1.0)]
setsamplers!(model, scheme)
## MCMC Simulations
sim = mcmc(model, surgical, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)
Iterations = 2502:10000
Thinning interval = 2
Chains = 1,2
Samples per chain = 3750
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
mu 2.550263247 0.151813688 0.001752993476 0.00352027397 1859.8115
pop_mean 0.073062651 0.010097974 0.000116601357 0.00022880854 1947.7088
s2 0.183080212 0.161182244 0.001861172237 0.00629499754 655.6065
p[1] 0.053571675 0.019444542 0.000224526231 0.00059140521 1080.9986
p[2] 0.103203725 0.022015796 0.000254216516 0.00049928289 1944.3544
p[3] 0.071050270 0.017662050 0.000203943787 0.00022426787 3750.0000
p[4] 0.059863573 0.008208359 0.000094781965 0.00033190971 611.6074
p[5] 0.052400628 0.013632842 0.000157418497 0.00052344316 678.3186
p[6] 0.069799258 0.014697820 0.000169715805 0.00026854081 2995.6099
p[7] 0.066927057 0.015703466 0.000181328008 0.00023888938 3750.0000
p[8] 0.122296440 0.023319424 0.000269269516 0.00086456417 727.5137
p[9] 0.070386216 0.014318572 0.000165336624 0.00019674962 3750.0000
p[10] 0.077978945 0.019775182 0.000228344129 0.00031544646 3750.0000
p[11] 0.101809685 0.017568938 0.000202868617 0.00037809859 2159.1405
p[12] 0.068534503 0.011579540 0.000133709009 0.00015162331 3750.0000
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
mu 2.878686485 2.637983754 2.540103608 2.455828143 2.269200544
pop_mean 0.053217279 0.066733498 0.073094153 0.079013392 0.093706084
s2 0.011890649 0.082427721 0.141031919 0.233947847 0.601137558
p[1] 0.017683165 0.039573948 0.052942713 0.067219732 0.092219913
p[2] 0.067482073 0.086541186 0.100283541 0.116521584 0.153268132
p[3] 0.040088134 0.059099562 0.070320337 0.080842394 0.110595078
p[4] 0.044965780 0.054145641 0.059405801 0.065149564 0.076754218
p[5] 0.028224773 0.042599948 0.051613569 0.061247876 0.079162530
p[6] 0.043031618 0.059518315 0.069130011 0.079935035 0.100122661
p[7] 0.038586653 0.055992241 0.066465865 0.076006492 0.101351899
p[8] 0.082604334 0.105852706 0.121313052 0.137724088 0.171202438
p[9] 0.044453383 0.060622475 0.070120723 0.079078759 0.101046045
p[10] 0.044439834 0.064694148 0.075946776 0.089236114 0.122899027
p[11] 0.071966143 0.088360531 0.100133691 0.112632331 0.139968412
p[12] 0.047466661 0.060611827 0.068267225 0.075638053 0.093309250
Magnesium: MetaAnalysis Prior Sensitivity¶
An example from OpenBUGS [44].
Number of events reported for treatment and control subjects in 8 studies is modelled as
where is the number of control group events, out of , in study ; is the number of treatment group events; and indexes differ prior specifications.
using Mamba
## Data
magnesium = Dict{Symbol, Any}(
:rt => [1, 9, 2, 1, 10, 1, 1, 90],
:nt => [40, 135, 200, 48, 150, 59, 25, 1159],
:rc => [2, 23, 7, 1, 8, 9, 3, 118],
:nc => [36, 135, 200, 46, 148, 56, 23, 1157]
)
magnesium[:rtx] = hcat([magnesium[:rt] for i in 1:6]...)'
magnesium[:rcx] = hcat([magnesium[:rc] for i in 1:6]...)'
magnesium[:s2] = 1 ./ (magnesium[:rt] + 0.5) +
1 ./ (magnesium[:nt]  magnesium[:rt] + 0.5) +
1 ./ (magnesium[:rc] + 0.5) +
1 ./ (magnesium[:nc]  magnesium[:rc] + 0.5)
magnesium[:s2_0] = 1 / mean(1 ./ magnesium[:s2])
## Model Specification
model = Model(
rcx = Stochastic(2,
(nc, pc) >
UnivariateDistribution[Binomial(nc[j], pc[i, j]) for i in 1:6, j in 1:8],
false
),
pc = Stochastic(2,
() > Uniform(0, 1),
false
),
rtx = Stochastic(2,
(nt, pc, theta) >
UnivariateDistribution[
begin
phi = logit(pc[i, j])
pt = invlogit(theta[i, j] + phi)
Binomial(nt[j], pt)
end
for i in 1:6, j in 1:8
],
false
),
theta = Stochastic(2,
(mu, tau) >
UnivariateDistribution[Normal(mu[i], tau[i]) for i in 1:6, j in 1:8],
false
),
mu = Stochastic(1,
() > Uniform(10, 10),
false
),
OR = Logical(1,
mu > exp(mu)
),
tau = Logical(1,
(priors, s2_0) >
Float64[
sqrt(priors[1]),
sqrt(priors[2]),
priors[3],
sqrt(s2_0 * (1 / priors[4]  1)),
sqrt(s2_0) * (1 / priors[5]  1),
sqrt(priors[6]) ]
),
priors = Stochastic(1,
s2_0 >
UnivariateDistribution[
InverseGamma(0.001, 0.001),
Uniform(0, 50),
Uniform(0, 50),
Uniform(0, 1),
Uniform(0, 1),
Truncated(Normal(0, sqrt(s2_0 / erf(0.75))), 0, Inf)
],
false
)
)
## Initial Values
inits = [
Dict(:rcx => magnesium[:rcx], :rtx => magnesium[:rtx],
:theta => zeros(6, 8), :mu => fill(0.5, 6),
:pc => fill(0.5, 6, 8), :priors => [1, 1, 1, 0.5, 0.5, 1]),
Dict(:rcx => magnesium[:rcx], :rtx => magnesium[:rtx],
:theta => zeros(6, 8), :mu => fill(0.5, 6),
:pc => fill(0.5, 6, 8), :priors => [1, 1, 1, 0.5, 0.5, 1])
]
## Sampling Scheme
scheme = [AMWG(:theta, 0.1),
AMWG(:mu, 0.1),
Slice(:pc, 0.25, Univariate),
Slice(:priors, [1.0, 5.0, 5.0, 0.25, 0.25, 5.0], Univariate)]
setsamplers!(model, scheme)
## MCMC Simulations
sim = mcmc(model, magnesium, inits, 12500, burnin=2500, thin=2, chains=2)
describe(sim)
Iterations = 2502:12500
Thinning interval = 2
Chains = 1,2
Samples per chain = 5000
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
tau[1] 0.55098858 0.35814901 0.0035814901 0.0221132365 262.31486
tau[2] 1.11557619 0.58886505 0.0058886505 0.0237788755 613.26606
tau[3] 0.83211110 0.49113676 0.0049113676 0.0222839957 485.75664
tau[4] 0.47864203 0.26258828 0.0026258828 0.0135868530 373.51920
tau[5] 0.48624861 0.35359386 0.0035359386 0.0215005369 270.46485
tau[6] 0.56841884 0.18877962 0.0018877962 0.0058505056 1041.17429
OR[1] 0.47784058 0.15389133 0.0015389133 0.0066922017 528.79852
OR[2] 0.42895913 0.32240192 0.0032240192 0.0081170895 1577.59150
OR[3] 0.43118350 0.18264467 0.0018264467 0.0064385836 804.69879
OR[4] 0.47587697 0.13947735 0.0013947735 0.0064893426 461.96170
OR[5] 0.48545299 0.14603013 0.0014603013 0.0083912319 302.85415
OR[6] 0.44554385 0.14121352 0.0014121352 0.0053818401 688.47941
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
tau[1] 0.050143630 0.28821905 0.49325476 0.73991793 1.42033542
tau[2] 0.326282292 0.70071249 0.98873505 1.39092470 2.65606546
tau[3] 0.136936046 0.49195419 0.74461372 1.06769203 2.03061696
tau[4] 0.091858771 0.28921566 0.44085251 0.61784984 1.10639558
tau[5] 0.028866916 0.23628679 0.42220429 0.65955432 1.37318750
tau[6] 0.214834142 0.43417486 0.56753402 0.70014948 0.94171060
OR[1] 0.206501871 0.37431326 0.47128600 0.57044051 0.80062151
OR[2] 0.107428346 0.27074745 0.38362516 0.52062237 0.99299623
OR[3] 0.145435475 0.30454141 0.41470065 0.53630024 0.83015778
OR[4] 0.231777387 0.38069803 0.47049713 0.56112444 0.76292805
OR[5] 0.207697044 0.38509504 0.48308284 0.59166588 0.75526778
OR[6] 0.208377218 0.34750042 0.43313882 0.52797553 0.76192141
Salm: ExtraPoisson Variation in a DoseResponse Study¶
An example from OpenBUGS [44] and Breslow [9] concerning mutagenicity assay data on salmonella in three plates exposed to six doses of quinoline.
Number of revertant colonies of salmonella are modelled as
where is the number of colonies in plate and dose .
using Mamba
## Data
salm = Dict{Symbol, Any}(
:y => reshape(
[15, 21, 29, 16, 18, 21, 16, 26, 33, 27, 41, 60, 33, 38, 41, 20, 27, 42],
3, 6),
:x => [0, 10, 33, 100, 333, 1000],
:plate => 3,
:dose => 6
)
## Model Specification
model = Model(
y = Stochastic(2,
(alpha, beta, gamma, x, lambda) >
UnivariateDistribution[
begin
mu = exp(alpha + beta * log(x[j] + 10) + gamma * x[j] + lambda[i, j])
Poisson(mu)
end
for i in 1:3, j in 1:6
],
false
),
alpha = Stochastic(
() > Normal(0, 1000)
),
beta = Stochastic(
() > Normal(0, 1000)
),
gamma = Stochastic(
() > Normal(0, 1000)
),
lambda = Stochastic(2,
s2 > Normal(0, sqrt(s2)),
false
),
s2 = Stochastic(
() > InverseGamma(0.001, 0.001)
)
)
## Initial Values
inits = [
Dict(:y => salm[:y], :alpha => 0, :beta => 0, :gamma => 0, :s2 => 10,
:lambda => zeros(3, 6)),
Dict(:y => salm[:y], :alpha => 1, :beta => 1, :gamma => 0.01, :s2 => 1,
:lambda => zeros(3, 6))
]
## Sampling Scheme
scheme = [Slice([:alpha, :beta, :gamma], [1.0, 1.0, 0.1]),
AMWG([:lambda, :s2], 0.1)]
setsamplers!(model, scheme)
## MCMC Simulations
sim = mcmc(model, salm, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)
Iterations = 2502:10000
Thinning interval = 2
Chains = 1,2
Samples per chain = 3750
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
s2 0.0690769709 0.04304237136 0.000497010494 0.001985478741 469.961085
gamma 0.0011250515 0.00034536546 0.000003987937 0.000025419899 184.590851
beta 0.3543443166 0.07160779229 0.000826855563 0.007122805926 101.069088
alpha 2.0100584321 0.26156942610 0.003020343571 0.027100522461 93.157673
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
s2 0.0135820292 0.039788136 0.0598307554 0.08740879470 0.17747300708
gamma 0.0017930379 0.001347435 0.0011325733 0.00090960998 0.00039271486
beta 0.2073237562 0.312747679 0.3582979574 0.40006921583 0.48789037325
alpha 1.5060295212 1.847111545 2.0062727893 2.16556036713 2.50547846044
Equiv: Bioequivalence in a CrossOver Trial¶
An example from OpenBUGS [44] and Gelfand et al. [29] concerning a twotreatment, crossover trial with 10 subjects.
Treatment responses are modelled as
where is the response for patient in period ; and is the treatment received.
using Mamba
## Data
equiv = Dict{Symbol, Any}(
:group => [1, 1, 2, 2, 2, 1, 1, 1, 2, 2],
:y =>
[1.40 1.65
1.64 1.57
1.44 1.58
1.36 1.68
1.65 1.69
1.08 1.31
1.09 1.43
1.25 1.44
1.25 1.39
1.30 1.52]
)
equiv[:N] = size(equiv[:y], 1)
equiv[:P] = size(equiv[:y], 2)
equiv[:T] = [equiv[:group] 3  equiv[:group]]
## Model Specification
model = Model(
y = Stochastic(2,
(delta, mu, phi, pi, s2_1, T) >
begin
sigma = sqrt(s2_1)
UnivariateDistribution[
begin
m = mu + (1)^(T[i, j]  1) * phi / 2 + (1)^(j  1) * pi / 2 +
delta[i, j]
Normal(m, sigma)
end
for i in 1:10, j in 1:2
]
end,
false
),
delta = Stochastic(2,
s2_2 > Normal(0, sqrt(s2_2)),
false
),
mu = Stochastic(
() > Normal(0, 1000)
),
phi = Stochastic(
() > Normal(0, 1000)
),
theta = Logical(
phi > exp(phi)
),
pi = Stochastic(
() > Normal(0, 1000)
),
s2_1 = Stochastic(
() > InverseGamma(0.001, 0.001)
),
s2_2 = Stochastic(
() > InverseGamma(0.001, 0.001)
),
equiv = Logical(
theta > Int(0.8 < theta < 1.2)
)
)
## Initial Values
inits = [
Dict(:y => equiv[:y], :delta => zeros(10, 2), :mu => 0, :phi => 0,
:pi => 0, :s2_1 => 1, :s2_2 => 1),
Dict(:y => equiv[:y], :delta => zeros(10, 2), :mu => 10, :phi => 10,
:pi => 10, :s2_1 => 10, :s2_2 => 10)
]
## Sampling Scheme
scheme = [NUTS(:delta),
Slice([:mu, :phi, :pi], 1.0),
Slice([:s2_1, :s2_2], 1.0, Univariate)]
setsamplers!(model, scheme)
## MCMC Simulations
sim = mcmc(model, equiv, inits, 12500, burnin=2500, thin=2, chains=2)
describe(sim)
Iterations = 2502:12500
Thinning interval = 2
Chains = 1,2
Samples per chain = 5000
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
s2_2 0.0173121833 0.014549568 0.00014549568 0.0007329722 394.02626
s2_1 0.0184397014 0.013837972 0.00013837972 0.0005689492 591.55873
pi 0.1874240524 0.086420302 0.00086420302 0.0032257037 717.76558
phi 0.0035569545 0.087590520 0.00087590520 0.0035141650 621.25503
theta 1.0002921934 0.088250458 0.00088250458 0.0036227671 593.40761
equiv 0.9751000000 0.155828169 0.00155828169 0.0036666529 1806.14385
mu 1.4387396416 0.042269208 0.00042269208 0.0013735876 946.96847
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
s2_2 0.0016375061 0.0056159514 0.013968228 0.024613730 0.053154674
s2_1 0.0017289114 0.0074958338 0.015849718 0.025963832 0.051967161
pi 0.3579631753 0.2432161807 0.187946319 0.130454165 0.014910965
phi 0.1723722017 0.0623573600 0.005681830 0.053144647 0.172913654
theta 0.8416658455 0.9395470702 0.994334281 1.054582177 1.188763456
equiv 1.0000000000 1.0000000000 1.000000000 1.000000000 1.000000000
mu 1.3552569594 1.4110400018 1.438593809 1.466525521 1.519643109
Dyes: Variance Components Model¶
An example from OpenBUGS [44], Davies [21], and Box and Tiao [8] concerning batchtobatch variation in yields from six batches and five samples of dyestuff.
using Mamba
## Data
dyes = Dict{Symbol, Any}(
:y =>
[1545, 1440, 1440, 1520, 1580,
1540, 1555, 1490, 1560, 1495,
1595, 1550, 1605, 1510, 1560,
1445, 1440, 1595, 1465, 1545,
1595, 1630, 1515, 1635, 1625,
1520, 1455, 1450, 1480, 1445],
:batches => 6,
:samples => 5
)
dyes[:batch] = vcat([fill(i, dyes[:samples]) for i in 1:dyes[:batches]]...)
dyes[:sample] = vcat(fill(collect(1:dyes[:samples]), dyes[:batches])...)
## Model Specification
model = Model(
y = Stochastic(1,
(mu, batch, s2_within) > MvNormal(mu[batch], sqrt(s2_within)),
false
),
mu = Stochastic(1,
(theta, batches, s2_between) > Normal(theta, sqrt(s2_between))
),
theta = Stochastic(
() > Normal(0, 1000)
),
s2_within = Stochastic(
() > InverseGamma(0.001, 0.001)
),
s2_between = Stochastic(
() > InverseGamma(0.001, 0.001)
)
)
## Initial Values
inits = [
Dict(:y => dyes[:y], :theta => 1500, :s2_within => 1, :s2_between => 1,
:mu => fill(1500, dyes[:batches])),
Dict(:y => dyes[:y], :theta => 3000, :s2_within => 10, :s2_between => 10,
:mu => fill(3000, dyes[:batches]))
]
## Sampling Schemes
scheme = [NUTS([:mu, :theta]),
Slice([:s2_within, :s2_between], 1000.0)]
scheme2 = [MALA(:theta, 50.0),
MALA(:mu, 50.0, eye(dyes[:batches])),
Slice([:s2_within, :s2_between], 1000.0)]
scheme3 = [HMC(:theta, 10.0, 5),
HMC(:mu, 10.0, 5, eye(dyes[:batches])),
Slice([:s2_within, :s2_between], 1000.0)]
scheme4 = [RWM(:theta, 50.0, proposal=Cosine),
RWM(:mu, 50.0),
Slice([:s2_within, :s2_between], 1000.0)]
## MCMC Simulations
setsamplers!(model, scheme)
sim = mcmc(model, dyes, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)
setsamplers!(model, scheme2)
sim2 = mcmc(model, dyes, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim2)
setsamplers!(model, scheme3)
sim3 = mcmc(model, dyes, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim3)
setsamplers!(model, scheme4)
sim4 = mcmc(model, dyes, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim4)
Iterations = 2502:10000
Thinning interval = 2
Chains = 1,2
Samples per chain = 3750
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
s2_between 3192.2614 4456.129102 51.45494673 487.44961486 83.57109
theta 1526.7186 24.549671 0.28347518 0.37724897 3750.00000
s2_within 2887.5853 1075.174260 12.41504296 76.89117959 195.52607
mu[1] 1511.4798 20.819711 0.24040531 0.52158448 1593.30921
mu[2] 1527.9087 20.344151 0.23491402 0.30199960 3750.00000
mu[3] 1552.6742 21.293738 0.24587891 0.70276515 918.08605
mu[4] 1506.6440 21.349176 0.24651905 0.61821290 1192.57616
mu[5] 1578.6636 25.512471 0.29459264 1.29216105 389.82685
mu[6] 1487.1934 24.693967 0.28514137 1.23710390 398.44592
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
s2_between 111.92351 815.9012 1651.4605 3269.5663 1.90261752x104
theta 1475.08243 1513.5687 1527.1763 1540.0550 1.57444106x103
s2_within 1566.61796 2160.3251 2654.9358 3324.8621 5.65161862x103
mu[1] 1469.68693 1498.0985 1511.3084 1525.7344 1.55281296x103
mu[2] 1486.15990 1514.8487 1527.6843 1541.2155 1.56770562x103
mu[3] 1512.72912 1537.7046 1552.1976 1566.7188 1.59498301x103
mu[4] 1463.91701 1492.3206 1506.9856 1521.3449 1.54854165x103
mu[5] 1528.52480 1562.1831 1579.2515 1596.0167 1.62731017x103
mu[6] 1440.27721 1470.7983 1486.1844 1502.9911 1.54464614x103
Stacks: Robust Regression¶
An example from OpenBUGS [44], Brownlee [14], and Birkes and Dodge [5] concerning 21 daily responses of stack loss, the amount of ammonia escaping, as a function of air flow, temperature, and acid concentration.
using Mamba
## Data
stacks = Dict{Symbol, Any}(
:y => [42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9,
15, 15],
:x =>
[80 27 89
80 27 88
75 25 90
62 24 87
62 22 87
62 23 87
62 24 93
62 24 93
58 23 87
58 18 80
58 18 89
58 17 88
58 18 82
58 19 93
50 18 89
50 18 86
50 19 72
50 19 79
50 20 80
56 20 82
70 20 91]
)
stacks[:N] = size(stacks[:x], 1)
stacks[:p] = size(stacks[:x], 2)
stacks[:meanx] = map(j > mean(stacks[:x][:, j]), 1:stacks[:p])
stacks[:sdx] = map(j > std(stacks[:x][:, j]), 1:stacks[:p])
stacks[:z] = Float64[
(stacks[:x][i, j]  stacks[:meanx][j]) / stacks[:sdx][j]
for i in 1:stacks[:N], j in 1:stacks[:p]
]
## Model Specification
model = Model(
y = Stochastic(1,
(mu, s2, N) >
UnivariateDistribution[Laplace(mu[i], s2) for i in 1:N],
false
),
beta0 = Stochastic(
() > Normal(0, 1000),
false
),
beta = Stochastic(1,
() > Normal(0, 1000),
false
),
mu = Logical(1,
(beta0, z, beta) > beta0 + z * beta,
false
),
s2 = Stochastic(
() > InverseGamma(0.001, 0.001),
false
),
sigma = Logical(
s2 > sqrt(2.0) * s2
),
b0 = Logical(
(beta0, b, meanx) > beta0  dot(b, meanx)
),
b = Logical(1,
(beta, sdx) > beta ./ sdx
),
outlier = Logical(1,
(y, mu, sigma, N) >
Float64[abs((y[i]  mu[i]) / sigma) > 2.5 for i in 1:N],
[1, 3, 4, 21]
)
)
## Initial Values
inits = [
Dict(:y => stacks[:y], :beta0 => 10, :beta => [0, 0, 0], :s2 => 10),
Dict(:y => stacks[:y], :beta0 => 1, :beta => [1, 1, 1], :s2 => 1)
]
## Sampling Scheme
scheme = [NUTS([:beta0, :beta]),
Slice(:s2, 1.0)]
setsamplers!(model, scheme)
## MCMC Simulations
sim = mcmc(model, stacks, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)
Iterations = 2502:10000
Thinning interval = 2
Chains = 1,2
Samples per chain = 3750
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
b[1] 0.836863707 0.13085145 0.0015109423 0.0027601754 2247.4171
b[2] 0.744454449 0.33480007 0.0038659382 0.0065756939 2592.3158
b[3] 0.116648437 0.12214077 0.0014103601 0.0015143922 3750.0000
b0 38.776564595 8.81860433 0.1018284717 0.0979006137 3750.0000
sigma 3.487643717 0.87610847 0.0101164292 0.0279025494 985.8889
outlier[1] 0.042666667 0.20211796 0.0023338572 0.0029490162 3750.0000
outlier[3] 0.054800000 0.22760463 0.0026281519 0.0034398827 3750.0000
outlier[4] 0.298000000 0.45740999 0.0052817156 0.0089200654 2629.5123
outlier[21] 0.606400000 0.48858046 0.0056416412 0.0113877443 1840.7583
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
b[1] 0.57218621 0.75741345 0.834874964 0.918345319 1.101502854
b[2] 0.16177144 0.52291878 0.714951465 0.933171533 1.476258382
b[3] 0.36401372 0.19028697 0.113463801 0.036994963 0.118538277
b0 56.70056875 44.11785905 38.698338454 33.409149788 21.453323631
sigma 2.17947513 2.86899865 3.348631697 3.953033535 5.592773118
outlier[1] 0.00000000 0.00000000 0.000000000 0.000000000 1.000000000
outlier[3] 0.00000000 0.00000000 0.000000000 0.000000000 1.000000000
outlier[4] 0.00000000 0.00000000 0.000000000 1.000000000 1.000000000
outlier[21] 0.00000000 0.00000000 1.000000000 1.000000000 1.000000000
Epilepsy: Repeated Measures on Poisson Counts¶
An example from OpenBUGS [44], Thall and Vail [92] Breslow and Clayton [10] concerning the effects of treatment, baseline seizure counts, and age on followup seizure counts at four visits in 59 patients.
Counts are modelled as
where are the counts on patient at visit , is a treatment indicator, is baseline seizure counts, is age in years, and is an indicator for the fourth visit.
using Mamba
## Data
epil = Dict{Symbol, Any}(
:y =>
[ 5 3 3 3
3 5 3 3
2 4 0 5
4 4 1 4
7 18 9 21
5 2 8 7
6 4 0 2
40 20 21 12
5 6 6 5
14 13 6 0
26 12 6 22
12 6 8 4
4 4 6 2
7 9 12 14
16 24 10 9
11 0 0 5
0 0 3 3
37 29 28 29
3 5 2 5
3 0 6 7
3 4 3 4
3 4 3 4
2 3 3 5
8 12 2 8
18 24 76 25
2 1 2 1
3 1 4 2
13 15 13 12
11 14 9 8
8 7 9 4
0 4 3 0
3 6 1 3
2 6 7 4
4 3 1 3
22 17 19 16
5 4 7 4
2 4 0 4
3 7 7 7
4 18 2 5
2 1 1 0
0 2 4 0
5 4 0 3
11 14 25 15
10 5 3 8
19 7 6 7
1 1 2 3
6 10 8 8
2 1 0 0
102 65 72 63
4 3 2 4
8 6 5 7
1 3 1 5
18 11 28 13
6 3 4 0
3 5 4 3
1 23 19 8
2 3 0 1
0 0 0 0
1 4 3 2],
:Trt =>
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1],
:Base =>
[11, 11, 6, 8, 66, 27, 12, 52, 23, 10, 52, 33, 18, 42, 87, 50, 18, 111, 18,
20, 12, 9, 17, 28, 55, 9, 10, 47, 76, 38, 19, 10, 19, 24, 31, 14, 11, 67,
41, 7, 22, 13, 46, 36, 38, 7, 36, 11, 151, 22, 41, 32, 56, 24, 16, 22, 25,
13, 12],
:Age =>
[31, 30, 25, 36, 22, 29, 31, 42, 37, 28, 36, 24, 23, 36, 26, 26, 28, 31, 32,
21, 29, 21, 32, 25, 30, 40, 19, 22, 18, 32, 20, 30, 18, 24, 30, 35, 27, 20,
22, 28, 23, 40, 33, 21, 35, 25, 26, 25, 22, 32, 25, 35, 21, 41, 32, 26, 21,
36, 37],
:V4 => [0, 0, 0, 1]
)
epil[:N] = size(epil[:y], 1)
epil[:T] = size(epil[:y], 2)
epil[:logBase4] = log(epil[:Base] / 4)
epil[:BT] = epil[:logBase4] .* epil[:Trt]
epil[:logAge] = log(epil[:Age])
map(key > epil[Symbol(string(key, "bar"))] = mean(epil[key]),
[:logBase4, :Trt, :BT, :logAge, :V4])
## Model Specification
model = Model(
y = Stochastic(2,
(a0, alpha_Base, logBase4, logBase4bar, alpha_Trt, Trt, Trtbar,
alpha_BT, BT, BTbar, alpha_Age, logAge, logAgebar, alpha_V4, V4,
V4bar, b1, b, N, T) >
UnivariateDistribution[
begin
mu = exp(a0 + alpha_Base * (logBase4[i]  logBase4bar) +
alpha_Trt * (Trt[i]  Trtbar) + alpha_BT * (BT[i]  BTbar) +
alpha_Age * (logAge[i]  logAgebar) +
alpha_V4 * (V4[j]  V4bar) + b1[i] +
b[i, j])
Poisson(mu)
end
for i in 1:N, j in 1:T
],
false
),
b1 = Stochastic(1,
s2_b1 > Normal(0, sqrt(s2_b1)),
false
),
b = Stochastic(2,
s2_b > Normal(0, sqrt(s2_b)),
false
),
a0 = Stochastic(
() > Normal(0, 100),
false
),
alpha_Base = Stochastic(
() > Normal(0, 100)
),
alpha_Trt = Stochastic(
() > Normal(0, 100)
),
alpha_BT = Stochastic(
() > Normal(0, 100)
),
alpha_Age = Stochastic(
() > Normal(0, 100)
),
alpha_V4 = Stochastic(
() > Normal(0, 100)
),
alpha0 = Logical(
(a0, alpha_Base, logBase4bar, alpha_Trt, Trtbar, alpha_BT, BTbar,
alpha_Age, logAgebar, alpha_V4, V4bar) >
a0  alpha_Base * logBase4bar  alpha_Trt * Trtbar  alpha_BT * BTbar 
alpha_Age * logAgebar  alpha_V4 * V4bar
),
s2_b1 = Stochastic(
() > InverseGamma(0.001, 0.001)
),
s2_b = Stochastic(
() > InverseGamma(0.001, 0.001)
)
)
## Initial Values
inits = [
Dict(:y => epil[:y], :a0 => 0, :alpha_Base => 0, :alpha_Trt => 0,
:alpha_BT => 0, :alpha_Age => 0, :alpha_V4 => 0, :s2_b1 => 1,
:s2_b => 1, :b1 => zeros(epil[:N]), :b => zeros(epil[:N], epil[:T])),
Dict(:y => epil[:y], :a0 => 1, :alpha_Base => 1, :alpha_Trt => 1,
:alpha_BT => 1, :alpha_Age => 1, :alpha_V4 => 1, :s2_b1 => 10,
:s2_b => 10, :b1 => zeros(epil[:N]), :b => zeros(epil[:N], epil[:T]))
]
## Sampling Scheme
scheme = [AMWG([:a0, :alpha_Base, :alpha_Trt, :alpha_BT, :alpha_Age,
:alpha_V4], 0.1),
Slice(:b1, 0.5),
Slice(:b, 0.5),
Slice([:s2_b1, :s2_b], 1.0)]
setsamplers!(model, scheme)
## MCMC Simulations
sim = mcmc(model, epil, inits, 15000, burnin=2500, thin=2, chains=2)
describe(sim)
Iterations = 2502:15000
Thinning interval = 2
Chains = 1,2
Samples per chain = 6250
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
s2_b 0.13523750 0.031819272 0.00028460022 0.0025401820 156.91007
s2_b1 0.24911885 0.073166731 0.00065442313 0.0044942066 265.04599
alpha_V4 0.09287934 0.083666872 0.00074833925 0.0051087325 268.21356
alpha_Age 0.45830900 0.394536219 0.00352883922 0.0310012419 161.96291
alpha_BT 0.24217000 0.190566444 0.00170447809 0.0163585849 135.70673
alpha_Trt 0.75931393 0.397734236 0.00355744316 0.0337796459 138.63592
alpha_Base 0.91104974 0.135354470 0.00121064718 0.0111438503 147.52807
alpha0 1.35617079 1.313240197 0.01174597740 0.1021568814 165.25442
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
s2_b 0.07155807 0.112591385 0.13626498 0.158032604 0.193715882
s2_b1 0.13817484 0.197135457 0.23761145 0.289655043 0.422805121
alpha_V4 0.25504531 0.148157017 0.09313603 0.036681431 0.072099011
alpha_Age 0.19666987 0.176356196 0.41608686 0.696647899 1.305075377
alpha_BT 0.09025014 0.108102968 0.22656174 0.358354280 0.657804969
alpha_Trt 1.63682212 1.011390894 0.75653998 0.480870874 0.016113397
alpha_Base 0.66318818 0.817700638 0.90268210 0.997417378 1.200619714
alpha0 4.16888778 2.157932918 1.26343143 0.436226488 0.866195785
Blocker: Random Effects MetaAnalysis of Clinical Trials¶
An example from OpenBUGS [44] and Carlin [15] concerning a metaanalysis of 22 clinical trials to prevent mortality after myocardial infarction.
Events are modelled as
where is the number of control group events, out of , in study ; and is the number of treatment group events.
using Mamba
## Data
blocker = Dict{Symbol, Any}(
:rt =>
[3, 7, 5, 102, 28, 4, 98, 60, 25, 138, 64, 45, 9, 57, 25, 33, 28, 8, 6, 32,
27, 22],
:nt =>
[38, 114, 69, 1533, 355, 59, 945, 632, 278, 1916, 873, 263, 291, 858, 154,
207, 251, 151, 174, 209, 391, 680],
:rc =>
[3, 14, 11, 127, 27, 6, 152, 48, 37, 188, 52, 47, 16, 45, 31, 38, 12, 6, 3,
40, 43, 39],
:nc =>
[39, 116, 93, 1520, 365, 52, 939, 471, 282, 1921, 583, 266, 293, 883, 147,
213, 122, 154, 134, 218, 364, 674]
)
blocker[:N] = length(blocker[:rt])
## Model Specification
model = Model(
rc = Stochastic(1,
(mu, nc, N) >
begin
pc = invlogit(mu)
UnivariateDistribution[Binomial(nc[i], pc[i]) for i in 1:N]
end,
false
),
rt = Stochastic(1,
(mu, delta, nt, N) >
begin
pt = invlogit(mu + delta)
UnivariateDistribution[Binomial(nt[i], pt[i]) for i in 1:N]
end,
false
),
mu = Stochastic(1,
() > Normal(0, 1000),
false
),
delta = Stochastic(1,
(d, s2) > Normal(d, sqrt(s2)),
false
),
delta_new = Stochastic(
(d, s2) > Normal(d, sqrt(s2))
),
d = Stochastic(
() > Normal(0, 1000)
),
s2 = Stochastic(
() > InverseGamma(0.001, 0.001)
)
)
## Initial Values
inits = [
Dict(:rc => blocker[:rc], :rt => blocker[:rt], :d => 0, :delta_new => 0,
:s2 => 1, :mu => zeros(blocker[:N]), :delta => zeros(blocker[:N])),
Dict(:rc => blocker[:rc], :rt => blocker[:rt], :d => 2, :delta_new => 2,
:s2 => 10, :mu => fill(2, blocker[:N]), :delta => fill(2, blocker[:N]))
]
## Sampling Scheme
scheme = [AMWG(:mu, 0.1),
AMWG([:delta, :delta_new], 0.1),
Slice([:d, :s2], 1.0)]
setsamplers!(model, scheme)
## MCMC Simulations
sim = mcmc(model, blocker, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)
Iterations = 2502:10000
Thinning interval = 2
Chains = 1,2
Samples per chain = 3750
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
s2 0.01822186 0.021121265 0.00024388736 0.0014150714 222.78358
d 0.25563567 0.061841945 0.00071408927 0.0040205781 236.58613
delta_new 0.25005767 0.150325282 0.00173580684 0.0050219145 896.03592
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
s2 0.0006855452 0.0041648765 0.0107615659 0.024442084 0.07735715
d 0.3734122953 0.2959169814 0.2581848849 0.218341380 0.12842580
delta_new 0.5385405488 0.3279958446 0.2557849252 0.177588413 0.07986060
Oxford: Smooth Fit to LogOdds Ratios¶
An example from OpenBUGS [44] and Breslow and Clayton [10] concerning the association between death from childhood cancer and maternal exposure to Xrays, for subjects partitioned into 120 age and birthyear strata.
Deaths are modelled as
where is the number of deaths among unexposed subjects in stratum , is the number among exposed subjects, and is the stratumspecific birth year (relative to 1954).
using Mamba
## Data
oxford = Dict{Symbol, Any}(
:r1 =>
[3, 5, 2, 7, 7, 2, 5, 3, 5, 11, 6, 6, 11, 4, 4, 2, 8, 8, 6, 5, 15, 4, 9, 9,
4, 12, 8, 8, 6, 8, 12, 4, 7, 16, 12, 9, 4, 7, 8, 11, 5, 12, 8, 17, 9, 3, 2,
7, 6, 5, 11, 14, 13, 8, 6, 4, 8, 4, 8, 7, 15, 15, 9, 9, 5, 6, 3, 9, 12, 14,
16, 17, 8, 8, 9, 5, 9, 11, 6, 14, 21, 16, 6, 9, 8, 9, 8, 4, 11, 11, 6, 9,
4, 4, 9, 9, 10, 14, 6, 3, 4, 6, 10, 4, 3, 3, 10, 4, 10, 5, 4, 3, 13, 1, 7,
5, 7, 6, 3, 7],
:n1 =>
[28, 21, 32, 35, 35, 38, 30, 43, 49, 53, 31, 35, 46, 53, 61, 40, 29, 44, 52,
55, 61, 31, 48, 44, 42, 53, 56, 71, 43, 43, 43, 40, 44, 70, 75, 71, 37, 31,
42, 46, 47, 55, 63, 91, 43, 39, 35, 32, 53, 49, 75, 64, 69, 64, 49, 29, 40,
27, 48, 43, 61, 77, 55, 60, 46, 28, 33, 32, 46, 57, 56, 78, 58, 52, 31, 28,
46, 42, 45, 63, 71, 69, 43, 50, 31, 34, 54, 46, 58, 62, 52, 41, 34, 52, 63,
59, 88, 62, 47, 53, 57, 74, 68, 61, 45, 45, 62, 73, 53, 39, 45, 51, 55, 41,
53, 51, 42, 46, 54, 32],
:r0 =>
[0, 2, 2, 1, 2, 0, 1, 1, 1, 2, 4, 4, 2, 1, 7, 4, 3, 5, 3, 2, 4, 1, 4, 5, 2,
7, 5, 8, 2, 3, 5, 4, 1, 6, 5, 11, 5, 2, 5, 8, 5, 6, 6, 10, 7, 5, 5, 2, 8,
1, 13, 9, 11, 9, 4, 4, 8, 6, 8, 6, 8, 14, 6, 5, 5, 2, 4, 2, 9, 5, 6, 7, 5,
10, 3, 2, 1, 7, 9, 13, 9, 11, 4, 8, 2, 3, 7, 4, 7, 5, 6, 6, 5, 6, 9, 7, 7,
7, 4, 2, 3, 4, 10, 3, 4, 2, 10, 5, 4, 5, 4, 6, 5, 3, 2, 2, 4, 6, 4, 1],
:n0 =>
[28, 21, 32, 35, 35, 38, 30, 43, 49, 53, 31, 35, 46, 53, 61, 40, 29, 44, 52,
55, 61, 31, 48, 44, 42, 53, 56, 71, 43, 43, 43, 40, 44, 70, 75, 71, 37, 31,
42, 46, 47, 55, 63, 91, 43, 39, 35, 32, 53, 49, 75, 64, 69, 64, 49, 29, 40,
27, 48, 43, 61, 77, 55, 60, 46, 28, 33, 32, 46, 57, 56, 78, 58, 52, 31, 28,
46, 42, 45, 63, 71, 69, 43, 50, 31, 34, 54, 46, 58, 62, 52, 41, 34, 52, 63,
59, 88, 62, 47, 53, 57, 74, 68, 61, 45, 45, 62, 73, 53, 39, 45, 51, 55, 41,
53, 51, 42, 46, 54, 32],
:year =>
[10, 9, 9, 8, 8, 8, 7, 7, 7, 7, 6, 6, 6, 6, 6, 5, 5, 5,
5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 3, 3, 3, 2,
2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2,
2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6,
6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 9, 9, 10],
:K => 120
)
oxford[:K] = length(oxford[:r1])
## Model Specification
model = Model(
r0 = Stochastic(1,
(mu, n0, K) >
begin
p = invlogit(mu)
UnivariateDistribution[Binomial(n0[i], p[i]) for i in 1:K]
end,
false
),
r1 = Stochastic(1,
(mu, alpha, beta1, beta2, year, b, n1, K) >
UnivariateDistribution[
begin
p = invlogit(mu[i] + alpha + beta1 * year[i] +
beta2 * (year[i]^2  22.0) + b[i])
Binomial(n1[i], p)
end
for i in 1:K
],
false
),
b = Stochastic(1,
s2 > Normal(0, sqrt(s2)),
false
),
mu = Stochastic(1,
() > Normal(0, 1000),
false
),
alpha = Stochastic(
() > Normal(0, 1000)
),
beta1 = Stochastic(
() > Normal(0, 1000)
),
beta2 = Stochastic(
() > Normal(0, 1000)
),
s2 = Stochastic(
() > InverseGamma(0.001, 0.001)
)
)
## Initial Values
inits = [
Dict(:r0 => oxford[:r0], :r1 => oxford[:r1], :alpha => 0, :beta1 => 0,
:beta2 => 0, :s2 => 1, :b => zeros(oxford[:K]),
:mu => zeros(oxford[:K])),
Dict(:r0 => oxford[:r0], :r1 => oxford[:r1], :alpha => 1, :beta1 => 1,
:beta2 => 1, :s2 => 10, :b => zeros(oxford[:K]),
:mu => zeros(oxford[:K]))
]
## Sampling Scheme
scheme = [AMWG([:alpha, :beta1, :beta2], 1.0),
Slice(:s2, 1.0),
Slice(:mu, 1.0),
Slice(:b, 1.0)]
setsamplers!(model, scheme)
## MCMC Simulations
sim = mcmc(model, oxford, inits, 12500, burnin=2500, thin=2, chains=2)
describe(sim)
Iterations = 2502:12500
Thinning interval = 2
Chains = 1,2
Samples per chain = 5000
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
beta2 0.005477119 0.0035675748 0.00003567575 0.00033192987 115.519023
beta1 0.043336269 0.0161754258 0.00016175426 0.00133361554 147.112695
alpha 0.565784774 0.0630050896 0.00063005090 0.00468384860 180.944576
s2 0.026238992 0.0307989154 0.00030798915 0.00302056007 103.967091
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
beta2 0.0010499046 0.0028489198 0.0056500394 0.0077473623 0.013630865
beta1 0.0745152363 0.0543180318 0.0434425931 0.0321216097 0.009920787
alpha 0.4438257884 0.5238801187 0.5675039159 0.6051427125 0.695968063
s2 0.0007134423 0.0033352655 0.0146737037 0.0397132522 0.118202258
LSAT: Item Response¶
An example from OpenBUGS [44] and Boch and Lieberman [6] concerning a 5item multiple choice test (32 possible response patters) given to 1000 students.
Item responses are modelled as
where is an indicator for correct response by student to questions .
using Mamba
## Data
lsat = Dict{Symbol, Any}(
:culm =>
[3, 9, 11, 22, 23, 24, 27, 31, 32, 40, 40, 56, 56, 59, 61, 76, 86, 115, 129,
210, 213, 241, 256, 336, 352, 408, 429, 602, 613, 674, 702, 1000],
:response =>
[0 0 0 0 0
0 0 0 0 1
0 0 0 1 0
0 0 0 1 1
0 0 1 0 0
0 0 1 0 1
0 0 1 1 0
0 0 1 1 1
0 1 0 0 0
0 1 0 0 1
0 1 0 1 0
0 1 0 1 1
0 1 1 0 0
0 1 1 0 1
0 1 1 1 0
0 1 1 1 1
1 0 0 0 0
1 0 0 0 1
1 0 0 1 0
1 0 0 1 1
1 0 1 0 0
1 0 1 0 1
1 0 1 1 0
1 0 1 1 1
1 1 0 0 0
1 1 0 0 1
1 1 0 1 0
1 1 0 1 1
1 1 1 0 0
1 1 1 0 1
1 1 1 1 0
1 1 1 1 1],
:N => 1000
)
lsat[:R] = size(lsat[:response], 1)
lsat[:T] = size(lsat[:response], 2)
n = [lsat[:culm][1]; diff(lsat[:culm])]
idx = mapreduce(i > fill(i, n[i]), vcat, 1:length(n))
lsat[:r] = lsat[:response][idx, :]
## Model Specification
model = Model(
r = Stochastic(2,
(beta, theta, alpha, N, T) >
UnivariateDistribution[
begin
p = invlogit(beta * theta[i]  alpha[j])
Bernoulli(p)
end
for i in 1:N, j in 1:T
],
false
),
theta = Stochastic(1,
() > Normal(0, 1),
false
),
alpha = Stochastic(1,
() > Normal(0, 100),
false
),
a = Logical(1,
alpha > alpha  mean(alpha)
),
beta = Stochastic(
() > Truncated(Flat(), 0, Inf)
)
)
## Initial Values
inits = [
Dict(:r => lsat[:r], :alpha => zeros(lsat[:T]), :beta => 1,
:theta => zeros(lsat[:N])),
Dict(:r => lsat[:r], :alpha => ones(lsat[:T]), :beta => 2,
:theta => zeros(lsat[:N]))
]
## Sampling Scheme
scheme = [AMWG(:alpha, 0.1),
Slice(:beta, 1.0),
Slice(:theta, 0.5)]
setsamplers!(model, scheme)
## MCMC Simulations
sim = mcmc(model, lsat, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)
Iterations = 2502:10000
Thinning interval = 2
Chains = 1,2
Samples per chain = 3750
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
beta 0.80404469 0.072990202 0.00084281825 0.0067491110 116.95963
a[1] 1.26236241 0.104040037 0.00120135087 0.0025355922 1683.61261
a[2] 0.48004293 0.069256031 0.00079969977 0.0013701533 2554.91753
a[3] 1.24206491 0.068338749 0.00078910790 0.0018426724 1375.42777
a[4] 0.16982672 0.072942222 0.00084226422 0.0012659899 3319.68982
a[5] 0.62957215 0.086601562 0.00099998871 0.0018787409 2124.79816
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
beta 0.678005795 0.75195190 0.79754709 0.85100547 0.96030934
a[1] 1.471693755 1.33040793 1.25998457 1.19317801 1.06168159
a[2] 0.347262397 0.43161040 0.48023957 0.52718291 0.61527668
a[3] 1.106413529 1.19669095 1.24105794 1.28858225 1.37451854
a[4] 0.023253916 0.11970853 0.17099598 0.21998896 0.30858397
a[5] 0.799988061 0.68755932 0.63052599 0.57168504 0.46028931
Bones: Latent Trait Model for Multiple Ordered Categorical Responses¶
An example from OpenBUGS [44], Roche et al. [82], and Thissen [93] concerning skeletal age in 13 boys predicted from 34 radiograph indicators of skeletal maturity.
Skeletal ages are modelled as
where is a discriminability parameter for indicator , is a threshold parameter, and is the cumulative probability that boy with skeletal age is assigned a more mature grade than .
using Mamba
## Data
bones = Dict{Symbol, Any}(
:gamma => reshape(
[ 0.7425, NaN, NaN, NaN, 10.2670, NaN, NaN, NaN,
10.5215, NaN, NaN, NaN, 9.3877, NaN, NaN, NaN,
0.2593, NaN, NaN, NaN, 0.5998, NaN, NaN, NaN,
10.5891, NaN, NaN, NaN, 6.6701, NaN, NaN, NaN,
8.8921, NaN, NaN, NaN, 12.4275, NaN, NaN, NaN,
12.4788, NaN, NaN, NaN, 13.7778, NaN, NaN, NaN,
5.8374, NaN, NaN, NaN, 6.9485, NaN, NaN, NaN,
13.7184, NaN, NaN, NaN, 14.3476, NaN, NaN, NaN,
4.8066, NaN, NaN, NaN, 9.1037, NaN, NaN, NaN,
10.7483, NaN, NaN, NaN, 0.3887, 1.0153, NaN, NaN,
3.2573, 7.0421, NaN, NaN, 11.6273, 14.4242, NaN, NaN,
15.8842, 17.4685, NaN, NaN, 14.8926, 16.7409, NaN, NaN,
15.5487, 16.8720, NaN, NaN, 15.4091, 17.0061, NaN, NaN,
3.9216, 5.2099, NaN, NaN, 15.4750, 16.9406, 17.4944, NaN,
0.4927, 1.3556, 2.3016, 3.2535, 1.3059, 1.8793, 2.4970, 3.2306,
1.5012, 1.8902, 2.3689, 2.9495, 0.8021, 2.3873, 3.9525, 5.3198,
5.0022, 6.3704, 8.2832, 10.4988, 4.0168, 5.1537, 7.1053, 10.3038],
4, 34)',
:delta =>
[2.9541, 0.6603, 0.7965, 1.0495, 5.7874, 3.8376, 0.6324,
0.8272, 0.6968, 0.8747, 0.8136, 0.8246, 0.6711, 0.978,
1.1528, 1.6923, 1.0331, 0.5381, 1.0688, 8.1123, 0.9974,
1.2656, 1.1802, 1.368, 1.5435, 1.5006, 1.6766, 1.4297,
3.385, 3.3085, 3.4007, 2.0906, 1.0954, 1.5329],
:ncat =>
[2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
3, 3, 3, 3, 3, 3, 3, 3, 4, 5, 5, 5, 5, 5, 5],
:grade => reshape(
[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,1,1,1,1,1,1,1,1,2,1,1,2,1,1,
2,1,1,1,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,3,1,1,1,1,1,1,1,1,3,1,1,2,1,1,
2,1,1,1,2,2,1,1,1,1,1,1,1,1,1,1,1,1,1,3,1,1,1,1,1,1,1,1,4,3,3,3,1,1,
2,1,1,1,2,2,1,1,1,1,1,1,NaN,1,1,1,1,1,1,3,1,1,1,1,1,1,1,1,4,5,4,3,1,1,
2,1,1,1,2,2,1,1,2,1,1,1,1,1,1,1,2,1,1,3,2,1,1,1,1,1,3,1,5,5,5,4,2,3,
2,1,1,1,2,2,1,2,1,1,1,1,1,2,1,1,2,NaN,1,3,2,1,1,1,1,1,3,1,5,5,5,5,3,3,
2,1,1,1,2,2,1,1,1,NaN,NaN,1,1,1,1,1,2,NaN,1,3,3,1,1,1,1,1,3,1,5,5,5,5,3,3,
2,1,2,2,2,2,2,2,1,NaN,NaN,1,2,2,1,1,2,2,1,3,2,1,1,1,1,1,3,1,5,5,5,5,3,4,
2,1,1,2,2,2,2,2,2,1,1,1,2,1,1,1,2,1,1,3,3,1,1,1,1,1,3,1,5,5,5,5,4,4,
2,1,2,2,2,2,2,2,2,1,1,1,2,2,2,1,2,NaN,2,3,3,1,1,1,1,1,3,1,5,5,5,5,5,5,
2,1,NaN,2,2,2,NaN,2,2,1,NaN,NaN,2,2,NaN,NaN,2,1,2,3,3,NaN,1,NaN,1,1,3,1,5,5,5,5,5,5,
2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,1,NaN,2,1,3,2,5,5,5,5,5,5,
2,2,2,2,2,2,2,2,2,2,NaN,2,2,2,2,2,2,2,2,3,3,3,NaN,2,NaN,2,3,4,5,5,5,5,5,5],
34, 13)',
:nChild => 13,
:nInd => 34
)
## Model Specification
model = Model(
grade = Stochastic(2,
(ncat, delta, theta, gamma, nChild, nInd) >
begin
p = Array{Float64}(5)
UnivariateDistribution[
begin
n = ncat[j]
p[1] = 1.0
for k in 1:(n  1)
Q = invlogit(delta[j] * (theta[i]  gamma[j, k]))
p[k] = Q
p[k + 1] = Q
end
Categorical(p[1:n])
end
for i in 1:nChild, j in 1:nInd
]
end,
false
),
theta = Stochastic(1,
() > Normal(0, 100)
)
)
## Initial Values
inits = [
Dict(:grade => bones[:grade], :theta => [0.5,1,2,3,5,6,7,8,9,12,13,16,18]),
Dict(:grade => bones[:grade], :theta => [1,2,3,4,5,6,7,8,9,10,11,12,13])
]
## Sampling Scheme
scheme = [MISS(:grade),
AMWG(:theta, 0.1)]
setsamplers!(model, scheme)
## MCMC Simulations
sim = mcmc(model, bones, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)
Iterations = 2502:10000
Thinning interval = 2
Chains = 1,2
Samples per chain = 3750
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
theta[1] 0.32603385 0.20640874 0.0023834028 0.0039448110 2737.81276
theta[2] 1.37861692 0.25824308 0.0029819342 0.0058663965 1937.82503
theta[3] 2.35227822 0.27998526 0.0032329913 0.0067161153 1737.93707
theta[4] 2.90165730 0.29713320 0.0034309987 0.0078730921 1424.33353
theta[5] 5.54427283 0.50242324 0.0058014839 0.0169090038 882.88350
theta[6] 6.70804782 0.57206890 0.0066056827 0.0221532973 666.83738
theta[7] 6.49138381 0.60154625 0.0069460578 0.0219158412 753.39330
theta[8] 8.93701249 0.73636136 0.0085027686 0.0336199950 479.71875
theta[9] 9.03585289 0.65172497 0.0075254717 0.0233182299 781.15561
theta[10] 11.93125529 0.69360918 0.0080091090 0.0282955741 600.88678
theta[11] 11.53686992 0.92271657 0.0106546132 0.0493587234 349.46912
theta[12] 15.81482824 0.54261736 0.0062656056 0.0210976666 661.48275
theta[13] 16.93028146 0.72458739 0.0083668145 0.0323302069 502.30161
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
theta[1] 0.11215555 0.19557824 0.33881555 0.45840506 0.7174563
theta[2] 0.91705346 1.19969433 1.36116575 1.53751273 1.9466119
theta[3] 1.78287586 2.17136780 2.35623189 2.53035766 2.9211580
theta[4] 2.32940825 2.69621746 2.89121336 3.10758151 3.4945343
theta[5] 4.59142954 5.20543314 5.53392246 5.86525435 6.5867245
theta[6] 5.56649066 6.30983797 6.70666338 7.09569168 7.8229872
theta[7] 5.38663728 6.07628064 6.46533033 6.88636840 7.7051374
theta[8] 7.47304526 8.43125608 8.96072241 9.45344704 10.2856733
theta[9] 7.80477915 8.60559136 9.01498109 9.46962522 10.3024722
theta[10] 10.64129157 11.48379528 11.89611699 12.37737647 13.3873043
theta[11] 9.83558611 10.88717498 11.49029895 12.15757004 13.4263451
theta[12] 14.79250437 15.45889470 15.79840132 16.15824313 16.9593310
theta[13] 15.61843069 16.42289725 16.90719268 17.41900248 18.3895761
Inhalers: Ordered Categorical Data¶
An example from OpenBUGS [44] and Ezzet and Whitehead [25] concerning a twotreatment, twoperiod crossover trial comparing salbutamol inhalation devices in 286 asthma patients.
Treatment responses are modelled as
where is a 4point ordinal rating of the device used by patient , and is the cumulative probability of the rating in treatment period being worse than category .
using Mamba
## Data
inhalers = Dict{Symbol, Any}(
:pattern =>
[1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4
1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4]',
:Ncum =>
[ 59 157 173 175 186 253 270 271 271 278 280 281 282 285 285 286
122 170 173 175 226 268 270 271 278 280 281 281 284 285 286 286]',
:treat =>
[ 1 1
1 1],
:period =>
[1 1
1 1],
:carry =>
[0 1
0 1],
:N => 286,
:T => 2,
:G => 2,
:Npattern => 16,
:Ncut => 3
)
inhalers[:group] = Array{Int}(inhalers[:N])
inhalers[:response] = Array{Int}(inhalers[:N], inhalers[:T])
i = 1
for k in 1:inhalers[:Npattern], g in 1:inhalers[:G]
while i <= inhalers[:Ncum][k, g]
inhalers[:group][i] = g
for t in 1:inhalers[:T]
inhalers[:response][i, t] = inhalers[:pattern][k, t]
end
i += 1
end
end
## Model Specification
model = Model(
response = Stochastic(2,
(a1, a2, a3, mu, group, b, N, T) >
begin
a = Float64[a1, a2, a3]
UnivariateDistribution[
begin
eta = mu[group[i], t] + b[i]
p = ones(4)
for j in 1:3
Q = invlogit((a[j] + eta))
p[j] = Q
p[j + 1] = Q
end
Categorical(p)
end
for i in 1:N, t in 1:T
]
end,
false
),
mu = Logical(2,
(beta, treat, pi, period, kappa, carry, G, T) >
[ beta * treat[g, t] / 2 + pi * period[g, t] / 2 + kappa * carry[g, t]
for g in 1:G, t in 1:T ],
false
),
b = Stochastic(1,
s2 > Normal(0, sqrt(s2)),
false
),
a1 = Stochastic(
a2 > Truncated(Flat(), 1000, a2)
),
a2 = Stochastic(
a3 > Truncated(Flat(), 1000, a3)
),
a3 = Stochastic(
() > Truncated(Flat(), 1000, 1000)
),
beta = Stochastic(
() > Normal(0, 1000)
),
pi = Stochastic(
() > Normal(0, 1000)
),
kappa = Stochastic(
() > Normal(0, 1000)
),
s2 = Stochastic(
() > InverseGamma(0.001, 0.001)
)
)
## Initial Values
inits = [
Dict(:response => inhalers[:response], :beta => 0, :pi => 0, :kappa => 0,
:a1 => 2, :a2 => 3, :a3 => 4, :s2 => 1, :b => zeros(inhalers[:N])),
Dict(:response => inhalers[:response], :beta => 1, :pi => 1, :kappa => 0,
:a1 => 3, :a2 => 4, :a3 => 5, :s2 => 10, :b => zeros(inhalers[:N]))
]
## Sampling Scheme
scheme = [AMWG(:b, 0.1),
Slice([:a1, :a2, :a3], 2.0),
Slice([:beta, :pi, :kappa, :s2], 1.0, Univariate)]
setsamplers!(model, scheme)
## MCMC Simulations
sim = mcmc(model, inhalers, inits, 5000, burnin=1000, thin=2, chains=2)
describe(sim)
Mice: Weibull Regression¶
An example from OpenBUGS [44], Grieve [42], and Dellaportas and Smith [22] concerning time to death or censoring among four groups of 20 mice each.
Time to events are modelled as
where is the time of death for mouse , and is a vector of covariates.
using Mamba
## Data
mice = Dict{Symbol, Any}(
:t =>
[12 1 21 25 11 26 27 30 13 12 21 20 23 25 23 29 35 NaN 31 36
32 27 23 12 18 NaN NaN 38 29 30 NaN 32 NaN NaN NaN NaN 25 30 37 27
22 26 NaN 28 19 15 12 35 35 10 22 18 NaN 12 NaN NaN 31 24 37 29
27 18 22 13 18 29 28 NaN 16 22 26 19 NaN NaN 17 28 26 12 17 26],
:tcensor =>
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 40 0 0
0 0 0 0 0 40 40 0 0 0 40 0 40 40 40 40 0 0 0 0
0 0 10 0 0 0 0 0 0 0 0 0 24 0 40 40 0 0 0 0
0 0 0 0 0 0 0 20 0 0 0 0 29 10 0 0 0 0 0 0]
)
mice[:M] = size(mice[:t], 1)
mice[:N] = size(mice[:t], 2)
## Model Specification
model = Model(
t = Stochastic(2,
(r, beta, tcensor, M, N) >
UnivariateDistribution[
begin
lambda = exp(beta[i] / r)
0 < lambda < Inf ?
Truncated(Weibull(r, lambda), tcensor[i, j], Inf) :
Uniform(0, Inf)
end
for i in 1:M, j in 1:N
],
false
),
r = Stochastic(
() > Exponential(1000)
),
beta = Stochastic(1,
() > Normal(0, 10),
false
),
median = Logical(1,
(beta, r) > exp(beta / r) * log(2)^(1 / r)
),
veh_control = Logical(
beta > beta[2]  beta[1]
),
test_sub = Logical(
beta > beta[3]  beta[1]
),
pos_control = Logical(
beta > beta[4]  beta[1]
)
)
## Initial Values
inits = [
Dict(:t => mice[:t], :beta => fill(1, mice[:M]), :r => 1.0),
Dict(:t => mice[:t], :beta => fill(2, mice[:M]), :r => 1.0)
]
## Sampling Scheme
scheme = [MISS(:t),
Slice(:beta, 1.0, Univariate),
Slice(:r, 0.25)]
setsamplers!(model, scheme)
## MCMC Simulations
sim = mcmc(model, mice, inits, 20000, burnin=2500, thin=2, chains=2)
describe(sim)
Leuk: Cox Regression¶
An example from OpenBUGS [44] and Ezzet and Whitehead [27] concerning survival in 42 leukemia patients treated with 6mercaptopurine or placebo.
Times to death are modelled using the Bayesian Cox proportional hazards model, formulated by Clayton [17] as
where is a counting process increment in time interval for patient ; is an indicator of whether the patient is observed at time ; is a vector of covariates; and is the increment in the integrated baseline hazard function during .
using Mamba
## Data
leuk = Dict{Symbol, Any}(
:t_obs =>
[1, 1, 2, 2, 3, 4, 4, 5, 5, 8, 8, 8, 8, 11, 11, 12, 12, 15, 17, 22, 23, 6,
6, 6, 6, 7, 9, 10, 10, 11, 13, 16, 17, 19, 20, 22, 23, 25, 32, 32, 34, 35],
:fail =>
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0],
:Z =>
[0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
0.5, 0.5],
:t => [1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 13, 15, 16, 17, 22, 23, 35]
)
leuk[:N] = N = length(leuk[:t_obs])
leuk[:T] = T = length(leuk[:t])  1
leuk[:Y] = Array{Int}(N, T)
leuk[:dN] = Array{Int}(N, T)
for i in 1:N
for j in 1:T
leuk[:dN][i, j] = leuk[:fail][i] * (leuk[:t_obs][i] == leuk[:t][j])
leuk[:Y][i, j] = Int(leuk[:t_obs][i] >= leuk[:t][j])
end
end
leuk[:c] = 0.001
leuk[:r] = 0.1
## Model Specification
model = Model(
dN = Stochastic(2,
(Y, beta, Z, dL0, N, T) >
UnivariateDistribution[
Y[i, j] > 0 ? Poisson(exp(beta * Z[i]) * dL0[j]) : Flat()
for i in 1:N, j in 1:T
],
false
),
mu = Logical(1,
(c, r, t) > c * r * (t[2:end]  t[1:(end  1)]),
false
),
dL0 = Stochastic(1,
(mu, c, T) > UnivariateDistribution[Gamma(mu[j], 1 / c) for j in 1:T],
false
),
beta = Stochastic(
() > Normal(0, 1000)
),
S0 = Logical(1,
dL0 > exp(cumsum(dL0)),
false
),
S_treat = Logical(1,
(S0, beta) > S0.^exp(0.5 * beta)
),
S_placebo = Logical(1,
(S0, beta) > S0.^exp(0.5 * beta)
)
)
## Initial Values
inits = [
Dict(:dN => leuk[:dN], :beta => 0, :dL0 => fill(1, leuk[:T])),
Dict(:dN => leuk[:dN], :beta => 1, :dL0 => fill(2, leuk[:T]))
]
## Sampling Scheme
scheme = [AMWG(:dL0, 0.1),
Slice(:beta, 3.0)]
setsamplers!(model, scheme)
## MCMC Simulations
sim = mcmc(model, leuk, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)
Iterations = 2502:10000
Thinning interval = 2
Chains = 1,2
Samples per chain = 3750
Empirical Posterior Estimates:
Mean SD Naive SE MCSE ESS
beta 1.552064427 0.424977799 0.00490722093 0.0121343026 1226.59967
S_treat[1] 0.983021837 0.014032858 0.00016203749 0.0006444873 474.09293
S_treat[2] 0.966246426 0.020819923 0.00024040776 0.0009762155 454.84854
S_treat[3] 0.956207131 0.024534430 0.00028329919 0.0011874261 426.91239
S_treat[4] 0.936432950 0.031636973 0.00036531230 0.0015140463 436.62796
S_treat[5] 0.913982241 0.037982886 0.00043858859 0.0017354084 479.04081
S_treat[6] 0.879372323 0.047653044 0.00055024996 0.0021437962 494.09936
S_treat[7] 0.867528817 0.051602796 0.00059585777 0.0024545959 441.96358
S_treat[8] 0.821750962 0.064689997 0.00074697574 0.0030256573 457.12479
S_treat[9] 0.805557685 0.068612901 0.00079227353 0.0032183437 454.51343
S_treat[10] 0.771831536 0.076942474 0.00088845517 0.0035776788 462.51903
S_treat[11] 0.735657019 0.085534730 0.00098766999 0.0040566359 444.58306
S_treat[12] 0.712600209 0.089598298 0.00103459203 0.0040292944 494.47179
S_treat[13] 0.691135357 0.094685927 0.00109333891 0.0044879053 445.12656
S_treat[14] 0.664423491 0.098857036 0.00114150273 0.0046673584 448.61405
S_treat[15] 0.636423003 0.102857125 0.00118769178 0.0047872478 461.63310
S_treat[16] 0.565616003 0.112893079 0.00130357699 0.0049341458 523.49276
S_treat[17] 0.471034334 0.120102602 0.00138682539 0.0050776645 559.47003
S_placebo[1] 0.927789861 0.050170096 0.00057931437 0.0023120443 470.86627
S_placebo[2] 0.859431833 0.067291385 0.00077701399 0.0030529117 485.83685
S_placebo[3] 0.820804570 0.074342778 0.00085843646 0.0034349005 468.43490
S_placebo[4] 0.748344200 0.085488930 0.00098714114 0.0038433955 494.75434
S_placebo[5] 0.671088723 0.090818803 0.00104868521 0.0036921538 605.05098
S_placebo[6] 0.564655613 0.097783399 0.00112910544 0.0038182815 655.83466
S_placebo[7] 0.532173306 0.099628193 0.00115040728 0.0042364979 553.03234
S_placebo[8] 0.418874416 0.097451300 0.00112527068 0.0038432054 642.96611
S_placebo[9] 0.383210551 0.095687796 0.00110490749 0.0036248738 696.83078
S_placebo[10] 0.317122087 0.091201496 0.00105310416 0.0031779306 823.59766
S_placebo[11] 0.256739191 0.086400289 0.00099766461 0.0029358329 866.09938
S_placebo[12] 0.223014153 0.082260904 0.00094986710 0.0026813211 941.21603
S_placebo[13] 0.195547225 0.079430544 0.00091718492 0.0029106845 744.70591
S_placebo[14] 0.164855445 0.074276591 0.00085767220 0.0026970503 758.44802
S_placebo[15] 0.137032754 0.068763778 0.00079401571 0.0026184621 689.64703
S_placebo[16] 0.083796000 0.054748991 0.00063218689 0.0020660044 702.24679
S_placebo[17] 0.040920336 0.037737842 0.00043575906 0.0012818385 866.73735
Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
beta 0.7524244645 1.259092877 1.540715015 1.829373993 2.41862572
S_treat[1] 0.9462141591 0.977211258 0.986955387 0.992864719 0.99836877
S_treat[2] 0.9150951603 0.955732194 0.970764331 0.981528195 0.99312428
S_treat[3] 0.8974993157 0.943194903 0.960671242 0.974341293 0.98958205
S_treat[4] 0.8599068641 0.919023530 0.941699242 0.959578177 0.98196246
S_treat[5] 0.8239755135 0.892499135 0.919841738 0.941506275 0.97166964
S_treat[6] 0.7738755125 0.850772834 0.885177533 0.913944652 0.95521105
S_treat[7] 0.7531822416 0.836159276 0.873267366 0.904856224 0.95037821
S_treat[8] 0.6801713643 0.780967298 0.828861810 0.868789540 0.92883388
S_treat[9] 0.6555328459 0.762836917 0.813277192 0.854664506 0.91930265
S_treat[10] 0.6026848684 0.723608005 0.779448002 0.826762737 0.90169913
S_treat[11] 0.5503082957 0.681881638 0.743447395 0.797449392 0.87990722
S_treat[12] 0.5233565529 0.654970518 0.720927823 0.777608202 0.86466162
S_treat[13] 0.4922868942 0.631085044 0.699537585 0.758920006 0.85262814
S_treat[14] 0.4597297883 0.598392933 0.670665550 0.735415731 0.83825032
S_treat[15] 0.4269750729 0.566537132 0.641228870 0.709931594 0.82208295
S_treat[16] 0.3438735907 0.487880437 0.566936069 0.644918756 0.77749139
S_treat[17] 0.2466734959 0.383395600 0.469857377 0.555064846 0.70280881
S_placebo[1] 0.7997904171 0.903185030 0.939595572 0.964924422 0.99089079
S_placebo[2] 0.7055958749 0.819313313 0.869352795 0.908385797 0.96318697
S_placebo[3] 0.6563988800 0.773124226 0.828159645 0.876023887 0.94337131
S_placebo[4] 0.5660793194 0.692444549 0.754818343 <