#
# BinomialExample1.R
#
# This file presents examples for how to use the binomial probability
# density and random number generator in R
#
# Please see the lecture notes for mathematical details on the
# binomial distribution.


# Recall that the binomial distribution makes use of the factorial
# operator (e.g., 4!).  In R, the n! is generated by gamma(n+1).

#Example of 4!
#by hand:
4*3*2*1

#4! using the R gamma command:
gamma(4+1)

#or more directly:
gamma(5)

#factorials grow very quickly in size:
gamma(10+1)

gamma(20+1)

gamma(40+1)

# n choose s in R is choose(n,s)

#EXAMPLE: Equation 100 in the notes (8 choose 1), first by and:
gamma(8+1)/(gamma(1+1)*gamma(8-1+1))

#or directly using the choose command:
choose(8,1)

#EXAMPLE: Equation 101 in the notes (69 choose 2), first by hand:
gamma(69+1)/(gamma(2+1)*gamma(69-2+1))

#or directly using the choose command:
choose(69,2)

#one can always get help by typing "help(command)":
help(choose)


#EXAMPLE: Equation 102---What's the probability of observing 2
#revolutions out of 69 chances if we use the binomial distribution with
#pi=2/69

#first by hand:
p  <- 2/69
choose(69,2)*p^2*(1-p)^(69-2)

#Equation 102 using the R command for the binomial PDF: dbinom()
#dbinom expects the following three arguments:
#1) "x": number of successes (e.g., revolutions)
#2) "size": the number of trials
#3) "prob": the probability of success in a single trial
# also see help(dbinom)

dbinom(2,69,p)

#EXAMPLE: Equation 103---What's the probability of observing
#9 revolutions out of 69 chances if we use the
#binomial distribution with pi=2/69

#by hand
p  <- 2/69
choose(69,9)*p^9*(1-p)^(69-9)

#using the binom pdf function dbinom()
dbinom(9,69,p)

#EXAMPLE: Equation 104--What's the probability of observing 2
#revolutions out of 69 chances if we use the binomial distribution with
#pi=1/8

#by "hand"
p  <- 1/8
choose(69,2)*p^2*(1-p)^(69-2)

#using the binom pdf function dbinom()
dbinom(2,69,p)

# Let's draw random numbers from the binomial random number function
# rbinom().  The function expects three arguments:
#1) "n": the number of random draws
#2) "size": the number of trials
#3) "prob": the probability of success in a single trial

#EXAMPLE: what proportion of random draws from a binomial distribution
#with a trial size of 69 and pi=2/69 will result in exactly 2
#successes?  This question is related to Equation 102

#let's take 1000 random draws
foo  <- rbinom(1000,69,2/69)

#let's look at the results
print(foo)

#how many of the results are exactly equal to 2 successes:
n2  <- sum(foo==2)
print(n2)

#what proportion of all draws is n2?
n2/1000

#but this differes slightly from Equation 102 which is .2746708
#Let's take 100,000 draws
n   <- 100000
foo <- rbinom(n,69,2/69)
n2  <- sum(foo==2)
n2/n

#now we are getting closer to the number from Equation 102.  We could
#keep increasing the number of random draws until we get exactly
#.2746708.  The difference between the probability given by the PDF
#and the proportion found by random draws shrinks as we make larger
#numbers of random draws.