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 3000 iterations; the first 1000 are used to tune the sampler and are discarded (n_adapts and discard_initial), leaving 2000 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[2], beta[3], alpha[3], beta[4], alpha[4], beta[5], │
│ alpha[5], beta[6], alpha[6], beta[7], alpha[7], beta[8], alpha[8], │
│ beta[9], alpha[9], beta[10], alpha[10], beta[11], alpha[11], │
│ beta[12], alpha[12], beta[13], alpha[13], beta[14], alpha[14], │
│ beta[15], alpha[15], beta[16], alpha[16], beta[17], alpha[17], │
│ beta[18], alpha[18], beta[19], alpha[19], beta[20], alpha[20], │
│ beta[21], alpha[21], beta[22], alpha[22], beta[23], alpha[23], │
│ beta[24], alpha[24], beta[25], alpha[25], beta[26], alpha[26], │
│ beta[27], alpha[27], beta[28], 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 ess_bulk ess_tail rhat … │
│ beta_tau 4.1600 1.5652 0.0428 1459.6478 1706.9372 0.9997 … │
│ beta_c 6.1861 0.1073 0.0020 2806.8557 2438.1797 0.9998 … │
│ alpha_tau 0.0050 0.0013 0.0000 3260.1726 2260.9901 0.9998 … │
│ alpha_c 242.7066 2.7080 0.0483 3167.8783 2013.6230 1.0001 … │
│ tau_c 0.0276 0.0041 0.0001 1546.6512 1804.5013 1.0002 … │
│ beta[1] 6.0634 0.2386 0.0040 3484.0411 2124.8881 1.0001 … │
│ alpha[1] 239.9379 2.6932 0.0452 3561.4566 2320.8877 1.0002 … │
│ beta[2] 7.0469 0.2603 0.0051 2618.9281 1765.2968 1.0007 … │
│ alpha[2] 247.8069 2.6443 0.0460 3319.8035 1919.0614 0.9997 … │
│ beta[3] 6.4832 0.2437 0.0039 3966.5561 2351.3502 1.0006 … │
│ alpha[3] 252.4140 2.6287 0.0442 3552.5734 2057.0314 1.0018 … │
│ beta[4] 5.3456 0.2548 0.0046 3093.3317 2024.0479 1.0019 … │
│ alpha[4] 232.5042 2.6994 0.0411 4309.6893 2358.9022 1.0007 … │
│ beta[5] 6.5736 0.2414 0.0043 3203.9874 2291.2986 1.0003 … │
│ alpha[5] 231.5822 2.5989 0.0409 4043.5471 2283.2571 1.0000 … │
│ beta[6] 6.1711 0.2354 0.0038 3910.0379 2177.9210 1.0028 … │
│ alpha[6] 249.7372 2.6493 0.0486 2971.9283 2206.0587 1.0009 … │
│ beta[7] 5.9792 0.2431 0.0045 2895.8654 2062.3430 1.0002 … │
│ alpha[7] 228.7119 2.5941 0.0466 3105.9915 2190.7383 1.0009 … │
│ beta[8] 6.4134 0.2414 0.0043 3206.3216 2262.0176 1.0009 … │
│ alpha[8] 248.3562 2.5897 0.0451 3302.0394 2055.2312 1.0016 … │
│ beta[9] 7.0526 0.2495 0.0046 2925.3491 2048.0961 0.9997 … │
│ alpha[9] 283.3009 2.6640 0.0466 3343.5696 2027.5865 0.9999 … │
│ beta[10] 5.8414 0.2462 0.0046 2809.4637 2067.5844 0.9997 … │
│ alpha[10] 219.2244 2.6685 0.0484 3052.7004 1935.4315 0.9997 … │
│ beta[11] 6.8045 0.2460 0.0045 2965.3322 2411.8745 1.0007 … │
│ alpha[11] 258.2240 2.5662 0.0389 4363.5508 2183.5659 1.0000 … │
│ beta[12] 6.1205 0.2367 0.0042 3252.9509 2213.4161 0.9998 … │
│ alpha[12] 228.0912 2.6275 0.0433 3684.5965 1704.4568 0.9999 … │
│ beta[13] 6.1659 0.2323 0.0040 3396.6531 2376.6685 1.0013 … │
│ alpha[13] 242.4278 2.6568 0.0451 3473.4221 2453.5960 0.9997 … │
│ beta[14] 6.6867 0.2437 0.0042 3442.2296 1968.3355 0.9998 … │
│ alpha[14] 268.1953 2.7047 0.0436 3840.7295 2259.2321 1.0000 … │
│ beta[15] 5.4235 0.2494 0.0048 2738.2956 2051.5517 1.0000 … │
│ alpha[15] 242.7645 2.5769 0.0422 3744.1764 2167.2506 0.9997 … │
│ beta[16] 5.9244 0.2440 0.0040 3679.2134 2134.9486 0.9998 … │
│ alpha[16] 245.3284 2.6826 0.0439 3733.4199 2329.8764 1.0000 … │
│ beta[17] 6.2737 0.2385 0.0038 4014.3519 2607.6651 1.0001 … │
│ alpha[17] 232.1710 2.6324 0.0443 3527.1297 2275.0776 1.0000 … │
│ beta[18] 5.8412 0.2487 0.0045 2991.9525 2137.9996 1.0001 … │
│ alpha[18] 240.5028 2.6340 0.0432 3726.4546 2260.9624 1.0000 … │
│ beta[19] 6.4022 0.2370 0.0037 4165.7489 2416.6153 0.9997 … │
│ alpha[19] 253.8298 2.6512 0.0471 3148.2038 2307.2869 1.0001 … │
│ beta[20] 6.0572 0.2373 0.0035 4496.3798 2128.2280 0.9999 … │
│ alpha[20] 241.6608 2.6016 0.0451 3343.8570 2248.3963 1.0007 … │
│ beta[21] 6.4115 0.2338 0.0042 3162.9775 2396.7413 1.0002 … │
│ alpha[21] 248.5867 2.6598 0.0469 3234.1170 1973.6722 0.9999 … │
│ beta[22] 5.8632 0.2375 0.0039 3677.5644 2282.9499 1.0000 … │
│ alpha[22] 225.2980 2.6515 0.0422 3947.0257 2201.0198 1.0001 … │
│ beta[23] 5.7516 0.2410 0.0042 3345.6429 2103.4170 1.0007 … │
│ alpha[23] 228.4433 2.6342 0.0457 3342.1569 2031.8502 0.9998 … │
│ beta[24] 5.8917 0.2363 0.0041 3294.0714 1831.4341 1.0000 … │
│ alpha[24] 245.0350 2.5296 0.0403 3939.3047 2487.3701 1.0031 … │
│ beta[25] 6.9075 0.2514 0.0047 2886.3491 2469.1145 1.0010 … │
│ alpha[25] 234.5209 2.6028 0.0442 3464.1181 2505.6424 1.0013 … │
│ beta[26] 6.5425 0.2369 0.0044 2892.6079 1917.1288 1.0031 … │
│ alpha[26] 254.0598 2.7081 0.0495 2989.1272 2530.4940 1.0008 … │
│ beta[27] 5.8976 0.2450 0.0046 2844.4121 2061.2707 0.9997 … │
│ alpha[27] 254.3538 2.7215 0.0447 3702.3326 2125.5234 1.0002 … │
│ beta[28] 5.8447 0.2501 0.0041 3772.1464 1926.9101 0.9999 … │
│ alpha[28] 242.9818 2.6780 0.0443 3638.8879 1943.4626 0.9998 … │
│ beta[29] 5.6677 0.2470 0.0046 2874.7861 1823.2879 1.0002 … │
│ alpha[29] 217.8905 2.6208 0.0442 3525.0192 1954.6036 1.0002 … │
│ beta[30] 6.1313 0.2408 0.0036 4557.2619 2194.6886 1.0000 … │
│ alpha[30] 241.3804 2.5367 0.0413 3757.4368 2196.6919 0.9997 … │
│ alpha0 106.6125 3.6076 0.0663 2939.2770 2193.7745 0.9999 … │
│ sigma 6.0652 0.4568 0.0118 1546.6512 1804.5013 1.0002 … │
╰──────────────────────────────────────────────────────────────────────────────╯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 2000 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.