# Additional calculations not in Stoplight.R stoplight <- read.csv(file = "C:\\data\\stoplight.csv") head(stoplight) alpha <- 0.05 n <- length(stoplight$vehicles) mu.hat <- mean(stoplight$vehicles) #Wald using a log() transformation – not in book or program exp(log(mu.hat) + qnorm(p = c(alpha/2, 1-alpha/2)) * sqrt(1/(mu.hat*n))) ######################################################################## #Probability of 9 or more prob9<-1 - ppois(q = 8, lambda = mu.hat) prob9 #Probability of happening at least once in an hour. 1 - dbinom(x = 0, size = 60, prob = prob9) library(car) #Use the g below instead factorial(1:8) #factorial() will not work with deltaMethod() and g so need to put actual numbers in g <- "1 - exp(-60*mu) * (1 + mu + mu^2/2 + mu^3/6 + mu^4/24 + mu^5/120 + mu^6/720 + mu^7/5040 + mu^8/40320)^60" names(mu.hat) <- "mu" #deltaMethod() requires it to have a name save.it <- deltaMethod(object = mu.hat, g = g, vcov. = mu.hat/n) #Don't need parameterNames = "mu" save.it$SE save.it$Estimate #Same as prob9 above save.it$Estimate + qnorm(p = c(0.025, 0.975))*save.it$SE #Wide #(0.3304907, 0.9915973) #Alternative - could specify parameterNames argument again deltaMethod(object = mu.hat, g = g, parameterNames = "mu", vcov. = mu.hat/n) #