# Drawing samples from a finite mixture of normal distributions?

After some Bayesian update steps, I am left with a posterior distribution of the form of a mixture of normal distributions,$$\Pr(\theta| \text{data} ) = \sum_{i=1}^k w_i N(\mu_i, \sigma^2).$$ That is, the parameter $\theta$ is drawn from a distribution whose PDF is given as a weighted mixture of normal PDFs, and is not a sum of normal RVs. I would like to draw samples $\theta\sim\Pr(\theta|\text{data})$ to use in an importance sampling approximation of this posterior. In practice, the sum over $i$ can have a large number of terms, so that it can be impractical to choose a term $i$ according to the weights $\{w_i\}$ and then draw $\theta\sim N(\mu_i, \sigma^2)$. Is there an efficient way of drawing samples from a posterior of this form?

• Have you actually tried the select then throw method? The selection can be made reasonably fast of O(k) steps go. Jan 27 '12 at 18:25
• If Barron's solution is really not correct, and you in fact mean a "mixture model", could you please use that term? Jan 27 '12 at 22:15
• Neil G: I'm not a statistician by trade, rather a physicist that sometimes needs to make use of statistics. As such, I didn't know the appropriate term to describe what I needed. I can go on and edit the question now, though, to make it more clear that the PDFs are being summed and not the RVs. Jan 27 '12 at 22:46
• @ChrisGranade: I wasn't trying to come down on you. I just wanted to make sure that's what you meant, and to suggest the edit. Jan 28 '12 at 0:57
• Why is it impractical to choose $i$ based on the weights $\{w_i\}$ and a sample from the uniform distribution on $[0,1]$, then sample $N(\mu_i,\sigma^2)$? This is only moderately more expensive than sampling a single normal distribution, the cost is independent of the number of mixed distributions $k$ and does not rely on those distributions being normal. Jan 28 '12 at 16:44

In principle one could preselect the number of samples to be drawn from each sub-distribution, then visit each sub distribution only once and draw than number of points.

That is

1. Find the random set $<n_1, n_2, \dots, n_k>$ such that $n = \sum_{i=1}^k n_i$ and respecting the weights.

I believe that you do this by drawing a Poisson distribution a multinomial distribution (see the comments) of mean $w_i * n$ for each sub-distribution and then normalizing the sum to $n$.

The work here is $\mathcal{O}(k) * \mathcal{O}(n)$

2. Then do

for (i=1; i<=k; ++i)
for (j=1; j<=n[i]; ++j)
theta ~ N(mu[i],sigma[i])


The work here is $\mathcal{O}(n)$

Though this means that you don't get the in random order. If random order is required you must then shuffle the draws (also big $\mathcal{O}(n)$).

It looks like the first step is dominate in run time and of the same order as the naive algorithm, but if you are sure that all $w_i * n \gg 1$ you could approximate the Poisson distributions with Normal distributions and speed up the first step.

• The distribution of the $n_i$ is not a Poisson distribution if $n$ is fixed, but a binomial distribution. Jan 27 '12 at 21:53
• @FrédéricGrosshans Uhm...here's where I admit my distressing weakness in probability. Looking I think you may be right. I don't have a link for throwing arbitrary binomial distributions, but wikipedia has some references. There is also a relationship between Poisson and Binomial which I'm going to claim was responsible for my uncertainty. Yeah, that's the ticket. Jan 27 '12 at 22:00
• @dmckee: Good answer for drawing from a mixture model, except that it should be a multinomial distribution rather than a Poisson distribution in step 1. Jan 27 '12 at 22:26

Note: The original version of this question asked about a "weighted sum of normal distributions" to which the following answer might be useful. However, after a good bit of discussion on this answer, the answer by @Geoff, and on the question itself, it became clear the question was really on sampling a "mixture of normal distributions" to which this answer is not applicable.

The sum of normal distributions is a normal distribution, so you could calculate the parameters of this single distribution and then simply draw samples from that. If we call that distribution $N(\mu_{sum},\sigma_{sum}^2)$ then,

$$\mu_{sum} = \sum_{i=1}^k w_i\mu_i$$

$$\sigma_{sum}^2=\sum_{i=1}^k w_i^2 \sigma_i^2$$

• To put it succinctly, Chris is summing probability density functions, not random variables. Jan 27 '12 at 21:25
• Chris wants a PDF that has (at least in principle) multiple bumps in it. That is, he was the sum of PDFs, not the PDF of a sum. Jan 27 '12 at 22:13
• It is true that the sum of normally distributed random variables is itself a normally distributed random variable. However, the sum of normal distributions is not a normal distribution. So if $X_{1} \sim N(\mu_{1},\sigma_{1}^2)$ and $X_{2} \sim N(\mu_{2}, \sigma_{2}^{2})$, it is true that $X_{1} + X_{2} \sim N(\mu_{1} + \mu_{2}, \sigma_{1}^{2} + \sigma_{2}^{2})$, but $PDF(X_{1} + X_{2}) \neq PDF(X_{1}) + PDF(X_{2})$. (Credit goes to @ChrisGranade for the explanation.) Jan 27 '12 at 22:15
• @dmckee: that's not a "weighted sum of normal distributions", that's a "mixture of normal distributions". Jan 27 '12 at 22:17
• @Barron comments aren't considered an essential part of the page. You should definitely edit your answer to include the gist of the comments so that readers who don't look at the comments don't get misled. Jan 28 '12 at 3:01

Update: This answer is incorrect, stemming from confusion in terminology (see the comment chain below for details); I'm only leaving it up as a guidepost so that people do not repost this answer (besides Barron). Please do not vote it up or down.

I'd just use properties of random variables to reduce it to a single normally distributed random variable. The sum of two independent, normally distributed random variables is itself a random variable, so if $X_{1} \sim N(\mu_{1}, \sigma_{1}^{2})$ and $X_{2} \sim N(\mu_{2}, \sigma_{2}^{2})$, then

$$X_{1} + X_{2} \sim N(\mu_{1} + \mu_{2}, \sigma_{1}^{2} + \sigma_{2}^{2}).$$

Also, if $w_{1} \in \mathbb{R}$, then

$$w_{1}X_{1} \sim N(w_{1}\mu_{1}, w_{1}^{2}\sigma_{1}^{2}).$$

Using these two results combined, then

$$Pr(\theta | \rm{data}) \sim N\big(\sum_{i=1}^{k}w_{i}\mu_{i}, \sum_{i=1}^{k}w_{i}^{2}\sigma_{i}^{2}\big).$$

So in this case, you'll only need to pull samples from a single distribution, which should be much more tractable.

• This is the solution to a different problem which can be seen from the fact that the original distribution is multi-modal and your suggestion is uni-modal. Jan 27 '12 at 18:49
• @ChrisFerrie: I believe you, but based on the notation, I'm confused as to why the distribution above would be multimodal, while the sum of two independent Gaussian random variables would not be. What am I missing here? Jan 27 '12 at 19:10
• I think the confusion is that we aren't looking at a sum of random variables, but a PDF that is the sum of many PDFs. These aren't always the same, since $p(X_1 + X_2)\ne p(X_1) + p(X_2)$. Instead, our PDF can be thought of as marginalizing over the random variable $i$. Jan 27 '12 at 19:21
• Ah, you're looking at sums of PDFs. Yes, that is a different beast entirely. Now that I read the question more closely, I see what you're saying, and I'm going to delete my response. Thanks! Jan 27 '12 at 19:24
• I have undeleted my previously deleted answer only to serve as a guidepost for others so that no one else answers this question like Barron and I did. Please do not up or down vote my answer anymore. Jan 27 '12 at 23:17