library(MASS)

p = 0.4
n = 10000

# Simulate from Bernoulli distribution
u = runif(n)

simBer = as.numeric(u < p)
hist(simBer)
sum(simBer)

# Simulate from Binomial distribution
N = 20
simBin = rbinom(n, size=N, prob=p)
simBin_manual = c()
for(i in 1:n){
  u = runif(N)
  simBin_manual[i] = sum(u < p)
}

par(mfrow=c(1,2))
hist(simBin, ylim=c(0, 0.25), freq=FALSE)
hist(simBin_manual, ylim=c(0, 0.25), freq=FALSE)
    
mean(simBin_manual)    
var(simBin_manual)