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?

$\endgroup$
  • 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 '19 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 '19 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 '19 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 '19 at 14:54
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.

$\endgroup$

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