Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I know there is an easy way to sample standard mixture models if the mixture components can be easily sampled. But, I have a slightly different scenario and I'm wondering if there is an analogous method (or any method at all that is nice). I'll first explain the standard scenario so there is no confusion, and then explain how my problem is different. Sorry if some of this is too remedial, I just want my question to be completely clear.

A standard mixture model has pdf

$$f_X(x) = \sum_i w_i f^i_X(x)$$

where the $w_i$'s are non-negative and sum to one and the $f^i$'s are pdfs. A pretty important mixture model is mixture of Gaussians, where the $f^i$'s are normal distributions with different means and variances, i.e.,

$$f_X(x) = \sum_i w_i \mathcal{N}(\mu_i, \sigma^2_i; x)$$

I can easily sample this distribution with two samples. I'll assume there is a latent variable that governs which of the components generates the value of the random variable. Thus, I can choose $i$ with probability $w_i$ and then sample $\mathcal{N}(\mu_i, \sigma^2_i; \cdot)$ to generate my sample $\bar{x}$.

In my case, I have a slightly different mixture model

$$f_X(x) = \sum_i w_i f^i_X(x) - \sum_j w_j f^j_X(x)$$

Basically, some of the weights are negative so the components $w_j f^j$ subtract from the density at any $x$. Assume that $f_X$ is a pdf, so it is non-negative and integrates to one. Of course, there are methods for generating samples from arbitrary distributions, but it seems that computing the cdf for this one may be hard.

What is the best way to draw a sequence of samples from this distribution?

Thanks in advance for any help you can provide.

share|cite|improve this question
In the standard model, $\sum_i w_i = 1$, and $\int f_i(x) dx =1 \forall i$. What is the corresponding normalization in your model? – highBandWidth Feb 20 '11 at 5:28
$\int f_i(x) dx = 1 ~\forall i$, but weights don't necessarily sum to one. Basically the weights are what they need to be so that f(x) is a pdf with respect to the $f_i$'s and the signs of the weighting coefficients. – RandomGuy Feb 20 '11 at 7:32
Integrating both sides, I think you need $ \sum_i w_i - \sum_j w_j = 1$. But this doesn't guarantee that total $f_X(x)>0 \forall x$. We'll just assumethat this is true. – highBandWidth Feb 20 '11 at 11:58
Yeah, you're right. I actually have this constraint written down in my notes. No useful constraints on the weights guarantee $f_X(x) > 0$ for all $x$, but its not too troublesome due to the way the distributions arise (which is not really relevant to the sampling). – RandomGuy Feb 20 '11 at 18:33
up vote 2 down vote accepted

Since $$\int f_X^i(x) dx = 1 \forall i$$ and $$\int f_X^j(x) dx = 1 \forall j,$$ it is also required that $$\sum_i w_i - \sum_j w_j = 1.$$ Let us call $\sum_i w_i = a$, so we know that $\sum w_j = a-1$.

Here is the algorithm:

  1. Pick $i$ based on weights $w_i$.
  2. Sample y from $f_X^i(x)$.
  3. Calculate $$p_a = \frac{\sum_i w_i f^i_X(x)- \sum_j w_j f^j_X(x)}{\sum_i w_i f_X^i(x)}$$
  4. Sample $u$ from $Uniform(0,1)$.
  5. If $u < p_a$, then emit y, else reject this sample and go back to 1 to try another.

This is called acceptance-rejection sampling. For more information see the wikipedia article.

share|cite|improve this answer
This doesn't work. Try a concrete example (say, sample $2f-g$ where $f$ is $\pm 1$ with equal probability and $g$ is always $1$). – mjqxxxx Feb 20 '11 at 13:34
I'm not sure I understand the intuition. Perhaps you could show how the pdf of the samples generated is (or approximates) $f_X$. – RandomGuy Feb 20 '11 at 18:36
@mjqxxxx, We need to have $2f-g \ge 0 $. Your example doesn't satisfy that. – highBandWidth Feb 21 '11 at 2:12
@RandomGuy, I have corrected my answer, hopefully the rejection sampling article should be clear. – highBandWidth Feb 21 '11 at 18:56
Very cool, seems to work perfectly in the case where the $f^i$'s and $f^j$'s are Gaussians. The paper and book linked on the Wikipedia seemed to be a bit difficult to download, do you have a readily-accessible reference that discusses the random variable transform? – RandomGuy Feb 22 '11 at 0:43

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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