set.seed(123)
simulate_GW_once <- function(p, n) {
X <- 1 #starting with a single cell
for (gen in 1:n) {
X <- rbinom(1, size = X, prob = p) #each cell either survives or dies
}
return(X)
}
#parameters
p <- 0.7
n <- 5
R <- 10000
results <- replicate(R, simulate_GW_once(p, n))
mean(results)
## [1] 0.1651
table(results) / R #empirical distribution
## results
## 0 1
## 0.8349 0.1651
#theoretical probability of Xn=1
theo <- p^n
theo
## [1] 0.16807
#Both the theoretical probability and emperical probability are almost similar