How to choose a constant for reject sampling

Question

When using a non-Markov Monte Carlo sampling method, for example acceptance-rejection sampling, we choose a density $\ h(x) $ and a known constant $\ c $ such that $\ ch(x) $ acts as a blanketing function for our target distribution $\pi(x) $.

Since we just have access to the target distribution up to the unknown normalizing constant, we don´t really know what the normalized target distribution looks like. So how is it ever possible to define a blanket function to cover all values? The only examples I have seen so far were basic simulations where the analytical form of $\pi(x) $ was already known.

Answer 1

Let $\pi(x) = M f(x)$, where $M$ is the normalizing constant. In many situations, only $f(x)$ is known and $M$ is unknown.

To implement rejection sampling, you want $c$ such that, for all $x$, $$\dfrac{\pi(x)}{h(x)} \leq c. $$

Then for all $x$,

$$\dfrac{f(x)}{h(x)} \leq \dfrac{c}{M} := c'. $$

You don't know $c$ or $M$, but you should be able to find $c'$, if you can play around with $f/h$. (This is more difficult to do for higher dimensional distributions)

The algorithm accepts when $\pi(x)/ch(x)$ is greater than a realization from a uniform random variable, which is the same as $f(x)/c'h(x)$ being greater than the realization. Thus, the algorithm can be implemented even though $M$ is not known.

Answer 2

Just as an example consider the density proportional to $e^{-(1+x^4)^\frac14},\quad-\infty<x<\infty \,.$

I currently have no idea what the normalizing constant is for that density. (Well I can tell it's not going to be terribly far from 1, so I have some notion, but being able to guess it to within an order of magnitude isn't much help)

It's not hard to see that $(1+x^4)^\frac14 > x$ for $x\geq 0$ (if it's not obvious, look at the fourth power of both sides, from which it is surely clear), and so by symmetry $(1+x^4)^\frac14 > |x|$ everywhere, and hence it's clear that $e^{-|x|}>e^{-(1+x^4)^\frac14}$ on the real line.

Consequently I can use the standard Laplace distribution scaled up by a factor of 2 as a majorizing function. (There are somewhat more efficient choices, but it suffices to use this simple case. Note that we now know for sure that the integral of that unscaled density must be less than 2.)

I still don't know quite what the density in question integrates to, but I can simulate from it.

Indeed, here's a histogram of (just over) a million values from it:

(And having just run code to simulate those values, I now have some idea what the normalizing constant is, since you can work an approximation out from the acceptance rate. The integral of the unscaled density comes out to about 1.397 -- checking that with the integrate function in R, it gives 1.396785 with absolute error <0.00011 -- so the normalizing constant would approximately be the reciprocal of that.)

Answer 3

More generally a principle to choose $M$ is

$$\inf_{\theta\in\Theta}\sup_{x\in\mathbb{R}}\frac{f(x)}{g_\theta(x)}$$

where $f$ is the normalized target and $g_\theta$ is the proposal density.

As an example, for $f(x)=(2\pi)^{-\frac{1}{2}}e^{-x^2/2}$ and $g(x)=\theta^{-1}e^{-\theta |x|}$

We will have $\sup_xf(x)/g_\theta(x)=\theta(2\pi)^{-\frac{1}{2}}e^{-\theta^2/2}$ and $\inf_\theta\theta(2\pi)^{-\frac{1}{2}}e^{-\theta^2/2}=(2\pi)^{-\frac{1}{2}}e^{-\frac{1}{2}}$

But if your target is unnormalized then you can only try by luck. (correct me if I am wrong)