Bayesian sampling method

1. Monte Carlo methods.
   • Direct sampling.
   • Rejection sampling.
   • Importance sampling.
   • Sampling-importance resampling.
2. Markov chain Monte Carlo (MCMC).
   • Metropolis-Hasting.
   • Gibbs sampling.

Limitation of Monte Carlo method

• Direct sampling
  • Often hard to compute the normalization factor \(Z\).
  • Hard to get rare events, especially in higher dimensional spaces.
• Rejection sampling, importance sampling.
  • Do not work well if the proposed distribution \(q^*(x)\) is very different from \(p(x)\).
  • Constructing a \(q^*(x)\) similar to \(p(x)\) can be difficult.
    • Making a good proposal usually requires knowledge of the analytic form of \(p(x)\) - but if we had that, we wouldn’t even need to sample!
• Solution: instead of a fixed proposed distribution \(q^*(x)\), we can use an adaptive proposal.

Motivation of Metropolis-Hastings

• Drawbacks of rejection sampling and SIR are that it is difficult to propose an efficient distribution \(q^*(x)\)
• For rejection sampling, it is also difficult to find \(M\).
• A smart idea is let the proposed distribution depends on the last accepted value.
• Instead of fixed \(q(x')\), we use \(q(x'|x)\) where \(x'\) is the new state being sampled and \(x\) is the previous sample.
• As \(x\) changes, \(q(x'|x)\) can also change (as a function of \(x'\).)

Comparison between importance sampling and Metropolis-Hastings

Monte Carlo vs Markov Chain Monte Carlo

• Monte Carlo methods: simulation/sampling methods.
  • Simulation
  • Rejection sampling
  • SIR
• Markov Chain Monte Carlo: A type of Monte Carlo method – next step sample depend on previous sample.
  • Metropolis-Hastings
  • Gibbs Sampling

Introduce MH algorithm

• 

  Comparison between importance sampling and Metropolis-Hastings

• Draw a sample \(x'\) from \(q(x'|x)\), where \(x\) is the previous sample.
• The new sample \(x'\) is accepted or rejected with some probability \(A(x'|x)\). \[A(x'|x) = \min \bigg(1, \frac{p(x')q(x|x')}{p(x)q(x'|x)} \bigg)\]

• \(\frac{p(x')}{q(x'|x)}/\frac{p(x)}{q(x|x')}\) is a ratio of importance sampling weights.

MH algorithm

• Initialize our parameters \(x^0\).
• Given an accepted value \(x^{t-1}\)
  1. Draw \(x \sim q(\cdot|x^{t-1})\).
  2. Accept \(x\) with probability \(A(x|x^{t-1}) = \min \bigg(1, \frac{p(x)q(x^{t-1}|x)}{p(x^{t-1})q(x|x^{t-1})} \bigg)\), and if accepted, set \(x^t = x\).
     • Draw \(u \sim UNIF(0,1)\).
     • If \(u < A(x|x^{t-1})\), accept.
  3. If we didn’t accept \(x\), set \(x^t = x^{t - 1}\)
• We repeat \(T\) times (\(t = 1, \ldots, T\)).

Implementation of MH algorithm

• target distribution: \(p(x) = \exp [ 0.4(x-0.4)^2 - 0.08x^4 ]\)
• sampling distribution: \(q(x^{(t)}) = dnorm(x^{(t)},x^{(t-1)},1)\)

p <- function(x, a=.4, b=.08){exp(a*(x-a)^2 - b*x^4)}
x <- seq(-4, 4, 0.1)
plot(x,p(x),type="l")

N <- 10000
x.acc5 <- rep.int(NA, N)
u <- runif(N)
acc.count <- 0
std <- 1 ## Spread of proposal distribution
xp <- 0; ## Starting value
for (ii in 1:N){
  set.seed(ii)
  xc <- rnorm(1, mean=xp, sd=std) ## xp previous sample; xc: current sample
  alpha <- min(1, (
      p(xc) /dnorm(xc, mean=xp,sd=std) /
      (p(xp) / dnorm(xp, mean=xc,sd=std))
                  )
               )
  x.acc5[ii] <- xp <- ifelse(u[ii] < alpha, xc, xp)
  ## find number of acccepted proposals:
  acc.count <- acc.count + (u[ii] < alpha)
}
## Fraction of accepted *new* proposals
acc.count/N

## [1] 0.7468

Check samples from MH

plot(x.acc5,type="l")

• Good convergence.
• Samples are focused at the underlying distribution (stationary distribution)

Vary some parameters in MH

• Burnin period: disgard intial samples since they may not be in the stationary distribution
  • initial x = 8?
• Play with variance? (doesn’t converge to the stationary distribution with limited sample size)
  • sd = 0.1
  • sd = 0.01

Trade off between acceptance rate and convergence.

• How many samples (empirically):
  • total 10,000 samples.
  • 500 burnin samples.

Remove initial values for burnin period

N <- 1000
x.acc5 <- rep.int(NA, N)
u <- runif(N)
acc.count <- 0
std <- 1 ## Spread of proposal distribution
xp <- 8; ## Starting value
for (ii in 1:N){
  xc <- rnorm(1, mean=xp, sd=std) ## xp previous sample; xc: current sample
  alpha <- min(1, (
      p(xc) /dnorm(xc, mean=xp,sd=std) /
      (p(xp) / dnorm(xp, mean=xc,sd=std))
                  )
               )
  x.acc5[ii] <- xp <- ifelse(u[ii] < alpha, xc, xp)
  ## find number of acccepted proposals:
  acc.count <- acc.count + (u[ii] < alpha)
}
## Fraction of accepted *new* proposals
acc.count/N

## [1] 0.73

plot(x.acc5,type="l")

Doesn’t converge

N <- 1000
x.acc5 <- rep.int(NA, N)
u <- runif(N)
acc.count <- 0
std <- 0.1 ## Spread of proposal distribution
xc <- 0; ## Starting value
for (ii in 1:N){
  xp <- rnorm(1, mean=xc, sd=std) ## proposal
  alpha <- min(1, (p(xp)/p(xc)) *
                 (dnorm(xc, mean=xp,sd=std)/dnorm(xp, mean=xc,sd=std)))
  x.acc5[ii] <- xc <- ifelse(u[ii] < alpha, xp, xc)
  ## find number of acccepted proposals:
  acc.count <- acc.count + (u[ii] < alpha)
}
## Fraction of accepted *new* proposals
acc.count/N

## [1] 0.971

plot(x.acc5,type="l")

Why MH converge?

Concepts:

• Assume \(x\) follows \(p(x)\)
• For any MCMC algorithm, the transition kernel is \(T(x'|x)\):
  • probablity of sampling \(x'\) when previous step is \(x\).
• Detailed balance condition:
  \[ p(x') = \int_{x} p(x)T(x'|x) dx\]
• If the transition kernel is \(T(x'|x)\) satisfies the detailed balance condition, the sampling distribution will be the underlying distribution.

Why MH converge?

• In MH we draw sample \(x'\) according to \(q(x'|x)\), then we accept/reject according to \(A(x'|x)\).
• The transition kernel is \(T(x'|x) = q(x'|x) A(x'|x)\).
• \(A(x'|x) = \min \bigg(1, \frac{p(x')q(x|x')}{p(x)q(x'|x)} \bigg)\)
• Proof:
  • If \(A(x'|x) \le 1\), then \(\frac{p(x')q(x|x')}{p(x)q(x'|x)} \le 1\), \(\frac{p(x)q(x'|x)}{p(x')q(x|x')} \ge 1\), \(A(x|x') = 1\).

\[ A(x'|x) = \frac{p(x')q(x|x')}{p(x)q(x'|x)}\] \[ p(x)q(x'|x) A(x'|x) = p(x')q(x|x') \] \[ p(x)q(x'|x) A(x'|x) = p(x')q(x|x') A(x|x')\] \[ p(x)T(x'|x) = p(x')T(x|x') \] The last line is called detailed balance condition.

\[ \int_{x} p(x)T(x'|x) dx = \int_{x} p(x')T(x|x') dx\] \[ p(x') = \int_{x} p(x)T(x'|x) dx\]

Since p(x) is the true distribution. MH algorithm will eventually converges to the true distribution.

Special cases for MH

• Metropolis algorithm:
  • The proposed distribution is symmetrical, e.g. \(q(x'|x) = q(x|x')\) for all pairs (x,x’). In this case the acceptance probability is \(A(x'|x) = \min (1, \frac{p(x')}{p(x)})\)
• Random-walk Metropolis:
  • A popular choice for proposal in a Metropolis algorithm is \(q(x'|x) = g(x-x')\) where g is symmetric.
• Independence sampler:
  • The proposed distribution \(q(x'|x) = q(x')\) doesn’t depend on \(x\). The acceptance probability becomes \(A(x'|x) = \min(1, \frac{p(x')}{p(x)} \frac{q(x)}{q(x')})\). This works well when \(q\) is a good approximation to \(p\).
• Gibbs sampling:
  • \(q(x'|x)\) is the conditional probability given all other variables. \(A(x'|x) = 1\).

Markov Chain Monte Carlo (MCMC)

• We always need to pay attention to:
  • Burn-in period
  • Monitor convergence
• Question: Since we know samples are correlated from Metropolis-Hastings algorithm, can we still estimate an estimator \(\theta\) by sample average?
  • Yes, if the underlying density function is also about \(\theta\).
  • Weak law of large number (WLLN) is still valid if samples are correlated, the only requirement is the samples are identical distributed.
• Another popular algorithm: Gibbs sampling.

Gibbs Sampling algorithm (Wikipedia)

• Gibbs sampling is named after the physicist Josiah Willard Gibbs, in reference to an analogy between the sampling algorithm and statistical physics.
• The algorithm was described by brothers Stuart and Donald Geman in 1984, some eight decades after the death of Gibbs.
• Josiah Willard Gibbs (February 11, 1839 - April 28, 1903) was an American scientist who made important theoretical contributions to physics, chemistry, and mathematics.

Gibbs Sampling motivating example

• Prior
  • \(P(I) = 0.5\), \(P(-I) = 0.5\)
• Likelihood (Generative process) for \(G\):
  • \(P(G|I) = 0.8\), \(P(-G|I) = 0.2\)
  • \(P(G|-I) = 0.5\), \(P(-G|-I) = 0.5\)
• Likelihood (Generative process) for \(S\):
  • \(P(S|I) = 0.7\), \(P(-S|I) = 0.3\)
  • \(P(S|-I) = 0.5\), \(P(-S|-I) = 0.5\)

Posterior distribution (via Bayes rule)

\[\begin{align*} P(I | G, S) &= \frac{P(I ,G, S)}{P(G, S)} \\ &= \frac{P(G, S|I)P(I)}{P(G, S|I)P(I) + P(G, S|-I)P(-I)} \\ &= \frac{P(G|I)P(S|I)P(I)}{P(G|I)P(S|I)P(I) + P(G|-I)P(S|-I)P(-I)} \\ &= \frac{0.8\times 0.7 \times 0.5}{0.8\times 0.7 \times 0.5 + 0.5\times 0.5 \times 0.5} \\ &= 0.69 \end{align*}\]

Similarly

• \(P(I | G, -S) = 0.49\)
• \(P(I | -G, S) = 0.36\)
• \(P(I | -G, -S) = 0.19\)

Gibbs Sampling algorithm

• Suppose the graphical model contains variable \(\theta_1, \ldots, \theta_p\).
• Initialize starting values for \(\theta_1, \ldots, \theta_p\).
• Do until convergence:
  • Pick an ordering of the \(p\) variables (can be fixed or random).
  • For each variable \(\theta_i\) in order:
    • Sample \(\theta\) from \(P(\theta_i| \theta_1, \ldots,\theta_{i-1}, \theta_{i+1},\ldots \theta_p,X)\), the conditional distribution of \(\theta_i\) given the current values of all other variables.
    • Update \(\theta_i \leftarrow \theta\)

Gibbs Sampling motivating example

1. Initialize I,G,S.
2. Pick an updating order. (e.g. I,G,S)
3. Update each individual variable given all other variables.

Implementation in R

I <- 1; G <- 1; S <- 1
pG <- c(0.5, 0.8)
pS <- c(0.5, 0.7)
pI <- c(0.19, 0.36, 0.49,0.69) ## see below for justification
for(g in c(0,1)){
  for(s in c(0,1)){
    cat("G:" ,g ,", S:", s, ", p(I):", pI[2*g + s + 1], "\n")
  }
}

## G: 0 , S: 0 , p(I): 0.19 
## G: 0 , S: 1 , p(I): 0.36 
## G: 1 , S: 0 , p(I): 0.49 
## G: 1 , S: 1 , p(I): 0.69

Frequentist and Bayesian philosophy

• Frequentist:
  • parameters \(\theta\) are fixed
  • data \(X\) are random variables
  • Goal: create procedures with frequency guarantees
• Bayesian:
  • parameters \(\theta\) are random variables
  • data \(X\) are random variables
  • analyze beliefs (distribution of \(\theta\) and \(X\))

An example of linear regression

• Underlying truth: \(n=100\), \(\alpha=0\), \(\beta=2\), \(\sigma^2=0.5\).
• Simulate iid samples: \(x_i \sim N(0,1)\), \(y_i \sim N(\alpha + \beta x_i,\sigma^2)\). (\(i = 1, \ldots, 100\))
• Purpose: want to infer \(\alpha\), \(\beta\), \(\sigma^2\).
• Set up priors: \(\alpha \sim N(\alpha_0, \tau_a^2)\), \(\beta \sim N(\beta_0, \tau_b^2)\), \(\sigma^2 \sim IG(\nu, \mu)\). \(\alpha_0 = 0\), \(\beta_0 = 0\), \(\tau_a^2 = 10\), \(\tau_b^2 = 10\), \(\nu=3\), \(\mu=3\).

Graphical model representation of the data generative process. Shaded nodes rerepsent abserved variables, dashed nodes represent priors.

• Posterior can be derived using Bayes rule:

  • \(var_\alpha \doteq 1/(1/\tau_a^2 + n/\sigma^2)\), \(var_\beta \doteq 1/(1/\tau_b^2 + \sum_{i=1}^n x_i^2/\sigma^2)\)

  • \(\alpha | x_1^n, y_1^n, \beta, \sigma^2 \sim N(var_\alpha(\sum_{i=1}^n (y_i - \beta x_i)/\sigma^2 + \alpha_0 / \tau_a^2), var_\alpha)\)

  • \(\beta | x_1^n, y_1^n, \alpha, \sigma^2 \sim N(var_\beta(\sum_{i=1}^n \big( (y_i - \alpha) x_i \big) /\sigma^2 + \beta_0 / \tau_b^2), var_\beta)\)

  • \(\sigma^2 | x_1^n, y_1^n, \alpha, \beta \sim IG(\nu + n/2, \mu + \sum_{i=1}^n (y_i - \alpha - \beta x_i)^2/2)\)

Demonstrate the R code

• MCMC simulation 1,000 times.
• Visualization of MCMC chain
• Remove burning 100.
• Distribution of parameter.
• Auto correlation. (ACF)
• Effective samples.
• Posterior mean.
• chain mixing.

Demonstrate the R code

set.seed(32611)
n = 100; alpha = 0; beta=2; sig2=0.5;true=c(alpha,beta,sig2)
x=rnorm(n)
y=rnorm(n,alpha+beta*x,sqrt(sig2))
# Prior hyperparameters
alpha0=0;tau2a=10;beta0=0;tau2b=10;nu0=3;mu0=3
# Setting up starting values
alpha=0;beta=0;sig2=1
# Gibbs sampler
M = 1000
draws = matrix(0,M,3)
draws[1,] <- c(alpha,beta,sig2)
for(i in 2:M){
  var_alpha = 1/(1/tau2a + n/sig2)
  mean_alpha = var_alpha*(sum(y-beta*x)/sig2 + alpha0/tau2a)
  alpha = rnorm(1,mean_alpha,sqrt(var_alpha))
  var_beta = 1/(1/tau2b + sum(x^2)/sig2)
  mean_beta = var_beta*(sum((y-alpha)*x)/sig2+beta0/tau2b)
  beta = rnorm(1,mean_beta,sqrt(var_beta))
  sig2 = 1/rgamma(1,(nu0+n/2),(mu0+sum((y-alpha-beta*x)^2)/2))
  draws[i,] = c(alpha,beta,sig2)
}

Check Gibbs sampling result

# Markov chain + marginal posterior
names = c('alpha','beta','sig2')
colnames(draws) <- names
ind = 101:M
par(mfrow=c(3,3))
for(i in 1:3){
  ts.plot(draws[,i],xlab='iterations',ylab="",main=names[i])
  abline(v=ind[1],col=4)
  abline(h=true[i],col=2,lwd=2)
  acf(draws[ind,i],main="")
  hist(draws[ind,i],prob=T,main="",xlab="")
  abline(v=true[i],col=2,lwd=2)
}

Posterior mean and median

• Posterior mean

colMeans(draws[101:M,])

##       alpha        beta        sig2 
## -0.06394806  1.93993716  0.54674891

• Posterior median

apply(draws[101:M,],2,median)

##       alpha        beta        sig2 
## -0.06487311  1.93895984  0.53787251

Conjugate prior in Gibbs

• The prior in the previous examples are all conjugate priors
• Conjugate priors will guarantee the posterior probablity (conditional distribution) has the same closed form formulation of the prior probablity
  • Try to use conjugate priors as much as possible.
• What if conjugate priors are not available, and the posterior probablity has no easy formulation, hard to sample?
  • Use MH to sample (HW).
  • Use MH within Gibbs interations for certain updates.

MCMC tradeoff for number of iterations

• Large number of iterations will tend to recover the underlying distribution in higher resolution.
• Large number of iterations will also add computational burden.

Burnin period 1

• If the initial value is within the range of stationary distribution, we won’t need burnin.
• If the initial value is out the range of stationary distribution, we need need to discard them.

A burnin example from my own research

Auto correlation

• MCMC chains always show autocorrelation (AC).
• AC means that adjacent samples in time are highly correlated.

\[R_x(k) = \frac {\sum_{t=1}^{n-k} (x_t - \bar{x})(x_{t+k} - \bar{x})} {\sum_{t=1}^{n-k} (x_t - \bar{x})^2} \]

• The first-order AC \(R_x(1)\) can be used to estimate the Sample Size Inflation Factor (SSIF): \[ s_x = \frac{1 + R_x(1)}{1 - R_x(1)} \]
• If we took n samples with SSIF, then the effective sample size is \(n/s_x\)

How to deal with autocorrelation?

High autocorrelation leads to smaller effective sample size.

• Design smarter algorithm to make auto-correlation smaller.
• Thining: Only take samples every 10 iterations.
• Just keep everything and make the Markov chain longer.

Converge by chain mixing

• Monitor convergence by plotting samples (of r.v.s) from multiple chains (multiple MCMC runs).
  • If the chains are well mixed (left), they are probably converged.
  • If the chains are poorly-mixed (right), we should continue burn-in.

Converge by likelihood

• How to monitor the convergence of the Markov chain.
  • Monitor the pattern of samples.
  • Monitor likelihood.

Why Gibbs is MH with \(A(x'|x) = 1\)?

• The Gibbs sampling proposal distribution is: \[q(x'_i, \textbf{x}_{-i} | x_i, \textbf{x}_{-i} ) = p(x'_i | \textbf{x}_{-i})\]
• Applying MH to this proposed distribution:

\[\begin{align*} A(x'_i, \textbf{x}_{-i} | x_i, \textbf{x}_{-i} ) &= \min \bigg( 1, \frac{p(x'_i,\textbf{x}_{-i}) q(x_i, \textbf{x}_{-i} | x'_i, \textbf{x}_{-i} ) } {p(x_i,\textbf{x}_{-i}) q(x'_i, \textbf{x}_{-i} | x_i, \textbf{x}_{-i} )} \bigg) \\ &= \min \bigg( 1, \frac{p(x'_i,\textbf{x}_{-i}) p(x_i| \textbf{x}_{-i} ) } {p(x_i,\textbf{x}_{-i}) p(x'_i, |\textbf{x}_{-i} )} \bigg) \\ &= \min \bigg( 1, \frac{p(x'_i| \textbf{x}_{-i} ) p(\textbf{x}_{-i}) p(x_i| \textbf{x}_{-i} ) } {p(x_i| \textbf{x}_{-i} ) p( \textbf{x}_{-i}) p(x'_i |\textbf{x}_{-i} )} \bigg) \\ &= \min (1,1)\\ &= 1 \end{align*}\]

• Gibbs sampling is a MH with acceptance rate 100%.

References

• Bayesian Data analysis