Exact Sampling from Improper Mixtures

Question

Suppose I want to sample from a continuous distribution $p(x)$. If I have an expression of $p$ in the form

$$p(x) = \sum_{i=1}^\infty a_i f_i(x)$$

where $a_i \geqslant 0, \sum_i a_i= 1$, and $f_i$ are distributions which can easily be sampled from, then I can easily generate samples from $p$ by:

1. Sampling a label $i$ with probability $a_i$
2. Sampling $X \sim f_i$

Is it possible to generalise this procedure if the $a_i$ are occasionally negative? I suspect I've seen this done somewhere - possibly in a book, possibly for the Kolmogorov distribution - so I'd be perfectly happy to accept a reference as an answer.

If a concrete toy example is helpful, let's say I'd like to sample from $$p(x,y) \propto \exp(-x-y-\alpha\sqrt{xy})\qquad x,y > 0$$ I'll then take $\alpha \in (0, 2)$ for technical reasons which should not matter too much, in the grand scheme of things.

In principle, I could then expand this as the following sum:

$$p(x,y) \propto \sum_{n=0}^\infty \frac{(-1)^n \alpha^n \left( \frac{n}{2} \right)! \left( \frac{n}{2} \right)!}{n!} \left( \frac{x^{n/2} e^{-x}}{\left( \frac{n}{2} \right)!}\right) \left( \frac{y^{n/2} e^{-y}}{\left( \frac{n}{2} \right)!}\right) .$$

The $(x,y)$-terms inside the sum can then be independently sampled from as Gamma random variates. My issue is evidently that the coefficients are "occasionally" negative.

Edit 1: I clarify that I am seeking to generate exact samples from $p$, rather than calculating expectations under $p$. For those interested, some procedures for doing so are alluded to in the comments.

Edit 2: I found the reference which includes a particular approach to this problem, in Devroye's 'Non-Uniform Random Variate Generation'. The algorithm is from 'A Note on Sampling from Combinations of Distributions', of Bignami and de Matteis. The method is effectively to bound the density from above by the positive terms of the sum, and then use rejection sampling based on this envelope. This corresponds to the method described in @Xi'an's answer.

Answer 1

I have puzzled over this question but never came with a satisfying solution.

One property that is of possible use is that, if a density writes $$f(x)=\frac{g(x)-\omega h(x)}{1-\omega}\qquad \omega>0$$ where $g$ is a density such that $g(x)\ge \omega h(x)$, simulating from $g$ and rejecting these simulations with probability $\omega h(x)/g(x)$ delivers simulations from $f$. In the current case, $g$ is the normalised version of the positive weight components $$g(x)=\sum_{\alpha_i>0} \alpha_i f_i(x) \big/ \sum_{\alpha_i>0} \alpha_i$$ and $\omega h$ is the remainder $$h(x)=\sum_{\alpha_i<0} \alpha_i f_i(x) \big/ \sum_{\alpha_i<0}\alpha_i$$ This is indeed found in the simulation bible of Devroye, Non-uniform random variate generation, Section II.7.4, but follows from a simple accept-reject reasoning.

A first computational drawback of this approach is that, despite simulating first from a chosen component $f_i$, the sums in both $g$ and $h$ must be computed for the rejection step. If the sums are infinite with no closed form version, this makes the accept-reject method impossible to implement.

A second difficulty is that, since both sums of weights are of the same order$$\sum_{\alpha_i>0}\alpha_i = 1 - \sum_{\alpha_i<0}\alpha_i$$ the rejection rate $$1-\varrho^\text{accept}=\sum_{\alpha_i<0}|\alpha_i| \Big/ \sum_i |\alpha_i|$$has no upper limit. Actually if the series associated with the $\alpha_i$'s is not absolutely converging, the acceptance probability is zero! And the method cannot be implemented in this situation.

In the case of a mixture representation, if $f$ can be written as $$f(x)=\sum_{i=1}^\infty \alpha_i \frac{g_i(x)-\omega_i h(x_i)}{1-\omega_i}\qquad \omega_i>0$$the component can be chosen first and then the method applied to the component. But this may be delicate to implement, identifying pairs $(g_i,h_i)$ that fit $g_i(x)-\omega_i h(x_i)>0$ from the possibly infinite sum being not necessarily feasible.

I think a more efficient resolution could come from the series representation itself. Devroye, Non-uniform random variate generation, Section IV.5, contains a large range of series methods. As for instance the following algorithm for an alternate series representation of the target $$f(x)=\kappa h(x)\{1-a_1(x)+a_2(x)-\cdots\}$$ when the $a_i(x)$'s converge to zero with $n$ and $h$ is a density:

The problem has been considered recently in the context of debiasing biased estimators for MCMC, as for instance in the Glynn-Rhee approach. And the Russian roulette estimator (with a connection with the Bernoulli factory problem). And the unbiased MCMC methodology. But there is no escape from the sign issue... Which makes its use challenging when estimating densities as in pseudo-marginal methods.

Upon further thinking, my conclusion is that there is no generic method to produce an actual simulation from this series [rather than mixture which turns out to be a misnomer], without imposing further > structure to the elements of the series, like the one in the above algorithm from Devroye's bible. Indeed, since most (?) densities allow for a series expansion of the kind above, this would otherwise imply the existence of a sort of universal simulation machine...

Answer 2

I have the draft of an idea that could work. It is not exact, but hopefully asymptotically exact. To turn it into a really rigorous method, where the approximation is controlled, or something about it can be proven, there's probably a lot of work needed.

First, as mentioned by Xi'an you can group the positive weights on one hand and the negative weights on the other hand, so that that finally the problem has only two distributions $g$ and $h$:

$$p=\lambda g - \mu h$$

with $\lambda-\mu=1$. Note that you have $\lambda\geq 1$.

My idea is the following. You want sample $N$ observations from $p$. Do:

• sample $\lambda N$ values from $g$ and store them in a list
• for each of $\mu N$ values sampled from $h$, remove their closest (remaining) neighbour from the list.

At the end you get $(\lambda-\mu)N=N$ points. It does not need to be exactly the closest neighbour, but just a point that is "close enough". The first step is like generating matter. The second step is like generating antimatter and let it collide and cancel with matter. This method is not exact, but I believe, under some conditions, it is asymptotically exact for large $N$ (to make it almost exact for small $n$ you need to use a large $N$ first and then take a small random part of the final list). I'm giving a very informal argument that is more an explanation than a proof.

Consider $x$ in the observation space and a small volume $v$ around $x$ with Lebesgue volume $\epsilon$. After sampling from $g$, the number of elements in the list that are also in $v$ is approximatively $\lambda Ng(x)\epsilon$. After the second step, approximatively $\mu Nh(x)\epsilon$ will be removed from it, and you have approximatively the desired number $Np(x)\epsilon$. For this you need to assume that the number of points in the volume is sufficiently large.

This method is very unlikely to resist to large dimension or some pathologies of $g$ and $h$ but might work in small dimension and sufficiently smooth, "sufficiently uniform" distributions.

Note about an exact method:

I first thought of this for discrete distributions, and clearly in that case the method is not exact, since it can generate samples that have probability 0. I have the strong intuition that an exact method is not possible with finite processing time, and that this impossibility could be proven, at least for discrete distributions. The rule of the game is that you are only allowed to use exact "oracle" samplers for $g$ and $h$ but don't know $g$ and $h$ as functions of $x$. For simplicity restrict to Bernoulli distributions. The non existence of an exact method is related to Bernoulli Factory theory: if you could create a $(\lambda p - \mu q)$-coin from a $p$-coin and a $q$-coin, then you could create a $\lambda p$-coin from a $p$-coin which is known to be impossible for $\lambda>1$.