Cinghio
Cinghio

Reputation: 332

How to know from which of 2 distributions the sample is drawn in R?

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.

Upvotes: 0

Views: 47

Answers (2)

Bryan Goggin
Bryan Goggin

Reputation: 2489

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.

Upvotes: 2

DatamineR
DatamineR

Reputation: 9618

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"

Upvotes: 1

Related Questions