# Exact Sampling from Improper Mixtures

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.

• Why can't you sample by just using the absolute value of $a_i$, and then negating your $X\sim f_i$ sample? In other words define $Z:=\sum_{i=1}^\infty |a_i|$ (assuming it's finite), and then renormalize your sum by $Z$. Feb 1 '18 at 21:35
• @AlexR. If I understand you, a version of this would be practical for calculating expectations under $p$, but still not for drawing exact samples from $p$. Certainly this is an answer to a relevant problem, though not quite what I'm looking for.
– πr8
Feb 1 '18 at 21:37
• It depends on what you intend to do with that sample. For the purpose of computing moments, for instance, it looks straightforward to generalize sampling from mixtures of densities by additionally flagging any point selected from a component with negative coefficient as being a "negative" point and weighting its contribution negatively in the moment estimate. Similarly you could construct a KDE with such negative weights, provided you can accept the possibility that some of its values will be negative! (cc @Xi'an)
– whuber
Feb 1 '18 at 21:42
• What would an "exact" sample of a distribution be? Again, whether and how you can exploit a mixture with negative weights comes down to how you intend to use the sample.
– whuber
Feb 1 '18 at 21:49
• @whuber 'Exact sample' in the sense of equality in distribution - I want a procedure which outputs a pair $(x,y)$, such that the law of the pair is given by the density $p(x,y)$. The sample is then to be fed into another algorithm (to be precise, in my problem, I would use $(x,y)$ as a momentum term in an HMC algorithm), and the only requirement I have is that my pair have the correct distribution. I acknowledge that there are other solutions to my problem (based on e.g. rejection sampling); I'm not currently interested in those.
– πr8
Feb 1 '18 at 21:59

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...

• Thank you! I appreciate the additional references also.
– πr8
Feb 1 '18 at 22:33
• Additional thanks for the very thorough response and references. I am happy to accept this answer as it succeeds in generating exact samples from $p$ in a.s. finite time. I'll likely continue to think about the problem to some extent though; the only additional idea I've had which seems promising is to view sampling from $p = \lambda g - \mu h$ as sampling $X \sim g$, conditional on $\lambda g \geqslant \mu h$, and that there may be some geometric insight which is useful for this characterisation (I'm thinking like a slice sampler on $\{(x,y): \mu h (x) < y < \lambda g(x) \}$). Cheers!
– πr8
Feb 3 '18 at 17:49
• I explained the conditional sampler quite poorly; the set-based characterisation is a bit clearer (in my opinion). My key point is that if you can sample $(x,y)$ uniformly from the two-dimensional set in the final line, it follows that the $x$-coordinate has the correct distribution. Whether this characterisation can be useful for longer sum-based improper mixtures remains to be seen.
– πr8
Feb 3 '18 at 18:03
• I was also thinking of a slice sampler, but this is not "exact" in a simulation sense. Feb 3 '18 at 18:23

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.

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$.

• I considered this but rejected it because my initial efforts to demonstrate it could work led to the realization that it will, at best, be an approximation and potentially a poor one. Yes, asymptotically it could work, but it won't meet the OP's request for "exact" sampling from the distribution.
– whuber
Feb 2 '18 at 18:54
• The efficiency of this method is exactly of the same order as the exact accept-reject method. Feb 3 '18 at 8:56
• Agreed. Yet they are quite different. The accept-reject method needs to compute $g$ and $h$ as functions of $x$. I focused on using only sampling from $g$ and $h$ as "oracle" samplers like in a true mixture. The more I think of it, the more I'm convinced an exact method based on sampling oracles cannot exist. Feb 3 '18 at 10:11
• I think that's generally correct, but there may be useful classes of special cases where such an exact method does exist. That's because (1) in some cases the computation of $g/(g+h)$ is easy and (2) you don't need to compute both $g$ and $h$--you only need to compute this ratio.
– whuber
Feb 3 '18 at 13:27
• @BenoitSanchez Thank you for your in-depth answer; I particularly appreciate the comments at the end about (potential) impossibility of exactness. I've come across Bernoulli Factories in the past and found them quite challenging; I'll try to revisit the topic and see if it provides any insights.
– πr8
Feb 3 '18 at 17:57