##
## Use simulation to estimate the probability of getting at least nine
## heads in ten tosses of a fair coin.
##


##
## Solution using loops.
##

## Number of simulation replications.
nrep = 1e4

## Number of times at least 9 flips give a head.
n9 = 0

## Simulation loop.
for (r in 1:nrep) {

  ## The results from one set of 10 flips.
  C = runif(10) < 0.5

  ## Check whether we got all heads.
  if (sum(C) >= 9) { n9 = n9 + 1 }
}

## The probability estimate is the proportion of the time we got at least
## nine heads.
p_est1 = n9/nrep



##
## Vectorized solution.
##

## Generate a nrep x 10 matrix of coin toss results.
C = runif(nrep * 10) < 0.5
C = matrix(C, nrow=nrep, ncol=10)

## Calculate the number of heads in each set.
n_head = apply(C, 1, sum)

## The probability estimate is the proportion of the time we got at least
## nine heads.
p_est2 = mean(n_head >= 9)



##
## For comparison, we can calculate the exact answer using the binomial 
## distribution.
##

## pbinom(q, n, p) is the probability of seeing q or fewer heads out
## of n independent attempts when the probability of each head is p.
p_exact1 = 1 - pbinom(8, 10, 0.5)


## dbinom(x, n, p) is the probability of seeing exactly x heads out of
## n independent attempts when the probability of each head is p.
##
## Here we are calculating: P(x >= 9) = P(X=9) + P(X=10)
p_exact2 = dbinom(9, 10, 0.5) + dbinom(10, 10, 0.5)


## lchoose(n, k) is the logarithm of "n choose k"
## Here we calculate the probability without using any built-in
## binomial functions.  Note that the binomial probabilities are 
## calculated on the log scale.
p_exact3 = exp(lchoose(10, 9) - 10*log(2)) + exp(lchoose(10, 10) - 10*log(2))