Getting Started
This tutorial takes you from a fresh Julia installation to a fitted Bayesian model. By the end you will have written a hierarchical growth-curve model in BUGS notation, drawn posterior samples from it, and read off posterior summaries that you can check against published results.
Setup
If you do not have Julia yet, download it from julialang.org/downloads; any recent version will do. Then start Julia and install the packages this tutorial uses:
using Pkg
Pkg.add(["JuliaBUGS", "AbstractMCMC", "AdvancedHMC", "ADTypes", "Mooncake", "FlexiChains"])JuliaBUGS compiles the model. The other packages do the sampling and bookkeeping: AdvancedHMC provides the NUTS sampler, AbstractMCMC runs it, ADTypes and Mooncake supply the gradients that NUTS needs, and FlexiChains stores the posterior draws. Load them all:
using JuliaBUGS
using AbstractMCMC, AdvancedHMC, ADTypes, Mooncake, FlexiChainsThe model
We will fit Rats, the first example in the classic BUGS examples (Volume 1), taken from Gelfand et al. (1990). Thirty young rats were weighed weekly for five weeks, so we have weights $Y_{ij}$ for rat $i = 1, \dots, 30$ at ages $x_j = 8, 15, 22, 29, 36$ days. Each rat grows roughly along a straight line, but rats differ, so each gets its own intercept $\alpha_i$ and slope $\beta_i$, and those in turn are drawn from population-level distributions — a normal hierarchical model:
\[\begin{aligned} Y_{ij} &\sim \text{Normal}\!\left(\alpha_i + \beta_i (x_j - \bar{x}),\ \tau_c\right) \\ \alpha_i &\sim \text{Normal}(\alpha_c,\ \alpha_\tau) \\ \beta_i &\sim \text{Normal}(\beta_c,\ \beta_\tau) \end{aligned}\]
Here $\bar{x} = 22$ is the mean age (centering the ages reduces correlation between intercepts and slopes), and — following the BUGS convention — normal distributions are written in terms of a precision (1/variance), not a variance.
In JuliaBUGS you write this model with the @bugs macro:
rats = @bugs begin
for i in 1:N
for j in 1:T
Y[i, j] ~ dnorm(mu[i, j], tau_c)
mu[i, j] = alpha[i] + beta[i] * (x[j] - xbar)
end
alpha[i] ~ dnorm(alpha_c, alpha_tau)
beta[i] ~ dnorm(beta_c, beta_tau)
end
tau_c ~ dgamma(0.001, 0.001)
sigma = 1 / sqrt(tau_c)
alpha_c ~ dnorm(0.0, 1.0e-6)
alpha_tau ~ dgamma(0.001, 0.001)
beta_c ~ dnorm(0.0, 1.0e-6)
beta_tau ~ dgamma(0.001, 0.001)
alpha0 = alpha_c - xbar * beta_c
endIf you have used WinBUGS, OpenBUGS, or JAGS, this should look familiar. Reading it line by line:
~means "is distributed as":Y[i, j] ~ dnorm(mu[i, j], tau_c)says the weight is normal with meanmu[i, j]and precisiontau_c. Distributions keep their BUGS names (dnorm,dgamma, and so on).=defines a deterministic quantity (BUGS uses<-for this):sigma = 1 / sqrt(tau_c)is the residual standard deviation, computed from the precision, andalpha0is the population intercept extrapolated back to birth (age zero).forloops express repetition over rats and over measurement times — the "plates" of the model. The loop boundsNandTwill come from the data.
The vague dgamma(0.001, 0.001) and dnorm(0.0, 1.0e-6) priors are the standard noninformative choices from the original example. The result, rats, is a model definition: it is not tied to any data yet.
Data
The data is just a NamedTuple whose names match the variables the model expects — here the ages x, their mean xbar, the counts N and T, and the 30×5 matrix of weights Y:
data = (
x = [8.0, 15.0, 22.0, 29.0, 36.0],
xbar = 22,
N = 30,
T = 5,
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]
)If your data lives in the R-style list() format used by the classic BUGS systems, it translates directly to a NamedTuple; see Coming from WinBUGS, OpenBUGS, and JAGS for the details.
Fit the model
Compiling the model is a single call: apply the definition to the data. Two of the three lines below are the essential steps — that call, and attaching the automatic-differentiation backend that NUTS needs. The line in between is an optional performance setting, explained in the note that follows:
model = rats(data)
model = JuliaBUGS.set_evaluation_mode(
model, JuliaBUGS.UseGeneratedLogDensityFunction()
)
model = JuliaBUGS.BUGSModelWithGradient(model, AutoMooncake(; config=nothing))rats(data) combines the model definition with the data and returns a compiled model, ready for inference. BUGSModelWithGradient then wraps that model with Mooncake automatic differentiation: NUTS works by following the gradient of the log posterior, so any gradient-based sampler needs this step.
Out of the box, JuliaBUGS computes a model's log density by walking the model's graph — that works for every model with no extra setup, and you can omit this line entirely. UseGeneratedLogDensityFunction() instead generates and compiles a dedicated Julia function for this model's log density. NUTS will evaluate the log density many thousands of times below, so for a model like Rats this makes sampling substantially faster. Its requirements — the default parameterization produced by compilation and a compatible AD backend such as Mooncake — are already met by the code above, and if JuliaBUGS cannot generate such a function for a model, it warns and keeps the default evaluator. See Evaluation Modes and Automatic Differentiation for the options and tradeoffs.
Now draw posterior samples with NUTS, the standard gradient-based MCMC sampler:
chain = AbstractMCMC.sample(
model, NUTS(0.8), 3000;
chain_type = VNChain,
n_adapts = 1000,
discard_initial = 1000,
)[ Info: Found initial step size 0.00078125This runs 4000 iterations. The first 1000 are used to tune the sampler and are discarded (n_adapts and discard_initial), and the 3000 after them are the posterior draws. chain_type = VNChain collects the draws into a chain object keyed by variable name. Expect the run to take a few minutes.
Read the results
summarystats prints the familiar table of posterior means, standard deviations, and convergence diagnostics for every quantity in the model:
summarystats(chain)╭─FlexiSummary (9 statistics) ─────────────────────────────────────────────────╮
│ iter collapsed │
│ chain collapsed │
│ ↓ stat = [mean, std, mcse, ess_bulk, ess_tail, rhat, q5, q50, q95] │
│ │
│ Parameters (67) ── VarName │
│ Float64 beta_tau, beta_c, alpha_tau, alpha_c, tau_c, beta[1], alpha[1] │
│ beta[2], …, alpha[28], beta[29], alpha[29], beta[30], alpha[30] │
│ alpha0, sigma │
│ │
│ Extras (13) │
│ Float64 lp, n_steps, is_accept, acceptance_rate, log_density │
│ hamiltonian_energy, hamiltonian_energy_error │
│ max_hamiltonian_energy_error, tree_depth, numerical_error │
│ step_size, nom_step_size, is_adapt │
│ │
│ Summary │
│ param mean std mcse … q5 q50 q95 │
│ beta_tau 4.1341 1.5785 0.0438 2.1954 3.8798 6.9072 │
│ beta_c 6.1851 0.1070 0.0022 6.0114 6.1838 6.3588 │
│ alpha_tau 0.0050 0.0014 0.0000 0.0030 0.0048 0.0074 │
│ alpha_c 242.5799 2.7605 0.0511 237.9989 242.6104 247.0876 │
│ tau_c 0.0273 0.0040 0.0001 … 0.0210 0.0272 0.0341 │
│ beta[1] 6.0710 0.2433 0.0042 5.6845 6.0691 6.4735 │
│ alpha[1] 239.9339 2.7070 0.0448 235.5233 239.8924 244.3822 │
│ beta[2] 7.0500 0.2540 0.0051 6.6181 7.0571 7.4651 │
│ alpha[2] 247.8240 2.6208 0.0386 243.5329 247.8112 252.1417 │
│ ⋮ ⋱ │
│ alpha[27] 254.2840 2.7152 0.0396 249.9098 254.3209 258.7936 │
│ beta[28] 5.8448 0.2517 0.0042 … 5.4421 5.8463 6.2593 │
│ alpha[28] 242.9489 2.7108 0.0455 238.5560 242.8934 247.4283 │
│ beta[29] 5.6724 0.2463 0.0044 5.2684 5.6705 6.0781 │
│ alpha[29] 217.9382 2.7475 0.0471 213.5799 217.9012 222.4586 │
│ beta[30] 6.1262 0.2390 0.0038 5.7385 6.1273 6.5153 │
│ alpha[30] 241.4604 2.6904 0.0467 … 236.9753 241.4897 245.9332 │
│ alpha0 106.5075 3.5901 0.0700 100.6062 106.5107 112.4450 │
│ sigma 6.1013 0.4587 0.0120 5.4118 6.0663 6.9042 │
╰──────────────────────────────────────────────────────────────────────────────╯The table has a row for each of the 30 intercepts alpha[i] and slopes beta[i], but the scientific questions concern the population-level rows:
beta_c, the average growth rate: about 6.2 grams per day (published value 6.186),alpha0, the average weight extrapolated back to birth: about 106.6 grams (published value 106.6),sigma, the residual standard deviation of the weight measurements: about 6.1 grams (published value 6.093).
Your numbers will not match the published values to every digit, and they will change slightly each time you run the sampler: posterior means estimated from 3000 draws carry a small Monte Carlo error (reported in the mcse column), so agreement to within that error is exactly what success looks like. As a quick health check, rhat should be very close to 1 for every row.
That is the whole workflow: write the model with @bugs, put the data in a NamedTuple, compile it, attach a gradient backend when the sampler needs one, sample to fit, and summarystats to read the results.
Where next
- The Example Gallery walks through more classic BUGS models, ready to run.
- The Seeds example treats a random-effects logistic regression in more depth, including supplying explicit initial values.
- Initial Values explains named, partial, array-valued, and sampler-specific initialization.
- Coming from WinBUGS, OpenBUGS, and JAGS maps your existing workflow — R
list()data, initial values, CODA summaries — onto JuliaBUGS. - DoodleBUGS lets you build JuliaBUGS models by drawing the graph, in the spirit of DoodleBUGS from WinBUGS.