Introduction to Bayesian Data Analysis
Chenzi Xu
MPhil DPhil (Oxon)
University of York
2022/12/04 (updated: 2022-12-12)
Outline
- Why Bayesian data analysis?
- Fundamentals
- Random variable and distribution
- Conditional probability
- Marginal likelihood
- Bayes’ rule
- Bayesian analysis
- Computational Bayesian models in R
Why Bayesian?
1 2 3 4 5
Frequentist Framework
1 2 3 4 5
(magic!)
The Central Limit Theorem
The sampling distribution of means is normal, provided that:
- the sample size is large enough
- μ and σ are defined for the probability density or mass function that generated the data
Illustration from Wikipedia
Frequentist Framework
1 2 3 4 5
Hypothetical repeated sampling
The Exponential distribution has a parameter λ.
Mean: λ;
Variance: 1/λ2
Frequentist Framework
1 2 3 4 5
Confidence Interval (CI)
It represents the range of plausible value of the μ parameter.
If you take samples repeatedly and compute the CI each time, 95% of those CIs will contain the true population mean μ.
Frequentist Framework
1 2 3 4 5
p-value
The probability that the null hypothesis is true, or the probability that the alternative hypothesis is false.
The probability of obtaining the observed sample statistic, or some value more extreme than that, conditional on the assumption that the null hypothesis is true.
Frequentist Framework
1 2 3 4 5
Errors and Power
- Type I, II, M, S errors
- Power (1-Type II error) is a function of:
- the effect size
- the standard deviation
- the sample size.
Null hypothesis significance testing (NHST) is only meaningful when power is high.
Why Bayesian?
1 2 3 4 5
- Friendly to limited data
- Intuitive and solid model
- Straightforward interpretation
- Model flexibility
- Not limited by optimisation constraints
Fundamentals
1 2 3 4 5
Random variable and distribution
1 2 3 4 5
Data: the underlying random variable (Y) producing the data
Discrete random variable
Examples:
- Acceptability ratings on a Likert scale
- Binary grammaticality judgements
Probability distribution p(Y):
- Probability Mass Function (PMF)
- PMF maps each possible outcome of Y to a value between [0, 1].
Continuous random variable
Examples:
- Reading times or reaction times (in ms)
- EEG signals (in microvolts)
- f0, F1, F2…
Probability distribution p(Y):
- Probability Density Function (PDF)
- PDF maps a range of values of possible outcomes of Y to a value between [0, 1].
Discrete random variable: Binomial distribution
1 2 3 4 5
PMF for binomial distribution:
f(k∣n,θ)=(nk)θk(1−θ)n−kHere, n represents the total number of trials, k the number of successes, and θ the probability of success. The term (nk), the number of ways in which one can choose k successes out of n trials, expands to n!k!(n−k)!.
Continuous random variable: Normal distribution
1 2 3 4 5
PDF for normal distribution Xiid∼Normal(μ,σ):
f(x∣μ,σ)=1√2πσ2exp(−(x−μ)22σ2)Assumption: Each data point in the vector of data y is independent of the others.
The kernel density of the normal PDF: g(x|μ,σ)=exp(−(x−μ)22σ2), and k∫g(x)dx=1.
Probability in a continuous distribution is the area under the curve P(X<u)=∫u−∞f(x)dx, and will always be zero at any point value P(X=x)=0.
Continuous random variable: Normal distribution
1 2 3 4 5
Expectation and Variance
1 2 3 4 5
Expectation E(x): the weighted mean of the possible outcomes, weighted by the respective probabilities of each outcome.
Variance is defined as: Var(X)=E[X2]−E[Y]2
Discrete: Binomial
Expectation: E[Y]=∑yy⋅f(y)=nθ
Variance: Var(X)=nθ(1−θ)
- Estimated ˆθ=kn (Maximum likelihood estimate of the true parameter θ)
Continuous: Normal
Expectation: E[X]=∫xf(x)dx=μ
Variance: Var(X)=σ2
- Estimated ˆμ=ˉμ=∑yn
- Estimated ˆσ2=∑(y−ˉy)2n−1 (MLE estimate)
Likelihood Function
1 2 3 4 5
The likelihood function L(θ∣k,n) refers to the PMF p(k|n,θ), treated as a function of θ.
Suppose that we record n=10 trials, and observe k=7 successes (heads in coin tosses).
L(θ∣k=7,n=10)=(107)θ7(1−θ)10−7
Marginal likelihood
1 2 3 4 5
The concept of“Integrating out a parameter”
Marginal likelihood: the likelihood computed by “marginalizing” out the parameter θ. It is a kind of weighted sum of the likelihood, weighted by the possible values of the parameter.
Bivariate distributions
1 2 3 4 5
Discrete case
We have a joint PMF pX,Y(x,y) for each possible pair of values of X and Y.
The marginal distributions: pX(x)=∑y∈SYpX,Y(x,y)
pY(y)=∑x∈SXpX,Y(x,y)The conditional distributions: pX∣Y(x∣y)=pX,Y(x,y)pY(y)
pY∣X(y∣x)=pX,Y(x,y)pX(x)Bivariate distributions
1 2 3 4 5
Continuous case
PDF and CDF (ρ=0):
The joint distribution of X and Y: (XY)∼N2((μXμY),(σ2XρXYσXσYρXYσXσYσ2Y))
In the variance-covariance matrix, ρXY is the correlation between X and Y; their covariance is Cov(X,Y)=ρXYσXσY.
The area under the curve sums to 1. ∬SX,YfX,Y(x,y)dxdy=1
The marginal distributions:
fX(x)=∫SYfX,Y(x,y)dyfY(y)=∫SXfX,Y(x,y)dxConditional probability
1 2 3 4 5
The probability of A given B: P(A|B)=P(A,B)P(B)=P(B|A)P(A)P(B) where P(B)>0
Since the probability A and B happening is the same as the probability of B and A happening P(B,A)=P(A,B), we have:
P(B,A)=P(B|A)P(A)=P(A|B)P(B)=P(A,B).Rearranging terms:
P(B|A)=P(A|B)P(B)P(A)This is Bayes’ rule.
Bayes’ rule
1 2 3 4 5
Bayes’ rule
1 2 3 4 5
If y is a vector of data points, we have:
f(θ∣y)=f(y∣θ)f(θ)f(y)Here f(⋅) is a PDF (continuous case) or a PMF (discrete case).
We can rewrite it as follows:
Posterior=Likelihood⋅PriorMarginalLikelihoodPrior f(θ): the initial probability distribution of the parameter(s), before seeing the data.
Posterior f(θ∣y): the probability distribution of the parameters conditional on the data.
Likelihood f(y∣θ): the PMF (discrete case) or the PDF (continuous case) expressed as a function of parameters θ.
Marginal Likelihood f(y): It standardizes the posterior distribution to ensure that the area under the curve of the distribution sums to 1.
Bayes’ rule
1 2 3 4 5
If y is a vector of data points, we have:
f(θ∣y)=f(y∣θ)f(θ)f(y)Here f(⋅) is a PDF (continuous case) or a PMF (discrete case).
We can rewrite it as follows:
Posterior=Likelihood⋅PriorMarginalLikelihoodMarginal Likelihood f(y): It functions as the “normalisating constant”.
f(y)=∫f(y,θ)dθ=∫f(y∣θ)f(θ)dθWithout it, we have: f(θ∣y)∝f(y∣θ)f(θ)
Posterior∝Likelihood×PriorNHST vs Bayes
1 2 3 4 5
Frequentist
P(Data∣Theory)
- No prior knowledge
- Quantifies long-run probability of finding a false positive
- Hard cut-off decisions
Bayesian
P(Theory∣Data)
- Incorporates prior knowledge
- Quantifies uncertainty around possible parameter values
- Gradual assessment of evidence
Bayesian analysis
1 2 3 4 5
Deriving the Posterior: Analytical
1 2 3 4 5
Binomial likelihood and Beta prior
1. Choosing a likelihood
The responses follow a binomial distribution:
f(x∣n,θ)=(xn)θx(1−θ)n−xDeriving the Posterior: Analytical
1 2 3 4 5
Binomial likelihood and Beta prior
2. Choosing a prior distribution
We need a distribution that can represent our uncertainty about the probability θ of success.
The beta distribution is commonly used as prior for proportions.
f(θ)={1B(a,b)θa−1(1−θ)b−1if 0<θ<10otherwiseWhere B(a,b)=∫10θa−1(1−θ)b−1dθ, a normalizing constant that ensures that the area under the curve sums to 1.
The expectation and variance of Beta PDF:
E[θ]=aa+b Var(θ)=ab(a+b)2(a+b+1)
Deriving the Posterior: Analytical
1 2 3 4 5
Binomial likelihood and Beta prior
3. Conjugate analysis based on Bayes’ rule
To get the posterior distribution for θ, we have (ignore the normalising constant (xn)):
f(θ∣x)∝f(x∣n,θ)f(θ)We multiply the likelihood and the prior: θx(1−θ)n−xθa−1(1−θ)b−1=θa+x−1(1−θ)b+n−x−1
We computed the posterior up to proportionality.
Thus, given a Binomial(n,x∣θ) likelihood, and a Beta(a,b) prior on θ, the posterior will be Beta(a+x,b+n−x).
Deriving the Posterior: Analytical
1 2 3 4 5
The conjugate case of the beta-binomial:
- There are other similar conjugate cases such as the gamma-poisson case.
- The posterior mean is a compromise between the prior mean and the sample mean.
- Given sparse data, informative priors will play an important role.
- Expert knowledge or gut feelings
- Theoretical predictions
- Existing data
- Meta-analysis
In realistic data-analysis settings, the likelihood function will be very complex, and many parameters will be involved.
Computational Bayesian models in R
1 2 3 4 5
Deriving the Posterior: Sampling
1 2 3 4 5
Sampling algorithms
MCMC (Markov Chain Monte Carlo) sampling
- Gibbs sampling
- Random walk Metropolis
- Hamiltonian Monte Carlo (HMC, by Stan)
Bayes’ rule: p(Θ|Y)=p(Y|Θ)⋅p(Θ)p(Y)p(Θ|Y)=p(Y|Θ)⋅p(Θ)∫Θp(Y|Θ)⋅p(Θ)dΘ
MCMC Sampling
1 2 3 4 5
The inversion method
It works when we know the closed form of the PDF we want to simulate from and can derive the inverse of that function.
- Sample one number u from Unif(0,1). Let u=F(z)=∫zLf(x)dx.
- z=F−1(u) is a draw from f(x).
MCMC Sampling
1 2 3 4 5
HMC demonstration
- High-dimensional (hierarchical) models
- Only works with continuous parameters
- Based on Hamiltonian dynamics Energy(θ,k)=U(θ)+KE(k)
In trace plot, we see a characteristic “fat hairy caterpillar” shape, and we say that the sampler has converged.
The first few samples (warm-ups) were dropped, because the initial samples may not yet be sampling from the posterior.
We usually run four parallel chains with different initial values chosen randomly. They all end up sampling from the same distribution and overlap visually (mixing).
Bayesian Analysis using Stan: brms
1 2 3 4 5
Data:
Suppose we have data from a single subject repeatedly pressing the space bar as fast as possible, without paying attention to any stimuli. The data are response times in milliseconds in each trial.
RQ: How long it takes to press a key for this subject?
Bayesian Analysis using Stan: brms
1 2 3 4 5
Bayesian Analysis using Stan: brms
1 2 3 4 5
library(brms) fit_press <- brm(rt ~ 1, data = df_spacebar, family = gaussian(), prior = c( prior(uniform(0, 60000), class = Intercept, lb = 0, ub = 60000), prior(uniform(0, 2000), class = sigma, lb = 0, ub = 2000) ), chains = 4, iter = 2000, warmup = 1000 )library(brms) fit_press <- brm(rt ~ 1, data = df_spacebar, family = gaussian(), prior = c( prior(uniform(0, 60000), class = Intercept, lb = 0, ub = 60000), prior(uniform(0, 2000), class = sigma, lb = 0, ub = 2000) ), chains = 4, iter = 2000, warmup = 1000 )library(brms) fit_press <- brm(rt ~ 1, data = df_spacebar, family = gaussian(), prior = c( prior(uniform(0, 60000), class = Intercept, lb = 0, ub = 60000), prior(uniform(0, 2000), class = sigma, lb = 0, ub = 2000) ), chains = 4, iter = 2000, warmup = 1000 )library(brms) fit_press <- brm(rt ~ 1, data = df_spacebar, family = gaussian(), prior = c( prior(uniform(0, 60000), class = Intercept, lb = 0, ub = 60000), prior(uniform(0, 2000), class = sigma, lb = 0, ub = 2000) ), chains = 4, iter = 2000, warmup = 1000 )
SAMPLING FOR MODEL 'f7a7fc24ef4802fd1ddd47ff5ebae03c' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.00019 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.9 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.047038 seconds (Warm-up)
Chain 1: 0.032751 seconds (Sampling)
Chain 1: 0.079789 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'f7a7fc24ef4802fd1ddd47ff5ebae03c' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 1.8e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.18 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.043865 seconds (Warm-up)
Chain 2: 0.037392 seconds (Sampling)
Chain 2: 0.081257 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'f7a7fc24ef4802fd1ddd47ff5ebae03c' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 1.5e-05 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.040491 seconds (Warm-up)
Chain 3: 0.042745 seconds (Sampling)
Chain 3: 0.083236 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'f7a7fc24ef4802fd1ddd47ff5ebae03c' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 1.6e-05 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.16 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.041945 seconds (Warm-up)
Chain 4: 0.036886 seconds (Sampling)
Chain 4: 0.078831 seconds (Total)
Chain 4: Sampling Specifications:
chains:the number of independent runs for sampling (default: 4)
iter: the number of iterations that the sampler makes to sample from the posterior distribution of each parameter (default: 2000)
warmup: the number of iterations from the start of sampling that are eventually discarded (default: half of iter)
Bayesian Analysis using Stan: brms
1 2 3 4 5
Bayesian Analysis using Stan: brms
1 2 3 4 5
# A draws_df: 3 iterations, 1 chains, and 4 variables
b_Intercept sigma lprior lp__
1 168 24 -19 -1683
2 169 24 -19 -1683
3 168 26 -19 -1683
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}[1] 168.6334 2.5% 97.5%
166.0657 171.2347 [1] 24.99585 2.5% 97.5%
23.29355 26.88645 Bayesian Analysis using Stan: brms
1 2 3 4 5
Sensitivity Analysis
- Flat, uninformative priors: E.g., μ∼Uniform(−1020,1020)
- Regularizing priors: E.g., μ∼Normal+(0,1000)
- Principled priors: E.g., μ∼Normal+(250,100)
- Informative priors: E.g., μ∼Normal+(200,20)
Bayesian Analysis using Stan: brms
1 2 3 4 5
This sensitivity analysis showed that the posterior is not overly affected by the choice of prior.
Bayesian Analysis using Stan: brms
1 2 3 4 5
Bayesian Workflow:
- Define the model (likelihood and prior)
- Check if priors are reasonable
- Prior predictive distribution
- Sensitivity analysis
- Fit model(s)
- Check if posterior prediction matches data
- Posterior predictive distribution
- p(Dpred∣y)=∫Θp(Dpred∣Θ)p(Θ∣y)dΘ
References and Resources
This talk is dedicated to the SMLP 2022. All the materials are from https://vasishth.github.io/.
- Bob Carpenter’s Probability and Statistics: a simulation-based introduction
- Lambert (2018) and Lynch (2007)
- Michael Betancourt’s introduction
- https://youtu.be/jUSZboSq1zg