Day 4 Afternoon · Accessible lecture notes

Bayesian statistics I

A screen-reader-friendly transcript of the lecture slides: all slide text, tables, figure descriptions, and speaker notes.

About these notes

Each section is one slide. Figures carry text descriptions where available. For the concepts with worked examples, also read the plain-language summary and lab for this session.

Slide 1

Introduction to Bayesian Analysis

Figure 1.1. 3D surface plot of a bivariate normal probability density over axes x1 and x2, a bell-shaped peak rising to about 0.15, with a semi-transparent horizontal reference plane cutting through it.

Many slides adapted from Matt Falcy ODFW, and SESYNC course by Chris Che-Castaldo, Mary B. Collins, N. Thompson Hobbs (with permission)

Figure 1.2. Portrait of Thomas Bayes over a blue background with the handwritten formula P(A|B) = P(B|A) P(A) / P(B).

Slide 2

“Frequentist” or classical inference
Probability is a long-run relative frequency.
Parameters are fixed and unknown constants.
Statistical procedures are judged by how they perform over hypothetical replications of data that are never actually observed.

Bayesian inference
Probability is a degree of belief.
Once data are observed, they are treated as fixed (not random) values.
Parameters are unknown, therefore treated as random variables
Rules of probability used to revise “prior” beliefs in light of data.

Frequentist and Bayesian inference

Slide 3

Figure 3.1. Portrait of Thomas Bayes over a blue background with the handwritten formula P(A|B) = P(B|A) P(A) / P(B).

Remember Bayes Rule?

We derived this from simple rules of conditional probability

P(A,B) = P(A|B)P(B)

Slide 4

1% of women at age forty who participate in routine screening have breast cancer. 80% of women with breast cancer will get positive mammographies. 9.6% of women without breast cancer will also get positive mammographies. A woman in this age group had a positive mammography in a routine screening. What is the probability that she actually has breast cancer?

[c | + ] =

[ +|c ] * [c]

[ + ]

[+ |c ] * [c]

[ +| c]*[c] + [+ | not c]*[not c]

0.8 * 0.01

0.8*0.01 + 0.096*0.99

= 0.0776

95 out of 100 doctors put it at 0.70 to 0.80 (Eddy 1982, with replicated results)!!

Slide 5

We had complete information about medical testing probabilities

In almost all practical scenarios, we will just have data, y, and then want to make inference about parameters, θ.

Slide 6

[θ | y] =

[y | θ] [θ]

[y]

What is the marginal probability of the data?

L(θ|y) = [y| θ] ≠ [θ | y]

Bayes Rule

Prior

Posterior

Slide 7

Figure 7.1. Slide 'What is [y]?' showing a Venn-style diagram partitioning a region into events B1, B2, B3 overlapping event A, illustrating the law of total probability: [A] = sum over i of [A|Bi][Bi], and the continuous form as an integral.

Integrating B out of the joint distribution

Slide 8

[θ | y] =

[y | θ] [θ]

[y]

What is the marginal probability of the data?

[y] = [y| θ] [θ] dθ

Figure 8.1. Large integral sign, symbolizing integration.

It’s the integral (or sum) of the numerator!
This “scales” the numerator so that it integrates to 1.

Figure 8.2. Animated reaction clip of a man in a suit holding his head with both hands, an overwhelmed "mind blown" reaction.

L(θ|y) = [y| θ] ≠ [θ | y]

We call this a normalizing constant

Slide 9

Figure 9.1. Slide titled 'Joint Distribution of Continuous Random Variables' showing a 3D surface of a multivariate normal density, a bell-shaped peak over two axes X and Y.

Figure 9.2. Slide 'Marginal Distribution of Continuous Random Variables': a scatter of a bivariate cloud with its two marginal densities, p(X) and p(Y), plotted along the top and right margins as histograms/curves.

P(X, Y) or [X,Y]

P(X) or [X], and P(Y) or [Y]

Sum or integrate over all values of the random variable you DON’T care about
Or if you’ve sampled from the joint distribution, just take the values of the variable of interest

Slide 10

Figure 10.1. Slide '[y] is critical to Bayes' with Poisson likelihood and gamma-prior posterior formulas, plus two small curves: Likelihood versus theta (peak near 9) and Posterior versus theta (peak near 9).

E.g. Number of chocolate chips observed in cookies

For simple problem, can use numerical integration (what would we integrate)

Figure 10.2. R plot of a gamma(0.0001, 0.0001) prior density versus x from 0 to 20; the curve is extremely high near x=0 and drops steeply toward zero, showing a vague prior.

Slide 11

Figure 11.1. Slide '[y] is critical to Bayes' with Poisson likelihood and gamma-prior posterior formulas, plus small Likelihood and Posterior curves versus theta, both peaking near theta = 9.

E.g. Number of chocolate chips observed in cookies

For simple problem, can use numerical integration (what would we integrate)

Figure 11.2. R plot of a gamma(0.0001, 0.0001) prior density versus x from 0 to 20; density is very high near zero then drops sharply, illustrating a vague prior.

fun = function(x){
dpois(5,x)*dpois(10,x)*dpois(11,x)*dpois(12,x)*dpois(14,x)*dpois(9,x)*dpois(8,x)*dpois(6,x) #The likelihood
*dgamma(x,0.0001,0.0001)} #The prior
integrate(fun,0,20) #[y]; Integral of the numerator = 8.010413e-14
integrate(fun,0,20)$value/8.010413e-14
x=seq(0,20,0.1)
plot(x,fun(x),type="l")
plot(x,fun(x)/8.010413e-14,type="l")

Slide 12

[y| θ] “likelihood”

Figure 12.1. Likelihood curve versus p from 0 to 1 (blue), a bell shape peaking at about 0.23 near p = 0.42, for 5 heads out of 12 flips.

Data: We get 5 heads out of 12 flips. Binomial process. Make inference on p.
How would we have done this with maximum likelihood?

Figure 12.2. Vague prior: a flat probability-density line at 1.0 across p from 0 to 1, the Beta(1,1) uniform prior.

[θ]
“vague prior”

Beta_pdf(1,1)

Binomial Example

Slide 13

Figure 13.1. Curve of the product [y|theta][theta] versus p (orange), a bell shape peaking at about 0.23 near p = 0.42; the unnormalized prior times likelihood.

How would we scale this so that it integrates to 1 and becomes a “real distribution”?

[y|θ][θ] Still doesn’t integrate to 1

Slide 14

Figure 14.1. Curve of the product [y|theta][theta] versus p (orange), peaking at about 0.23 near p = 0.42, before normalizing to a distribution.

The integral is 0.077.
So lets divide every number on the y-axis by 0.077

Divide [y|θ][θ] by the are under it’s curve ([y])
AKA “integral”

Slide 15

Figure 15.1. Posterior curve versus p (orange), a bell shape peaking near 3.0 at about p = 0.42, now integrating to 1; note the rescaled y-axis.

> fun1<-function(x){dbinom(5,12,x)}
> integrate(fun1,0,1)
0.07692308 with absolute error < 8.5e-16

> fun2<-function(x){dbinom(5,12,x)*dbeta(x,1,1)/0.07692308}
> integrate(fun2,0,1)
1 with absolute error < 1.1e-14

Notice new scale on y-axis. NOW this integrates to 1 and is a proper distribution

Divide by marginal distribution of the data,
Which is just a “normalizing constant”

Slide 16

Figure 16.1. Vintage black-and-white photo of one young girl whispering into the ear of another, who looks surprised.

The Bayesian inferential paradigm requires subjective prior beliefs!

Slide 17

Figure 17.1. Posterior curve versus p (red), a bell peaking near 3.0 around p = 0.42 (vague Beta(1,1) prior case).

Figure 17.2. Two posteriors versus p: magenta curve peaking near 3.9 around p = 0.47 and a low cyan curve, comparing informative Beta(3,2) versus vague priors.

Figure 17.3. Curve versus p (blue) with legend Likelihood; a low bell peaking near 0.24 at p = 0.42, shown against the Beta(3,2) prior on the companion plot.

The integral of this thing is 0.092308

The integral of this thing is 1

We have 5 heads out of 12 flips – 5/12 = 0.417. Now let’s use an informative prior, Beta(3,2).

Why is posterior mean larger than 0.417?

Question: What is the mean of a Beta(3,2)? What are we saying with this?

E[X] = α / (α + β)

var[X] = αβ / ((α + β)² (α + β + 1))

Slide 18

Figure 18.1. Posterior curve versus p (magenta), a tall narrow bell peaking near 11 around p = 0.42, from 50 heads out of 120 flips; a low cyan curve lies near zero.

Figure 18.2. Product [y|theta][theta] curve versus p (cyan), a narrow bell peaking near 0.09 around p = 0.42, for 50 out of 120 flips (much lower than for 5 of 12).

Instead of 5 heads out of 12 flips, lets see what happens with 50 heads out of 120 flips (50/120 = 0.417)

The integral of this thing is 0.010036

The integral of this thing is 1

Figure 18.3. Plot versus p with legend: Likelihood (blue, nearly flat near zero) and Prior Beta(3,2) (red), a skewed bell peaking near 1.8 around p = 0.67.

The more data you have, the lower the likelihood.
(Think here. Why? Likelihood of datum vs. data)

50/120=0.41667

Look at uncertainty

Much closer to data mean than posterior mean now

Slide 19

Think about what we just did

We multiplied the prior and likelihood
Each value was divided by the integral of the prior*likelihood, to “normalize” or scale so that it integrates to 1 (thereby satisfying a definition of probability).
This gave us a posterior for p to make probabilistic inference

[θ|y] = [y|θ] [θ] / ∫θ [y|θ] [θ] dθ
(posterior = likelihood × prior ÷ marginal distribution of the data)

Slide 20

Figure 20.1. Cartoon of a man in top hat and riding boots beating a dead horse lying on the ground.

Slide 21

Figure 21.1. Slide 'The components of Bayes theorem' showing equation 12: posterior equals likelihood times prior over the marginal (integral of likelihood times prior over theta).

Slide 22

Figure 22.1. Slide 'So what?' with equation 15 for [theta|y] and a plot over theta from -6 to 6 of the prior PDF (green), observed likelihood (red), and their product/numerator (dashed blue); prior peaks left of the likelihood.

Slide 23

Figure 23.1. Slide 'So what?' with equation 17 (posterior = likelihood times prior over normalization) and a plot of prior PDF (green), likelihood (red), and posterior (blue) versus theta; the posterior sits between prior and likelihood.

If we transform this into many conditional distributions, then we can sample from it

Slide 24

Figure 24.1. Directed acyclic graph figure (Ruprecht 2021, Fig. 3): a hybrid spatial mark-resight plus spatial capture-recapture model, with parameter boxes feeding detection-data ovals (y.marking, y.genetics, y.resight, telemetry) via solid and dashed arrows.

The conditional distributions can be represented by a Directed Acyclic Graph (DAG) based on the dependencies of the process

Variables with arrow heads should be on LHS of conditionality. They depend on tails of arrows, which go on RHS of conditionality. Variables without arrow pointing to them need a prior.

Slide 25

Complex models are now in widespread use in fisheries and wildlife
E.g. Joel Ruprecht’s 2021 density estimation paper using spatial capture recapture and spatial mark resight. Can incorporate GPS data, marking process, integrate multiple datasets, etc.

Figure 25.1. Directed acyclic graph (DAG) of a hybrid gSMR + SCR model: parameter boxes and photos (trap, camera, GPS collar, tracks) connect by solid and dashed arrows to detection-data ovals for marking, genetics, resight, and telemetry.

Slide 26

But let’s start simpler!
Back to coin flip problem

Slide 27

Figure 27.1. Posterior curve versus p (magenta), a tall narrow bell peaking near 11 around p = 0.42, from 50 heads out of 120 flips.

Figure 27.2. Product [y|theta][theta] curve versus p (cyan), a narrow bell peaking near 0.09 around p = 0.42, for 50 out of 120 flips.

Instead of 5 heads out of 12 flips, lets see what happens with 50 heads out of 120 flips.

The integral of this thing is 0.010036

The integral of this thing is 1

Figure 27.3. Plot versus p with legend: Likelihood (blue, near zero) and Prior Beta(3,2) (red), a skewed bell peaking near 1.8 around p = 0.67.

The more data you have, the lower the likelihood.
(Think here. Why? Likelihood of datum vs. data)

50/120=0.41667

Look at uncertainty

Much closer to data mean than posterior mean now

Slide 28

Funny thing is that the beta distribution is the “conjugate prior” to the binomial distribution.

Figure 28.1. Laughing cartoon smiley-face emoji pointing outward.

Note: 3 and 2 come from the beta prior

“conjugate priors” have the special property that the posterior follows the same distribution as the prior. This gives us a ‘closed-form’ solution for the posterior

Slide 29

Mean of the NEW Beta distribution is

The new Beta distribution for the posterior of p

Now we can see precisely how much information the prior adds!
A Beta(1,1) prior is like observing 1 more success and 1 more failure because y increases by 1 and n increases by 2
A Beta(3,2) prior is like observing 3 more success and 2 more failures because y increases by 3 and n increases by 5
Large datasets typically swamp priors, but how would you minimize information in this prior?
That is called the Haldane prior – an “improper prior”. What is its functional form here?

Sampling from a known posterior is called “Gibbs sampling”

αnew = y + α, βnew = β + n − y

αnew / (αnew + βnew) = (y + α) / (α + β + n)

Slide 30

Beta-Binomial Walkthrough

Slide 31

Figure 31.1. Slide 'The components of Bayes theorem' showing equation 12: posterior equals likelihood times prior over the marginal (integral of likelihood times prior over theta).

To reiterate: Conjugacy

“The beta is the conjugate prior for the binomial distribution”
“The gamma is the conjugate prior for the Poisson distribution”
So what does that mean about our Poisson problem from earlier?
Let’s derive this together

Slide 32

Figure 32.1 — needs description This slide contains an image not yet described in text. Instructor: add a description.

Slide 33

Figure 33.1 — needs description This slide contains an image not yet described in text. Instructor: add a description.

Slide 34

Figure 34.1. Slide 'So what?' with equation 18 factoring a complex joint distribution into conditional distributions, plus a components-of-Bayes-theorem diagram.

Slide 35

Figure 35.1. Conjugacy relationships diagram: arrows from Geometric, Negative Binomial, Binomial, and Bernoulli into Beta (labeled p); and from Exponential, Lognormal, Poisson, Normal into Gamma (labeled lambda, tau, mu).

Conjugacy relationships

Note this is the “precision” or inverse of the variance

Slide 36

The elegance of conjugate priors is not available except for simple models.
Symbolic calculus usually not possible when extended into multiple dimensions (models with many parameters) and more complex probability distributions.
The numerical integrations we did earlier…
can’t be done for high-dimensional (many parameters) problems
http://math.stackexchange.com/questions/425782/why-does-monte-carlo-integration-work-better-than-naive-numerical-integration-in
https://en.wikipedia.org/wiki/Numerical_integration

Figure 36.1. Sad-face emoji with downturned mouth and teary blue eyes.

Unfortunately…

> fun1<-function(x){dbinom(5,12,x)}
> integrate(fun1,0,1)
0.07692308 with absolute error < 8.5e-16

(We did them to illustrate concepts. Not advice on solving real Bayesian problems)

Slide 37

Because we couldn’t solve high-dimensional integrals, Thomas Bayes’ (1702-1761) approach to inferential logic languished for two centuries!
In the mean time, Ronald Fisher and others developed an alternative approach that did not use priors and did not attempt to characterize parameters probabilistically.

Slide 38

Figure 38.1. Photograph of Monte Carlo, Monaco: a hillside city on a Mediterranean coastline.

Figure 38.2. Historical black-and-white portrait of mathematician Andrey Markov.

Figure 38.3. Photograph of a metal chain of linked steel loops, a visual pun on 'Markov chain'.

Markov chain Monte Carlo (MCMC)

Then comes WWII and some mathematical physicists at Los Alamos needed to evaluate high dimensional integrals.

Note: like Newton’s invention of calculus to study gravity/orbits

“Markov” means dependence only on the previous sample

“Monte Carlo” means estimation of properties of a distribution by random sampling from the distribution. Apparently named for the Grand Casino in Monte Carlo

Slide 39

We’re going to sample from the numerator

Figure 39.1. Product/numerator curve [y|theta][theta] versus p from 0 to 1 (cyan), a bell shape peaking near 0.30 at about p = 0.47; the unnormalized posterior we sample from.

Pick values of p. Keep them in the chain at a rate proportional to the relative value on the y-axis.

The Metropolis Algorithm

Figure 39.2. Black-and-white portrait of physicist Nicholas Metropolis (1915-1999).

Nicholas Metropolis
(1915-1999)

We say: ”The posterior is proportional to the likelihood times the prior”

Slide 40

For example

The symmetry requirement can be problematic if θ is constrained (e.g. standard deviations and proportions must be positive).

If this is too big, then you propose lots of values that don’t get accepted, which is inefficient. If too small there is too much autocorrelation and it takes forever to sample the full posterior

The Metropolis Algorithm

Figure 40.1. Black-and-white portrait of physicist Nicholas Metropolis (1915-1999).

Nicholas Metropolis
(1915-1999)

A small modification:
The “Metropolis-Hastings” algorithm allows for asymmetric proposals

Slide 41

Figure 41.1. Scatter plot of MCMC draws, sig.K versus mu.K for 100 iterations; early points trace a path from the corner then cluster around mu near 9.7 and sigma near 0.8.

Time series of an MCMC

Slide 42

Figure 42.1. Scatter plot of MCMC draws, sig.K versus mu.K for 1000 iterations; points form a dense cloud centered near mu = 9.8, sigma = 0.9 after the initial burn-in path.

Time series of an MCMC

Slide 43

Figure 43.1. Scatter plot of three MCMC chains (blue, green, red), sig.K versus mu.K over 1000 iterations; all three overlap in one dense cloud near mu = 9.8, sigma = 0.9, showing convergence.

Time series of 3 MCMCs

Slide 44

count

Figure 44.1. 3D histogram of the joint posterior over parameters Mu and Sigma; counts (color-scaled to about 4000) form a tall spike over the high-density region, near-zero elsewhere.

count

Here is the posterior of the joint distribution (both parameters) from which we also have the marginal distribution for each

Slide 45

Figure 45.1. Two-panel animation still: left, an empty 'Markov chains' plane of sigma versus mu labeled Iteration 1; right, a flat 'Posterior density' surface over mu and sigma before sampling begins.

https://mbjoseph.github.io/posts/2018-12-25-animating-the-metropolis-algorithm/

Slide 46

Figure 46.1. Slide 'So what?' with six small DAG examples (I-VI) of two- and three-node directed graphs, each with the corresponding joint-probability factorization using the chain rule.

Slide 47

https://chi-feng.github.io/mcmc-demo/

There are now many algorithms to solve this problem

Click on an algorithm below to view interactive demo:
Random Walk Metropolis Hastings
Adaptive Metropolis Hastings [1]
Hamiltonian Monte Carlo [2]
No-U-Turn Sampler [2]
Metropolis-adjusted Langevin Algorithm (MALA) [3]
Hessian-Hamiltonian Monte Carlo (H2MC) [4]
Gibbs Sampling
Stein Variational Gradient Descent (SVGD) [5]
Nested Sampling with RadFriends (RadFriends-NS) [6]

Slide 48

BUGS, JAGS, Nimble, STAN, and INLA oh my!

In 2008, my first Bayesian professor, Herbie Lee, thought WinBUGS shouldn’t be taught because it would make it too easy so that incompetent people would implement poorly specified models
Now BUGS (Bayesian Analysis Using Gibbs Sampling) is very common
OpenBUGS and WinBUGS needs Windows. JAGS was designed to work on multiple platforms.
Unlike OpenBUGS and WinBUGS, which have a graphical user interface, JAGS works on the command line and interfaces nicely with R. It can be run on UNIX clusters like OSU’s CQLS.
Nimble is similar to JAGS, but is more flexible and faster for many problems in Ecology!
STAN uses the very fast Hamiltonian MCMC algorithm, but doesn’t use BUGS language
INLA is extremely fast because it approximates the posterior instead of sampling. Used for a restricted subset of problems.

Slide 49

MSY = rK/4

Greatest population growth rate

Greatest per-capita growth rate

“Negative density dependence”

Estimating parameters in logistic population growth

dN/dt = rN(1 − N/K)

(1/N) dN/dt = r(1 − N/K) = r − rN/K

(1/N) dN/dt

Slide 50

JAGS PRIMER
estimating parameters in logistic population growth
From Hobbs and Che-Castaldo’s SESYNC course

Figure 50.1. Line graph of per-capita growth rate dN/(N dt) versus population size N: a straight red line starting at r on the y-axis and declining with slope -r/K to zero at N = K.

Observations of per-capita growth rate paired with population size
Ah ha, our scientific model is a line! We know what to do with lines!
Let’s try it with our standard tools.

(1/N)(dN/dt) = r − (r/K) N

Slide 51

reg=lm(GrowthRate~PopulationSize,data=Logistic)

Good Ole Regression Approach

Figure 51.1. R lm() summary for GrowthRate ~ PopulationSize: intercept 0.2018, slope -1.647e-04 (both p < 2e-16), residual SE 0.02789, R-squared 0.744.

Figure 51.2. R console output: coef(reg) gives intercept 0.2018 and slope -0.0001647; derived rmax = 0.2018 and K = 1224.99.

But what can we say about K, K/2 or MSY?

(1/N)(dN/dt) = r − (r/K) N

Slide 52

JAGS PRIMER
estimating parameters in logistic population growth

Joint posterior distribution of parameters is proportional to likelihood*prior
We want positive continuous values for r and K. Bring on the gamma!
The gamma is the conjugate prior for the precision term in Normal distribution
But why gamma(0.001, 0.001)

μi = r − r xi / K
[r, K, τ | y] ∝ ∏i=1n [yi | μi, τ] [r] [K] [τ]
[r, K, τ | y] ∝ ∏i=1n normal(yi | μi, τ) · gamma(K | .001, .001) · gamma(τ | .001, .001) · gamma(r | .001, .001)

Slide 53

Figure 53.1. Bayesian network / DAG for the logistic model: x_i points (dashed) to y_i, with parameters r, K and tau pointing into y_i; the joint factors as [y|r,K,tau][r][K][tau], with rules for reading a DAG.

“Bayesian Network” or
”Directed Acyclic Graph”

Slide 54

JAGS PRIMER
estimating parameters in logistic population growth

Figure 54.1. Equations: mu_i = r - r x_i / K and tau = 1/sigma^2, with the joint posterior proportional to product of normal likelihoods times uniform priors K~U(0,4000), sigma~U(0,5), r~U(0,2).

JAGS parses this code to set up MCMC sampling, including proposal distribution, tuning parameters, and returns MCMC chain for each parameter so we can approximate posterior distributions and find associated statistics (e.g., means, medians, standard deviations, quantiles).

Figure 54.2. Numbered list of the three parts of a JAGS model: 1. code for priors, 2. code for the deterministic model, 3. code for the likelihood(s).

In general, we always need

#Use “for loop” to loop over all the data

Normal distribution in JAGS is parameterized with precision, not variance or standard deviation

Slide 55

Some other priors

sink("LogisticJAGS.R")
cat("
model {
# priors
K ~ dunif(0, 4000)
r ~ dunif (0, 2)
sigma ~ dunif(0, 50)
tau <- 1/sigma^2
# likelihood
for(i in 1:n) {
mu[i] <- r - r/K * x[i]
y[i] ~ dnorm(mu[i], tau)
}
} ",fill = TRUE)
sink()

Use “sink” to keep all your code in one R script, but can also make JAGS code separately

Normal distribution in JAGS is parameterized with precision, not variance or standard deviation

μi = r − r xi / K
τ = 1 / σ²
[r, K, σ | y] ∝ ∏i=1n normal(yi | μi, τ) · uniform(K | 0, 4000) · uniform(σ | 0, 5) · uniform(r | 0, 2)

Slide 56

JAGS PRIMER
estimating parameters in logistic population growth

Unlike with least squares or maximum likelihood, with Bayesian statistics it is easy to derive chains for other quantities of interest and their posterior distributions.
These are called “derived quantities”
For example, K/2, rK/4, N as a function of time, dN/dt as a function of N with uncertainty, or incorporate errors in the observations of population growth rate or population size.

This is a BFD!

Slide 57

JAGS PRIMER

Figure 57.1. R code screenshot using sink() to write a JAGS logistic model file: uniform priors K~dunif(0,4000), r~dunif(0,2), sigma~dunif(0,50), tau=1/sigma^2, and a likelihood loop with y[i]~dnorm(mu[i],tau).

For small models, use the sink() function to write your JAGS model to file for input

Figure 57.2. Text: 'Notice that the for loop replaces the product from i=1 to n in the likelihood.'

JAGS is not a procedural language like R. It is not a set of instructions to be executed sequentially. It is a way to translate a mathematical statement into syntax that determine the MCMC sampler that JAGS uses behind the scenes

Remember another sink() to write the file

fill = TRUE to print on separate lines

Slide 58

JAGS PRIMER – slightly different approaches for rjags, jagsUI, and R2jags packages

Figure 58.1. R code screenshot: sorting the Logistic data, three lists of initial values for K, r, sigma, the data list, MCMC settings (n.adapt, n.update, n.iter = 1000/10000/10000), and the jags.model() call.

Figure 58.2. R code screenshot: update(jm, n.iter = n.update); zm = coda.samples(jm, variable.names = c('K','r','sigma','tau'), n.iter, n.thin = 1); MCMCsummary(zm).

Use jags.model() to run JAGS. Specify the data, initial values, and MCMC parameters

Use update() to add iterations to an existing model. Use coda.samples() to access the MCMC chains

JAGS uses this to optimize the MCMC parameters

The burn in

Posterior draws we keep

The burn in

The samples

Slide 59

JAGS PRIMER

Figure 59.1. R code screenshot: update(jm, n.iter = n.update); zm = coda.samples(jm, variable.names = c('K','r','sigma','tau'), n.iter, n.thin = 1); MCMCsummary(zm).

Use update() to add iterations to an existing model. Use coda.samples() to access the MCMC chains.
Use MCMCsummary() to summarize the parameters and derived quantities (if they exist)

Figure 59.2. MCMCsummary output table: posterior mean, sd, and 2.5/50/97.5% quantiles for K (~1238), r (~0.201), sigma (~0.0287), tau (~1258); all Rhat = 1, indicating convergence.

Note that Rhat near 1 means that there is good convergence. This is for model diagnostics.

Slide 60

Nimble PRIMER

Figure 60.1. Nimble model code screenshot: nimbleCode with uniform priors K~dunif(0,4000), r~dunif(0,2), sigma~dunif(0,2), tau=1/sigma^2, and a likelihood loop y[i]~dnorm(mu[i],tau).

BUGS Code is written in a nimbleCode() function.
Models compiled into C++ and are more customizable
Can be much faster than JAGS
Can produce MCMC output in one line with nimbleMCMC()
No separate adaption phase to tune MCMC
Can also compile model first and then run nimbleMCMC

Figure 60.2. R code: samples = nimbleMCMC(code=logistic, constants=data, inits=inits, monitors=c('K','r','sigma','tau','mu'), nchains=3, nburnin=1000).

Slide 61

Logistic Growth in Nimble. Visualize and interpret the output