example: 10 random coin flips (1 = heads, 0 = tails) from Bern(rho)

n <- 10
rho <- 0.5
# note: Bernoulli is a special case of Binomial: Bern(rho)=Binom(1,rho) 
# generate sample of n random trials of Binom(1,rho) using
x_sample <- rbinom(n, size=1, prob=rho)  
heads <- sum(x_sample)

x_sample =1, 1, 1, 1, 1, 0, 0, 0, 1, 1 and heads = 7

Create approximate sampling distribution for number of heads from N trials:

N <- 10**5  # number of samples 
heads <- numeric(N) # create a vector to hold number of heads for each sample
# loop to create sample, count number of heads, store in array
for (i in seq(1,N)) {   
  x_sample <- rbinom(n, size=1, prob=rho)
  heads[i] <- sum(x_sample)
}   
head(heads)
## [1] 5 2 5 4 5 4

Plot a histogram of the approximate sampling distribution

xplot = seq(0,n) # possible outcomes for number of heads
# for histogram specify bins centered on possible # of heads
bin_edges <- c(xplot-0.5,n+0.5)
# use freq=TRUE to return the counts in each bin (here saved as histout)
histout <- hist(heads,breaks=bin_edges,freq=TRUE,main="")

head_counts <- histout$counts

head_counts = 91, 918, 4290, 11902, 20540, 24751, 20548, 11497, 4440, 930, 93

Find approximate pmf from the head counts (p_hat, blue o) for each sample
and compare to the exact result (ptheory, red +)

p_hat = head_counts/N   # approximate pmf for simulated sampling dist
plot(xplot,p_hat,type="p",col="blue",ylim=c(0,0.3))
ptheory = dbinom(xplot,n,rho)
lines(xplot,ptheory,type="p",col="red",pch=3)

comparison of mean and std error

# for theory (see notes)
mean_theory <- n*rho 
stderr_theory <- sqrt(n*rho*(1-rho))
# for simulation
mean_approx <- mean(heads)
stderr_approx <- sd(heads)
type mean std dev
theory 5 1.5811388
approx 4.99706 1.5729705

Note: as N increases, approximation approaches theoretical (exact) results.