# topic 11 h
# Load the functions we will need in this script
source( "../assess_normality.R")
source( "../pop_sd.R")
#
# Let us look at the binomial probabilities.
# We did this earlier in Topic 11b but there
# we used the function pbinom. So, if we had
# an problem where the probability of a single
# success is 0.68, then we could ask for
# the probability on 36 trials of getting
# 23 or fewer successes.
pbinom( 23, 36, 0.68)
#
# Now, we will try this as an experiment
# and we will do the experiment 10,000 times.
L1 <- 1:10000
for (i in 1:10000) {
# do a 36 trial experiment
how_many <- 0
for (j in 1:36 ) {
this_single <- sample(1:2,1,prob=c(0.68,0.32))
if( this_single == 1 ){
how_many <- how_many+1}
}
# now record the result as the number of successes
L1[i] <- how_many
}
# then let us look at this population of
# 36 trial experiments
summary( L1 )
boxplot( L1, horizontal=TRUE )
hist( L1, breaks=seq(10,36,1), xaxp=c(10,36,13) )
assess_normality( L1 )
#
# so this looks like a normal distribution
mean(L1 )
pop_sd(L1 )
# let us try to use this to solve our
# earlier problem. We know that the mean
# of a binomial is n*p and the standard
# deviation of the number of successes is
# sqrt(n*p*(1-p))
pnorm( 23, mean=36*0.68,
sd=sqrt(36*0.68*(1-0.68)))
# compare this to the earlier result
pbinom( 23, 36, 0.68)
# We can make an adjustment here because we
# really want to use 23.5
pnorm( 23.5, mean=36*0.68,
sd=sqrt(36*0.68*(1-0.68)))
# But it seems silly to use this normal
# approximation when we have pbinom() which
# just gives the computed value. The normal
# approximation was useful years ago, before
# we had calculators and computers.
#
# The steps in doing a normal approximation
# for the binomial distribution are always the
# same. So we might as well capture those
# steps in a function, pnormbino()
pnormbino <- function( num_success, p, n, lower.tail=TRUE)
{
if(p*n < 10) {return("n*p < 10, will not compute this")}
if(n*(1-p) < 10)
{return("n*(1-p)<10, will not compute this")}
nbsd <- sqrt( n*p*(1-p))
prob <- pnorm(num_success, n*p, nbsd, lower.tail)
return( prob )
}
pnormbino( 23.5,0.68,36 )
# we should probably note that pnormbino is NOT
# in our list of functions in our parent
# directory so source will not work for it.
# However, we do have pbinom() so pnormbino()
# is not really needed.