Gibbs sampling to produce posterior pdf

Question

Suppose we have the following classical normal linear regression model:

$$y_i = \beta_1 x_{1i} + \beta_2x_{2i} + \beta_3x_{3i} + e_i$$

where $e_{i} \sim iid.N(0, \sigma^2)$ for all $i = 1, 2, \cdots, n$ and $x_{1i} = 1$ for all $i = 1, 2, \cdots, n$.

Assume that we have known data values for both $x_{2i}$ and $x_{3i}$ for all $i = 1, 2, \cdots, n$. Defining $\boldsymbol{\beta} = (\beta_1, \beta_2, \beta_3)'$ and assuming a non-informative prior of the form $p(\boldsymbol{\beta}, \sigma) \propto \frac{1}{\sigma}$, then we can show that the conditional posterior pdf for $\boldsymbol{\beta}$ and $\sigma$, that is, $p(\boldsymbol{\beta}\mid\sigma, \mathbf{y})$ and $p(\sigma\mid\boldsymbol{\beta},\mathbf{y})$ are normal and inverted gamma, respectively.

The question is: Use a Gibbs Sampler and estimate the posterior pdf of the parameter function: $\displaystyle \psi = \frac{\beta_2 + \beta_3}{\sigma^2}$.

Now I have run a Gibbs sampler (in R) and after a burn in period of 100 draws, I have obtained 1000 draws of $(\boldsymbol{\beta}^{(i)}, \sigma^{(i)})$, that is, I have a sample of $(\beta_1^{(i)}, \beta_2^{(i)}, \beta_3^{(i)}, \sigma^{(i)})$ for $i = 1, 2, \cdots, 1000$, how can I use these draws to produce an estimate of the posterior pdf of $\psi$? In other words, how can I estimate $p(\psi\mid\mathbf{y})$?

EDIT: I'm new to MCMC logarithms. I do understand what you mean but I am still not too sure how to use it in the context of this question. From what I've learnt so far, say we have $\boldsymbol{\theta} = (\theta_1, \theta_2)'$ and $p(\boldsymbol{\theta}\mid\mathbf{y}) = p(\theta_1, \theta_2\mid\boldsymbol{y})$ is the joint posterior, then the marginal posterior of $\theta_1$ is given by $p(\theta_1\mid\mathbf{y}) = \int_{\theta_1} p(\theta_1\mid\theta_2,\mathbf{y})p(\theta_2\mid\mathbf{y})d\theta_2$, now say I have a sample of $M$ draws of $(\theta_1^{(i)}, \theta_2^{(i)})$ from $p(\theta_1, \theta_2\mid\boldsymbol{y})$, then to estimate $p(\theta_1\mid\mathbf{y})$, we use the following sample mean: $\widehat{p(\theta_1\mid\mathbf{y})} = \frac{1}{M} \sum_{i=1}^M p(\theta_1\mid\theta_2^{(i)}, \mathbf{y}) $, that is, we estimate the marginal density by averaging over the conditional densities. How can I implement that here?

Answer 1

Perhaps I've misunderstood how it works, but if you have your samples, then I think you can just create $\psi$ for each $i$. For example, on the Wikipedia page for Gibbs sampling it says:

The samples then approximate the joint distribution of all variables. Furthermore, the marginal distribution of any subset of variables can be approximated by simply examining the samples for that subset of variables, ignoring the rest.

Here's a short R example to illustrate what I mean. Instead of the random number generation in the beginning you could just replace that with your draws.

beta1 <- rpois(1000, 3)
beta2 <- rnorm(1000, 2, beta1 + 1)
sigma2 <- runif(1000, 2, 3)

psi <- (beta1+beta2)/sigma2
plot(psi, type = "n", xlim = c(-5, 10), ylim = c(0, 0.4), xlab = "Density", ylab = "Psi")
lines(density(psi), col = "red", lwd = 3)