0

I am writing a pop simulation using R. The "growth" parameter in this simulation is sampled with probability p from two normal distributions with different means and sd. The code I am using to do so is the following:

R_0 <- sample(c(sample(rnorm(5000, mean = 3.5, sd = 0.7), 1),sample(rnorm(5000, mean = 1.5, sd = 0.3), 1)),1, prob = c(1 - p, p))

The two distributions correspond to a benign and stressful environment (B and S for simplicity). My problem is that later in the code I need to modify the growth rate R_0 according to a fixed quantity. This quantity depends on whether the R_0 coefficient was sampled from the B or the S environments. Is there a way to know from which of the two distributions R_0 is sampled?

I guess I could use and if else statement but I am not quite sure how to initialize it.

Thanks in advance.

Cinghio
  • 302
  • 1
  • 4
  • 14

2 Answers2

2

Try this. I have a feeling the output is not exactly in the form you would like, but without more context I cannot know what would work best for your case:

set.seed(9087)

R_0 <- function(p){
  pick<- sample(c("B","S"), 1, prob = c(p, 1-p))
  if(pick == "B") samp<-sample(rnorm(5000, mean = 3.5, sd = 0.7), 1)
  if(pick == "S") samp<-sample(rnorm(5000, mean = 1.5, sd = 0.3), 1)
  return(paste("Value of", samp, "drawn from", pick, sep = " "))
}

R_0(.5)

[1] "Value of 5.31171245776231 drawn from B"

If you cant figure out how to make this output work for your situation, let me know.

Bryan Goggin
  • 2,449
  • 15
  • 17
  • Hi Bryan thanks for the suggestion. I think DatamineR approach (see below) may work better for what I need to do. – Cinghio Jun 25 '16 at 21:44
1

First save the first sampling result:

set.seed(123)
p <- .5
my_sample <- c(sample(rnorm(5000, mean = 3.5, sd = 0.7), 1),sample(rnorm(5000, mean = 1.5, sd = 0.3), 1))
R_0 <- sample(my_sample,1, prob = c(1 - p, p))

Then you can check from which (first/A or second/B distribution the sample came from):

my_sample == R_0
[1] FALSE  TRUE

There is a tiny probability that the two samples values are identical. So you could also use a flag which tellst you from which sample the final sample came from:

p <- .5
set.seed(123)
sample(c(paste("A", sample(rnorm(5000, mean = 3.5, sd = 0.7), 1), sep = " "),
    paste("B", sample(rnorm(5000, mean = 1.5, sd = 0.3), 1), sep = " ")),1, prob = c(1 - p, p)) 
[1] "B 1.45624041440113"
DatamineR
  • 10,428
  • 3
  • 25
  • 45
  • Hi DatamineR, thanks for the suggestion. The code you provided seem to work fine. My simulation runs smoothly now. – Cinghio Jun 25 '16 at 21:45