1
$\begingroup$

What is the full-conditional distribution for $[\sigma|\textbf{y},\mu]$ given the following hierarchical structure?:

$y_i \sim N(\mu,\sigma^2)$

$\mu \sim N(\mu_0, \sigma^2_0)$

$log(\sigma) \sim N(\mu_\sigma,\tau_\sigma^2)$

My work: $[\sigma|\cdot] \propto (\sigma^2)^{-n/2}exp(\frac{1}{-2\sigma^2}\Sigma^n_{i=1}(y_i-\mu)^2)(\sigma^2)^{-1/2}exp(\frac{-(ln(\sigma)-\mu_\sigma)^2}{2\tau^2_\sigma})$

$\propto (\sigma^2)^{-(\frac{n-1}{2}+1)}exp[\frac{1}{-2\sigma^2}\Sigma^n_{i=1}(y_i-\mu)^2-\frac{1}{2\tau^2_\sigma}(ln^2(\sigma)-2ln(\sigma)\mu_\sigma)]$

From here, I recognize that this somewhat resembles an inverse gamma distribution, where you let $q=\frac{n-1}{2}$. However, I need to find the $r$. I know that the $exp$ term should follow the format $exp[-\frac{1}{r\sigma^2}]$, but I do not know how to manipulate my last line to obtain that $r$. Do you have any suggestions?

Edit:

I now am getting

$[\sigma|\cdot] \propto (\sigma)^{-n-1+\frac{\mu_\sigma}{2\tau_\sigma^2}}exp(\frac{1}{-2\sigma^2}\Sigma^n_{i=1}(y_i-\mu)^2-\frac{ln(\sigma)^2}{2\tau_\sigma^2})$

This more closely resembles an inverse gamma distribution, but I still have the last pesky exponential term that is getting in the way of this being the IG kernel.

$\endgroup$
  • 1
    $\begingroup$ Is there a reason why you cannot use the conjugate prior setup where the variance is a priori inverse-gamma? This is a standard setup and the posterior distribution of the variance will also be inverse-gamma. $\endgroup$ – bayes Feb 7 at 12:10
  • $\begingroup$ @bayes I have used that conjugate setup. 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. I have been trying to obtain a recognizable kernel from the posterior for a while, and I cannot figure it out. Yet he stated that it is possible to implement the posterior with an MCMC algorithm. Is it possible to implement this hierarchical model in an MCMC algorithm without knowing the posterior distribution? Do I have to use M-H ratio or Gibbs updates? $\endgroup$ – Ron Snow Feb 7 at 17:50
  • 1
    $\begingroup$ What makes you think the kernel must be recognisable? $\endgroup$ – Earlien Feb 9 at 23:29
  • $\begingroup$ @Earlien I would assume that it had to be recognizable, since I need to implement the lognormal prior for $\sigma$ in an MCMC algorithm. Is there a way where I can implement this hierarchical model in an MCMC algorithm without having to identify the full-conditional distribution? $\endgroup$ – Ron Snow Feb 12 at 2:17
  • $\begingroup$ In most real-world applications, conjugacy is not tenable or preferred, and thus most full conditionals are not recognisable. This isn't a problem. This is why have MCMC methods - to sample from tricky non-recognisable kernals. Metropolis-Hastings, slice sampling, ARS, Hamiltonian MC - whatever algorithm you want. Being recognisable is helpful and simplifies the procedure to simple Monte Carlo samples, but is not required. You have already "identified" the full conditional. Now you need to do the sampling. $\endgroup$ – Earlien Feb 12 at 5:42
2
+50
$\begingroup$

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.


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

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ No, that would require multiplying by the prior for $\mu$ as well. This is the full conditional for $\log\sigma$. $\endgroup$ – Earlien Feb 12 at 8:01
  • 1
    $\begingroup$ No, that is constant with respect to the parameter in question and thus irrelevant. Perhaps I should use “proportional to” on the first line. $\endgroup$ – Earlien Feb 12 at 9:26
  • 1
    $\begingroup$ Just a few follow-up questions: So, I would be obtaining posterior updates for $log(\sigma)$. Then I simply take the exponent to obtain the posterior updates for $\sigma$, right? Also, why is it so much harder to go down the path I was pursuing than yours? I figure that they would both lead to the same answer. I plan on using a Gibbs loop to obtain the posterior estimates for the other priors. I should have this implemented by the weekend, and I will update you with me progress. $\endgroup$ – Ron Snow Feb 13 at 6:13
  • 1
    $\begingroup$ Yes, to get the posterior for $\sigma$, you simply exponentiate each value in the MCMC chain for $\log\sigma$. To find the posterior for $\sigma$ analytically, you need to account for the Jacobian since you are transforming, that is, $\sigma$ has a log-Normal distribution. It would be like finding the posterior for $\sigma^2$ analytically. Its simpler to work with the parameter that you have the prior on rather than a transformation of it. Good luck with the sampling. $\endgroup$ – Earlien Feb 13 at 21:52
  • 1
    $\begingroup$ @Edison I think there is some misunderstanding here. Most real world examples will not be conjugate (and conjugate isn't always desirable). This is one example. MCMC sampling was designed specifically to sample from these unrecognisable distributions (as an extension of Monte Carlo (MC) methods for conditional distributions). You could code an MCMC algorithm from scratch to implement this model. But why not use freely available software like BUGS or JAGS? If by "work with the distribution" you mean estimate the posterior, I don't understand your aversion. $\endgroup$ – Earlien Feb 19 at 23:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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