1

I'm using the Metorpolis-Hastings algorithm in a setting where the acceptance function is essentially of the form $$\alpha(x,y)=1\wedge\frac{u(x,y)}{v(x,y)},$$ where $$u(x,y)=p+(1-p)\prod_{i=1}^mu_i(x,y)$$ and $$v(x,y)=p+(1-p)\prod_{j=1}^nv_j(x,y).$$

My problem is that, if $m$ or $n$ is large enough (in my application less than $15$ is sufficient), it may happen that the product in the definition of $u(x,y)$ or $v(x,y)$ is computed to $0$ (due to insufficient floating point precision) even when each $u_i(x,y)$ or $v_i(x,y)$ is strictly positive (each $u_i(x,y)$ or $v_i(x,y)$ is essentially the value of a normal distribution density).

What can I do to remedy this issue?

EDIT: The concrete shape of $u$ and $v$ is $u(x,y)=w_{d_1}(\varphi(x),\varphi(y))$ and $v(x,y)=w_{d_2}(\psi(x),\psi(y))$, where $\varphi$ and $\psi$ are transformations onto $[0,1)^{d_1}$ and $[0,1)^{d_2}$, respectively, and $$w_d(x',y'):=p+\frac{1-p}{\sqrt{2\pi\sigma^2}}\prod_{i=1}^d\sum_{k\in\mathbb Z}e^{-\frac{\left(k+y'_i-x'_i\right)^2}{2\sigma^2}}.$$

0xbadf00d
  • 539
  • 3
  • 8
  • 5
    Use the log-exp-sum trick, there are many posts on this site dealing with the issue. – hejseb Mar 13 '20 at 18:45
  • 1
    What is the issue? The underflow of a product does not change the precision with which $\alpha$ is computed and for this purpose you need very little precision anyway--maybe three decimal places would suffice. – whuber Mar 13 '20 at 19:22
  • @whuber The issue is that both $u(x,y)$ and $v(x,y)$ might be very small (and hence being computed to $0$), while the quotient $u(x,y)/v(x,y)$ is not. – 0xbadf00d Mar 13 '20 at 19:26
  • 1
    When they are "very small," then unless $p$ itself is comparably small, it makes no difference. The only values of $p$ that would make a real difference would have to be less than $10^{-300}.$ – whuber Mar 13 '20 at 19:29
  • @whuber I don't get what you mean. For example: $1\text{e-}100/1\text{e-}99=0.1$ which should be significantly large to matter. (In any case, let me note that I've omitted the multiplication of the target distribution density in the nominator and denominator of the quotient. $u$ and $v$ only correspond to the density of the proposal.) – 0xbadf00d Mar 13 '20 at 19:35
  • 1
    There's a mismatch between your example and the formulas in your question: in the question, each of $u$ and $v$ starts with the value $p,$ to which another term is added. You are concerned that the additional term is rounded down to zero. That makes no difference unless the amount of rounding is appreciable relative to the size of $p.$ However, no additive term comparable to $p$ appears in your comment: what gives? – whuber Mar 13 '20 at 19:47
  • @whuber I see your point (while I'm uncertain whether I might forget about an important detail in my application which might still turn this into a problem even when $p\gg0$), but please assume that the solution should work for the choice $p=0$ as well (any value in $[0,1]$ should be possible for $p$). – 0xbadf00d Mar 13 '20 at 20:00
  • @hejseb Maybe I'm missing something, but I don't see how I would need to apply this in my situation. Please take note of my edit. I've described the actual shape of $u$ and $v$ in more detail. – 0xbadf00d Mar 13 '20 at 20:07
  • 3
    If $p>0$, you have no problem for the reason given by whuber. If $p=0$, compute the products instead as sums in log-space, then subtract one from the other and take the exp() to get your $\alpha$. – Creosote Mar 13 '20 at 21:27
  • @Creosote (a) It is still a problem when $p$ is close to $0$, isn't it? (b) The estimator (I'm using the [importance sampling estimator](https://arxiv.org/abs/1805.07174), still requires the division by $v(x,y)$. How should I handle that? – 0xbadf00d Mar 14 '20 at 07:55
  • 1
    Potential issues: (1) the product can go to zero when you work in floating point; or (2) $p$ and the product are of greatly different size, leading to loss of precision in the smaller when added to the larger. Case (2) does not affect $\alpha$ significantly. Case (1) affects $\alpha$ only if $p=0$, for which log-space is the answer. – Creosote Mar 14 '20 at 08:24
  • We would need much more details to understand how $u(\cdot,\cdot)$ and $v(\cdot,\cdot)$ can be actual densities, while enjoying a constant part $p$ that integrates to infinity over an unbounded space. – Xi'an Mar 14 '20 at 08:41
  • @Xi'an They are integrated with respect to the Lebesgue measure on $[0,1)^d$. So, the constant is totally fine. – 0xbadf00d Mar 14 '20 at 09:48
  • In that case, I completely second previous comments that having $u(x,y)$ equal to zero or almost to zero should not matter in the least. – Xi'an Mar 14 '20 at 10:28

0 Answers0