How can I implement Gibbs sampler for the posterior distribution, and estimating the marginal posterior distribution by making histogram?


closed as off-topic by Xi'an, mdewey, Juho Kokkala, kjetil b halvorsen, Michael Chernick Jun 7 '18 at 23:18

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – Juho Kokkala, Michael Chernick
If this question can be reworded to fit the rules in the help center, please edit the question.



Suppose $Y \sim \text{N}(\text{mean} = \mu, \text{Var} = \frac{1}{\tau})$.

Based on a sample, obtain the posterior distributions of $\mu$ and $\tau$ using the Gibbs sampler.


$ \mu$ = population mean

$ \tau$ = population precision (1 / variance)

$n$ = sample size

$\bar{y}$ = sample mean

$s^2$ = sample variance

Gibbs sampler

[Casella, G. & George, E. I. (1992). Explaining the Gibbs Sampler. The American Statistician, 46, 167–174.]

At iteration $i$ ($i = 1, \dots, N$):

  • sample $\mu^{(i)}$ from $f(\mu \,|\, \tau^{(i - 1)}, \text{data})$ (see below)
  • sample $\tau^{(i)}$ from $f(\tau \,|\, \mu^{(i)}, \text{data})$ (see below)

The theory ensures that after a sufficiently large number of iterations, $T$, the set $\{(\mu^{(𝑖)}, \tau^{(𝑖)} ) : i = T+1, \dots, 𝑁 \}$ can be seen as a random sample from the joint posterior distribution.


$f(\mu, \tau) = f(\mu) \times f(\tau)$, with

$f(\mu) \propto 1$

$f(\tau) \propto \tau^{-1}$

Conditional posterior for the mean, given the precision $$(\mu \,|\, \tau, \text{data}) \sim \text{N}\Big(\bar{y}, \frac{1}{n\tau}\Big)$$

Conditional posterior for the precision, given the mean $$(\tau \,|\, \mu, \text{data}) \sim \text{Gam}\Big(\frac{n}{2}, \frac{2}{(n-1)s^2 + n(\mu - \bar{y})^2} \Big)$$

(quick) R implementation

# summary statistics of sample
n    <- 30
ybar <- 15
s2   <- 3

# sample from the joint posterior (mu, tau | data)
mu     <- rep(NA, 11000)
tau    <- rep(NA, 11000)
T      <- 1000    # burnin
tau[1] <- 1  # initialisation
for(i in 2:11000) {   
    mu[i]  <- rnorm(n = 1, mean = ybar, sd = sqrt(1 / (n * tau[i - 1])))    
    tau[i] <- rgamma(n = 1, shape = n / 2, scale = 2 / ((n - 1) * s2 + n * (mu[i] - ybar)^2))
mu  <- mu[-(1:T)]   # remove burnin
tau <- tau[-(1:T)] # remove burnin

$$ $$


enter image description here

  • 2
    $\begingroup$ +1 This is the way to answer an apparently code-specific question: provide the theoretical explanation (which makes it on-topic here) and then (to satisfy the OP) give the code as well. Nice to see you back! $\endgroup$ – whuber Mar 10 '17 at 15:06
  • $\begingroup$ Thanks for the explanation. Is there a way I can learn coding for all statistical algorithm. I am a student and I want to learn some coding for R so that I can solve my problem. $\endgroup$ – sabuj bhowmick Mar 16 '17 at 10:56

Not the answer you're looking for? Browse other questions tagged or ask your own question.