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.0763 1.4674 0.0286 2764.3872 2311.7847 0.9998 … │
│ beta_c 6.1849 0.1070 0.0018 3634.9853 2419.2716 1.0004 … │
│ alpha_tau 0.0050 0.0014 0.0000 4435.7530 2015.5669 1.0015 … │
│ alpha_c 242.7358 2.6774 0.0436 3800.0426 2214.4307 0.9999 … │
│ tau_c 0.0274 0.0042 0.0001 1805.0123 2340.5866 1.0001 … │
│ beta[1] 6.0627 0.2486 0.0036 4819.6371 2134.9486 0.9999 … │
│ alpha[1] 239.8982 2.8007 0.0394 5051.4672 2041.0419 0.9997 … │
│ beta[2] 7.0544 0.2520 0.0040 4037.0020 2329.9783 1.0000 … │
│ alpha[2] 247.7774 2.6981 0.0401 4542.9372 1997.1315 1.0005 … │
│ beta[3] 6.4849 0.2368 0.0034 4897.7329 2062.8128 1.0025 … │
│ alpha[3] 252.4837 2.5571 0.0385 4372.6522 2302.2178 0.9997 … │
│ beta[4] 5.3420 0.2532 0.0041 3862.8424 2090.9726 1.0001 … │
│ alpha[4] 232.6132 2.6912 0.0368 5311.8530 1971.2013 1.0019 … │
│ beta[5] 6.5732 0.2423 0.0037 4382.8942 2289.0888 1.0042 … │
│ alpha[5] 231.5901 2.6377 0.0369 5106.6277 2310.0675 1.0021 … │
│ beta[6] 6.1707 0.2388 0.0035 4717.4411 1879.8724 1.0004 … │
│ alpha[6] 249.7049 2.8493 0.0442 4149.0738 1965.8556 0.9997 … │
│ beta[7] 5.9767 0.2387 0.0033 5188.9819 1599.3278 1.0003 … │
│ alpha[7] 228.7023 2.6885 0.0430 3900.2319 1668.9370 1.0003 … │
│ beta[8] 6.4183 0.2448 0.0035 4826.7815 2015.1855 0.9997 … │
│ alpha[8] 248.3580 2.6886 0.0380 5070.9457 2145.4261 1.0009 … │
│ beta[9] 7.0531 0.2566 0.0041 3946.9531 2068.9243 0.9999 … │
│ alpha[9] 283.2932 2.7772 0.0397 4949.2927 1679.2674 0.9998 … │
│ beta[10] 5.8439 0.2486 0.0035 5004.7374 2059.9734 1.0003 … │
│ alpha[10] 219.3511 2.7293 0.0425 4138.6106 2146.5444 1.0001 … │
│ beta[11] 6.7957 0.2558 0.0038 4411.5178 2472.7518 1.0010 … │
│ alpha[11] 258.2061 2.6611 0.0400 4434.9896 2132.4447 1.0005 … │
│ beta[12] 6.1149 0.2437 0.0033 5278.6529 2065.6060 1.0026 … │
│ alpha[12] 228.1469 2.6724 0.0411 4235.1292 1842.6006 0.9997 … │
│ beta[13] 6.1641 0.2464 0.0035 4863.6129 2283.6616 1.0003 … │
│ alpha[13] 242.4670 2.7356 0.0400 4629.9326 2036.6151 1.0005 … │
│ beta[14] 6.6852 0.2529 0.0039 4233.7505 2360.1791 0.9999 … │
│ alpha[14] 268.3029 2.6738 0.0382 4920.4403 2219.5541 0.9999 … │
│ beta[15] 5.4142 0.2505 0.0037 4506.4125 2057.6051 1.0009 … │
│ alpha[15] 242.7567 2.6767 0.0381 4974.4580 1993.5883 1.0007 … │
│ beta[16] 5.9222 0.2441 0.0038 4185.3295 2237.8969 1.0001 … │
│ alpha[16] 245.3079 2.7111 0.0400 4603.7536 1781.7200 1.0000 … │
│ beta[17] 6.2716 0.2363 0.0033 5313.6255 2214.6970 0.9997 … │
│ alpha[17] 232.1187 2.6883 0.0392 4708.7337 2257.9634 1.0006 … │
│ beta[18] 5.8473 0.2519 0.0037 4712.2375 1837.0237 0.9997 … │
│ alpha[18] 240.4560 2.5488 0.0397 4110.7285 1837.9277 0.9998 … │
│ beta[19] 6.4050 0.2418 0.0034 4984.2494 2237.5389 0.9999 … │
│ alpha[19] 253.8487 2.7483 0.0398 4763.3881 1967.4210 1.0002 … │
│ beta[20] 6.0558 0.2361 0.0034 4819.2653 2091.2681 1.0017 … │
│ alpha[20] 241.6459 2.6407 0.0375 4906.7386 2173.4245 1.0041 … │
│ beta[21] 6.4042 0.2419 0.0035 4748.5811 2167.3161 0.9998 … │
│ alpha[21] 248.6335 2.7157 0.0387 4886.8832 1858.0954 1.0042 … │
│ beta[22] 5.8663 0.2511 0.0038 4474.4861 2179.6111 1.0020 … │
│ alpha[22] 225.1966 2.6509 0.0391 4587.1230 1903.4376 1.0005 … │
│ beta[23] 5.7473 0.2483 0.0040 3849.3431 2085.3684 1.0004 … │
│ alpha[23] 228.5482 2.7326 0.0380 5177.4237 2349.8163 1.0001 … │
│ beta[24] 5.8909 0.2442 0.0036 4469.3510 1972.6378 1.0000 … │
│ alpha[24] 245.0887 2.7120 0.0387 4703.6279 2081.8760 1.0011 … │
│ beta[25] 6.9127 0.2465 0.0040 3824.1887 2166.8360 1.0001 … │
│ alpha[25] 234.4694 2.6552 0.0355 5600.6588 2140.1585 1.0006 … │
│ beta[26] 6.5441 0.2489 0.0038 4377.9315 2045.1557 1.0002 … │
│ alpha[26] 254.0078 2.6180 0.0396 4376.8205 2043.1839 1.0000 … │
│ beta[27] 5.9026 0.2402 0.0033 5300.3202 2173.1442 1.0024 … │
│ alpha[27] 254.3232 2.6874 0.0391 4718.0001 2116.8228 0.9999 … │
│ beta[28] 5.8452 0.2463 0.0037 4487.6363 1913.4700 0.9999 … │
│ alpha[28] 242.9436 2.6546 0.0402 4373.9207 1916.1474 1.0004 … │
│ beta[29] 5.6736 0.2474 0.0035 4908.5770 1975.7489 1.0006 … │
│ alpha[29] 217.8982 2.5989 0.0402 4226.0067 1654.3743 1.0007 … │
│ beta[30] 6.1269 0.2467 0.0035 5100.0730 2297.9521 1.0000 … │
│ alpha[30] 241.4257 2.6816 0.0399 4550.0167 2168.6457 1.0007 … │
│ alpha0 106.6670 3.6044 0.0606 3548.1964 2375.1571 0.9999 … │
│ sigma 6.0933 0.4667 0.0110 1805.0123 2340.5866 1.0001 … │
╰──────────────────────────────────────────────────────────────────────────────╯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.