# 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.