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