5
$\begingroup$

I want to sample from mixed normal distribution, first one is $N(1,2)$, second one is $N(5,4)$. I used rnorm(100, c(mean=c(1,5), sd=c(2,4))). Is this correct?

The problem I am trying to solve is sampling from the 2 distribution above, first one with 75%, second one with 25%. Am I on the right track?

Edit: I will rewrite the problem for clearance, with easier numbers. :)

I want to sample from $N(0,1)$ with 70% probability, and $N(100,10)$ with 30% probability. Of course, that's just for sake of discussion, the actually distribution I am working with is n(21, 3.3), n(26,4).

$\endgroup$
6
  • 3
    $\begingroup$ how about c(rnorm(75, 1, 5), rnorm(25, 2, 4) $\endgroup$ Dec 20 '11 at 20:26
  • 1
    $\begingroup$ @David - you might want to make an answer out of the comment, so you can get "answer" credit for it. Also, it should be rnorm(75,1,2) and rnorm(25,5,4). $\endgroup$
    – jbowman
    Dec 20 '11 at 20:51
  • $\begingroup$ ^ +1. If you want to permute this result, install gregmisc and use the permute function. $\endgroup$
    – Arun
    Dec 20 '11 at 20:54
  • 1
    $\begingroup$ @David: The 75, and 25 will give me exactly 75 and 25 samples for each distribution. But the 75% and 25% are just the probabilities getting each, not exactly those numbers. $\endgroup$
    – mike
    Dec 20 '11 at 21:15
  • $\begingroup$ @Arun: can I avoid using any package? $\endgroup$
    – mike
    Dec 20 '11 at 21:15
5
$\begingroup$

If you want to sample unequally (with probability 0.7 and 0.3) from two gaussians with parameters $(\mu_1,\sigma_1^2)$ and $(\mu_2,\sigma_2^2)$, then you can probably try something like that:

n <- 100
yn <- rbinom(n, 1, .7)
# draw n units from a mixture of N(0,1) and N(100,3^2)
s <- rnorm(n, 0 + 100*yn, 1 + 2*yn)

In fact, this is one of the illustrations provided in Modern Applied Statistics with S, by Venables and Ripley (Springer, 2002; §5.2, pp. 110-111).

With different parameters, you can use an ifelse expression to select the mean and SD according to the binomial sequence given in yn, e.g. rnorm(n, mean=ifelse(yn, 21, 26), sd=ifelse(yn, 3.3, 4)). (No need to cast yn to a logical with as.logical.)

$\endgroup$
15
  • 1
    $\begingroup$ elegant. impressed! :) $\endgroup$
    – Arun
    Dec 20 '11 at 21:47
  • $\begingroup$ @chl: I am looking at it, I don't see that page. $\endgroup$
    – mike
    Dec 20 '11 at 21:47
  • $\begingroup$ @mike My first example is exactly the same as the code you wrote: as the two theoretical distributions have the same mean but different variance, they only varied the SD, which defaults to 3 when the $i$th Bernoulli event (which has probability 0.05 of occurring) equals 1, 1 otherwise. $\endgroup$
    – chl
    Dec 20 '11 at 21:50
  • $\begingroup$ @chl: I am a little confused. Could you explain what the ifelse() does in here? I am guessing it matches up mu1 with sd1, mu2 with sd2? $\endgroup$
    – mike
    Dec 20 '11 at 21:55
  • 1
    $\begingroup$ @mike If you want to get the variance of the combined sample, you need to account for the different location parameters (or their estimates) as well. $\endgroup$
    – chl
    Dec 21 '11 at 12:07
2
$\begingroup$

To accomplish the goal of sampling from an uneven mixture of distributions, the most straightforward approach is to sample separately, in proportion to the desired ratio:

 p <- 0.70 #P(from N(mu1, sd1)) 
 n.samps <- 10000
 mu1 <- 0
 sd1 <- 1
 mu2 <- 100
 sd2 <- 10

 x <- vector()
 for(i in 1:n.samps){
    b <- runif(1, 0, 1)
    if(b < p){
       x[i] <- rnorm(1, mu1, sd1)
     } else { 
       x[i] <- rnorm(1, mu2, sd2)
     }
   }

this can be done ~50 x faster:

 binary <- runif(n.samps, 0, 1) > p
 x <- c(rnorm(sum(binary), 1, 2), rnorm(sum(!binary), 5, 4)

then to draw a sample:

sample(x, 1)

or to reshuffle:

x <- sample(x, n.samp)
$\endgroup$
3
  • $\begingroup$ I edit the question, hopefully, it's clear now. I have seen packages doing mixture normal, my guess is you can't do it with one line in rnorm()? Or maybe it is very complicated? $\endgroup$
    – mike
    Dec 20 '11 at 21:19
  • 1
    $\begingroup$ A defect of this method is that x is ordered, you would need to shuffle it with sample(x, length(x)). The value of sum(binary) could be generated directly as s <- rbinom(1,n,p) and of course the value of sum(!binary) is n-s. $\endgroup$
    – Elvis
    Dec 21 '11 at 8:09
  • $\begingroup$ @elvis thanks for pointing that out; the for loop output is shuffled, I just tried to vectorize after seeing chi's answer. $\endgroup$ Dec 21 '11 at 16:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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