Skip to main content
Cross Validated

Return to Answer

Added link to external file for further reading.
Source Link
Earlien
Earlien
  • 758
  • 6
  • 22

Some additional commentary on the practicality of obtaining a posterior estimate:

  1. Most real world examples will not have conjugate priors (and conjugate isn't always desirable). Consequently, the full conditionals will not always be recognisable. That is, you cannot expect to compare the resulting full conditional with the kernel of another recognisable distribution and expect the parameters to align. If they did, this would certainly make sampling easier. This would allow Gibbs sampling, and by that I mean, you can essentially use Monte Carlo methods (direct sampling from the distribution). Most software like R has capability to draw from such distributions, e.g. rnorm, rbinom, rgamma etc.

  2. What happens when one or more full conditionals are not recognisable? Monte Carlo methods (like naive accept-reject methods) can also sample from these unfamiliar kernels. The problem lies in the fact that the full conditionals are not independent. This problem was first solved in the 1950s with the advent of the Metropolis algorithm (and with the later extension to the Metropolis-Hastings algorithm in 1970). These algorithms extend Monte Carlo methods by using Markov Chain theory to sample from inter-dependent distributions, hence the name, Markov Chain Monte Carlo. More modern algorithms include Hamiltonian Monte Carlo (1987; extended 2014), adaptive rejection-sampling (1992; extended 1995; 2010) and slice samplimg (1997; 2003) to name a few. Sequential Monte Carlo is another, slightly different method.

  3. So how do you obtain a posterior in practice? First, to recap, your model specified in terms of the two full conditionals are:

$$ \begin{align} p(\mu | \boldsymbol{y}, \sigma) & \propto \frac{1}{\sigma_0\sigma^n}\exp\left[-\frac{1}{2}\left(\frac{(\mu - \mu_0 )^2}{\sigma_0^2} +\frac{\sum_{i=1}^n(y_i-\mu )^2}{\sigma^2}\right)\right]\\ p(\log\sigma | \boldsymbol{y}, \mu) & \propto \frac{1}{\sigma^n}\exp\left[-\frac{1}{2}\left(\frac{(\log\sigma-\mu_\sigma )^2}{\tau_\sigma^2} +\frac{\sum_{i=1}^n(y_i-\mu )^2}{\sigma^2}\right)\right] \end{align} $$

You could code an MCMC algorithm from scratch to implement this model: you sample from $p(\mu | \boldsymbol{y}, \sigma)$ using an initial value for $\sigma$ then sample from $p(\log\sigma | \boldsymbol{y}, \mu)$ using an initial value for $\mu$. And repeat. As for what particular algorithm you use is up to you. E.g. you could use Metropolis-Hasting within Gibbs or Hamiltonian within Gibbs, etc. (The "within Gibbs" expression simply means you are making use of the full conditionals, as opposed to performing Metropolis-Hastings or Hamiltonian on the joint distribution $p(\mu, \sigma | \boldsymbol{y})$. Using the joint distribution rather than the full conditionals is also possible, but is often too slow computationally, especially for sparse, high-dimensional models, hence the usefulness of Gibbs sampling).

You don't have to implement an MCMC algorithm from scratch though (unless your adviser has specifically requested that you do). There are several freely available software that can easily implement a model like this, for example, WinBUGS/OpenBUGS, JAGS, and Stan. This will provide you with a posterior estimate of the two parameters. You can then compare the results via sample statistics or visually using plots of the posterior to alternate model specifications (like your inverse-gamma model).

For more about deriving the full conditionals and why having full conditionals that are unrecognisable is not a "problem", please see the slides from a presentation I gave a few years ago: https://bragqut.files.wordpress.com/2018/04/deriving-the-full-conditionals.pdf

Some additional commentary on the practicality of obtaining a posterior estimate:

  1. Most real world examples will not have conjugate priors (and conjugate isn't always desirable). Consequently, the full conditionals will not always be recognisable. That is, you cannot expect to compare the resulting full conditional with the kernel of another recognisable distribution and expect the parameters to align. If they did, this would certainly make sampling easier. This would allow Gibbs sampling, and by that I mean, you can essentially use Monte Carlo methods (direct sampling from the distribution). Most software like R has capability to draw from such distributions, e.g. rnorm, rbinom, rgamma etc.

  2. What happens when one or more full conditionals are not recognisable? Monte Carlo methods (like naive accept-reject methods) can also sample from these unfamiliar kernels. The problem lies in the fact that the full conditionals are not independent. This problem was first solved in the 1950s with the advent of the Metropolis algorithm (and with the later extension to the Metropolis-Hastings algorithm in 1970). These algorithms extend Monte Carlo methods by using Markov Chain theory to sample from inter-dependent distributions, hence the name, Markov Chain Monte Carlo. More modern algorithms include Hamiltonian Monte Carlo (1987; extended 2014), adaptive rejection-sampling (1992; extended 1995; 2010) and slice samplimg (1997; 2003) to name a few. Sequential Monte Carlo is another, slightly different method.

  3. So how do you obtain a posterior in practice? First, to recap, your model specified in terms of the two full conditionals are:

$$ \begin{align} p(\mu | \boldsymbol{y}, \sigma) & \propto \frac{1}{\sigma_0\sigma^n}\exp\left[-\frac{1}{2}\left(\frac{(\mu - \mu_0 )^2}{\sigma_0^2} +\frac{\sum_{i=1}^n(y_i-\mu )^2}{\sigma^2}\right)\right]\\ p(\log\sigma | \boldsymbol{y}, \mu) & \propto \frac{1}{\sigma^n}\exp\left[-\frac{1}{2}\left(\frac{(\log\sigma-\mu_\sigma )^2}{\tau_\sigma^2} +\frac{\sum_{i=1}^n(y_i-\mu )^2}{\sigma^2}\right)\right] \end{align} $$

You could code an MCMC algorithm from scratch to implement this model: you sample from $p(\mu | \boldsymbol{y}, \sigma)$ using an initial value for $\sigma$ then sample from $p(\log\sigma | \boldsymbol{y}, \mu)$ using an initial value for $\mu$. And repeat. As for what particular algorithm you use is up to you. E.g. you could use Metropolis-Hasting within Gibbs or Hamiltonian within Gibbs, etc. (The "within Gibbs" expression simply means you are making use of the full conditionals, as opposed to performing Metropolis-Hastings or Hamiltonian on the joint distribution $p(\mu, \sigma | \boldsymbol{y})$. Using the joint distribution rather than the full conditionals is also possible, but is often too slow computationally, especially for sparse, high-dimensional models, hence the usefulness of Gibbs sampling).

You don't have to implement an MCMC algorithm from scratch though (unless your adviser has specifically requested that you do). There are several freely available software that can easily implement a model like this, for example, WinBUGS/OpenBUGS, JAGS, and Stan. This will provide you with a posterior estimate of the two parameters. You can then compare the results via sample statistics or visually using plots of the posterior to alternate model specifications (like your inverse-gamma model).

For more about deriving the full conditionals and why having full conditionals that are unrecognisable is not a "problem", please see the slides from a presentation I gave a few years ago: https://bragqut.files.wordpress.com/2018/04/deriving-the-full-conditionals.pdf

Bounty Ended with 50 reputation awarded by Ron Snow
Source Link
Earlien
Earlien
  • 758
  • 6
  • 22

Initially I wondered why you wanted to derive the full conditional for $\sigma$. But in the comments you say:

My research adviser, in preparation for working with him, has given me the task to compare the a priori inverse-gamma with the log normal stated in the post.

To this end, would it suffice to simply have the two analytical forms of the priors? Perhaps the end goal is actually comparing the posteriors though. Your adviser also said

...it is possible to implement the [model to determine the] posterior with an MCMC algorithm

which is true. In which case, I suggesting working with the full-conditional for $\log\sigma$ since this is how you specify it in the model. Once you have a posterior estimate for $\log\sigma$, it is then trivial to get a posterior estimate fro $\sigma$ which you can then compare to your results from your atlernate model with the inverse-gamma prior.

Assuming that $\mu_0$, $\sigma_0^2$, $\mu_\sigma$, and $\tau_\sigma^2$ are known a priori (since you don't give them priors), then the full conditional for $\log\sigma$ is:

$$ \begin{align} p(\log\sigma | \boldsymbol{y}, \mu) & = p(\log\sigma) \times p(\boldsymbol{y} | \mu, \sigma^2) \\ & = \frac{1}{\sqrt{2\pi \tau_\sigma^2}}\exp\left(-\frac{1}{2}\frac{(\log\sigma-\mu_\sigma )^2}{\tau_\sigma^2} \right) \times \prod_{i=1}^n \frac{1}{\sqrt{2\pi \sigma^2}}\exp\left(-\frac{1}{2}\frac{(y_i-\mu )^2}{\sigma^2} \right) \\ & \propto \exp\left(-\frac{1}{2}\frac{(\log\sigma-\mu_\sigma )^2}{\tau_\sigma^2} \right) \times \frac{1}{\sigma^n} \exp\left(-\frac{1}{2}\sum_{i=1}^n\frac{(y_i-\mu )^2}{\sigma^2} \right) \\ & = \frac{1}{\sigma^n}\exp\left[-\frac{1}{2}\left(\frac{(\log\sigma-\mu_\sigma )^2}{\tau_\sigma^2} +\frac{\sum_{i=1}^n(y_i-\mu )^2}{\sigma^2}\right)\right] \end{align} $$

Together with the other two full conditionals, you could implement this using an MCMC algorithm (e.g. Metropolis Hastings, adaptive rejection sampling, Hamiltonian Monte Carlo, slice sampling, etc.)

If you are only interested in obtaining an MCMC posterior estimate, then you could implement this in software like BUGS or JAGS which doesn't even require you to determine the full conditionals.