# sampling from a mixture of two Gamma distributions

Assuming that all the mixture parameters are known, how can one sample from a mixture of $\texttt{Gamma}(\alpha,\beta)$ distributions: $$\theta \sim \pi \texttt{Gamma}(\alpha_1,\beta_1)+(1-\pi)\texttt{Gamma}(\alpha_2,\beta_2)$$ Is it done by sampling a whole bunch of $\theta$s from either of the distributions and putting them in a pool and at the end computing their corresponding probability according to the mixture, and then sampling proportional to the accumulated computed probabilities?

• – user230053 Oct 12 '19 at 19:48
• How can I apply the EM algorithm to the simulated data (a mixture of two Gamma distributions) in order to estimate the final parameters like p, beta, etc.? alphas are assumed to be fixed and known. – Mostafa Mar 16 at 22:08

Actually, it is much simpler. Assuming that you have $K$-component mixture of gamma distributions the algorithm to draw sample of size $N$ is as follows:

Repeat $N$ times:
1. draw $k$ from categorical distribution parametrized by vector $\pi$,
2. draw single value from $k$-th gamma distribution parametrized by $\alpha_k$, $\beta_k$.

You can use similar algorithm to draw samples from mixture of any distributions. Notice that this follows exactly from the definition of mixture distribution:

mixture distribution is the probability distribution of a random variable that is derived from a collection of other random variables as follows: first, a random variable is selected by chance from the collection according to given probabilities of selection, and then the value of the selected random variable is realized.

If you know R, this translates to the following example:

# density
dmixgamma <- function(x, pi, alpha, beta) {
k <- length(pi)
n <- length(x)
rowSums(vapply(1:k, function(i) pi[i] * dgamma(x, alpha[i], beta[i]), numeric(n)))
}

# random generation
rmixgamma <- function(n, pi, alpha, beta) {
k <- sample.int(length(pi), n, replace = TRUE, prob = pi)
rgamma(n, alpha[k], beta[k])
}

set.seed(123)

pi <- c(4/10, 6/10)
alpha <- c(20, 15)
beta <- c(10, 25)

hist(rmixgamma(1e5, pi, alpha, beta), 100, freq = FALSE)
xx <- seq(0, 10, by = 0.001)
lines(xx, dmixgamma(xx, pi, alpha, beta), col = "red")


• Ahhh why didn't I thought about it this way. Sweet! – user3639557 Aug 2 '16 at 10:07
• A secondary question: is there stat package for c++ that you would recommend? – user3639557 Aug 2 '16 at 10:08
• @user3639557 sorry but I'm using C++ only accompanied with R via Rcpp, but this is probably not what you are asking ;) – Tim Aug 2 '16 at 10:11
• @user3639557 but if you wanted Rcpp example of sampling from mixture, then here you are: github.com/twolodzko/extraDistr/blob/master/src/… (if you accidentally spot any bugs, posting comments and rising issues are welcomed) – Tim Aug 2 '16 at 10:15