3
$\begingroup$

Let's say I have $n$ observations of a random variable, $X_1, \dotsm, X_n \sim \mathcal{N}(0, \sigma^2)$. I also assume $\sigma^2$ has a Gamma(1,1) prior distribution, $\pi(x) = \exp(-x)$.

I'm now attempting to use Metropolis-Hastings to sample from the posterior distribution (which I believe is):

$$f(\sigma | X) \propto \frac{1}{(2\pi\sigma^2)^{n/2}}\exp{-\left( \sigma^2 + \sum\limits_{i=1}^n\frac{X_i^2}{2\sigma^2}\right)}$$

However, for larger $n$, this unnormalized density produces usually either quite small or quite large values, making it difficult to get the Markov chain to mix well.

My question is: in general (and for the MH algorithm), what are my options for attempting to sample with such an unwieldy distribution?

Improve this question
$\endgroup$
5
  • 3
    $\begingroup$ This is a generalised inverse Gaussian distribution. Not particularly unwieldy. If the issue is numerical (underflows and overflows), then work (a) with log-densities and (b) with differences from the maximal value of the log-density. $\endgroup$
    – Xi'an
    Dec 9, 2019 at 17:17
  • $\begingroup$ @Xi'an - does using the log-density not break any of the assumptions of MH or other algorithms used for sampling from densities? And can one just take the exponential of the samples obtained via MH to get transform them back to samples from the original density? $\endgroup$
    – CLL
    Dec 10, 2019 at 9:49
  • $\begingroup$ And good point - I guess by "unwieldy", I meant more unwieldy than most distributions that I've had to deal with. In my case, the $\sum X_i^2$ term is quite large, as is $n$, and it's made for a lot of underflows and overflows $\endgroup$
    – CLL
    Dec 10, 2019 at 9:51
  • 1
    $\begingroup$ Please see stats.stackexchange.com/… for threads related to computing with numbers that might overflow or underflow a likelihood or probability. The current reply by @Xi'an nicely illustrates a standard procedure. $\endgroup$
    – whuber ♦
    Dec 10, 2019 at 14:54
  • $\begingroup$ out of curiosity, why use a gamma prior when inverse gamma is conjugate? If you desire a gamma variate, why not parameterize the normal in terms of precision? $\endgroup$
    – Faydey
    Sep 7, 2021 at 18:06

1 Answer 1

Reset to default
6
$\begingroup$

does using the log-density not break any of the assumptions of MH or other algorithms used for sampling from densities? And can one just take the exponential of the samples obtained via MH to get transform them back to samples from the original density?

The issue with underflows and upperflows can be treated through logarithms to some extent, without jeopardising the validity of the Metropolis algorithm. The acceptance step $$u_t\le\dfrac{f(y_{t})q(y_t,x_{t-1})}{f(x_{t-1})q(x_{t-1},y_t)}$$ can be replaced with $$\log u_t\le \log\dfrac{f(y_{t})q(y_t,x_{t-1})}{f(x_{t-1})q(x_{t-1},y_t)}$$ and \begin{align*} \log\dfrac{f(y_{t})q(y_t,x_{t-1})}{f(x_{t-1})q(x_{t-1},y_t)} &= \log f(y_{t})-\log f(x_{t-1})+\log q(y_t,x_{t-1})-\log q(x_{t-1},y_t)\\ &= \{\log f(y_{t}) -\max_z \log f(z)\}-\{\log f(x_{t-1})-\max_z \log f(z)\}\\ &\qquad\qquad +\log q(y_t,x_{t-1})-\log q(x_{t-1},y_t)\\ \end{align*} which assuming $\max_z \log f(z)$ does not produce an overflow means a wider range of values of $y_t$ can be explored.

Here is a (codegolfed) rendition in R:

d<-function(y,x,S=1e6,n=1e4)n*log(y^2/x^2)/2+y^2-x^2+S/y^2/2-S/x^2/2
m<-function(T,y=rnorm(1))ifelse(rep(T>1,T),
   c(y*{d({z<-m(T-1)}[1],y+z[1])<rexp(1)}+z[1],z),y)

where the exponential generator appears because $-\log(u_t)$ is an exponential variate.

Improve this answer
$\endgroup$

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