# Please load the packages. 
library(ggplot2)
library(ipsRdbs)
# Start with a clean workspace 
rm(list=ls())

1 Modelling of computer failure data

summary(cfail)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    1.00    3.00    3.75    5.00   17.00
 # 95% Confidence interval 
 c(3.75-1.96 * 3.381/sqrt(104), 3.75+1.96*3.381/sqrt(104)) # =(3.10,4.40).
## [1] 3.100193 4.399807
 x <- cfail 
 n <- length(x)
 h <- qnorm(0.975) 
 # 95% Confidence interval Using quadratic inversion 
 mean(x) + (h*h)/(2*n) + c(-1, 1) * h/sqrt(n) * sqrt(h*h/(4*n) + mean(x))
## [1] 3.395836 4.141101
 # Modelling 
 # Observed frequencies 
 obs_freq <- as.vector(table(x))
 # Obtain unique x values 
 xuniques <- sort(unique(x))
 lam_hat <- mean(x)
 fit_freq <- n * dpois(xuniques, lambda=lam_hat)
 fit_freq <- round(fit_freq, 1)
 # Create a data frame for plotting 
 a <- data.frame(xuniques=xuniques, obs_freq = obs_freq, fit_freq=fit_freq)
 barplot(rbind(obs_freq, fit_freq), args.legend = list(x = "topright"), 
 xlab="No of weekly computer failures",  
 names.arg = xuniques,  beside=TRUE, col=c("darkblue","red"), 
 legend =c("Observed", "Fitted"), 
 main="Observed and Poisson distribution fitted frequencies 
 for the computer failure data: cfail")

2 Modelling the number of bombhits data

summary(bombhits)
##    numberhit         freq      
##  Min.   :0.00   Min.   :  1.0  
##  1st Qu.:1.25   1st Qu.: 14.0  
##  Median :2.50   Median : 64.0  
##  Mean   :2.50   Mean   : 96.0  
##  3rd Qu.:3.75   3rd Qu.:181.5  
##  Max.   :5.00   Max.   :229.0
 # Create a vector of data 
 x <- c(rep(0, 229), rep(1, 211), rep(2, 93), rep(3, 35), rep(4, 7), 5)
 y <- c(229, 211, 93, 35, 7, 1) # Frequencies 
 rel_freq <- y/576
 xbar <- mean(x)
 pois_prob <- dpois(x=0:5, lambda=xbar)
 fit_freq <- pois_prob * 576
  #Check 
  cbind(x=0:5, obs_freq=y, rel_freq =round(rel_freq, 4),  
  Poisson_prob=round(pois_prob, 4), fit_freq=round(fit_freq, 1))
##      x obs_freq rel_freq Poisson_prob fit_freq
## [1,] 0      229   0.3976       0.3950    227.5
## [2,] 1      211   0.3663       0.3669    211.3
## [3,] 2       93   0.1615       0.1704     98.1
## [4,] 3       35   0.0608       0.0528     30.4
## [5,] 4        7   0.0122       0.0122      7.1
## [6,] 5        1   0.0017       0.0023      1.3
  obs_freq <- y
  xuniques <- 0:5
  a <- data.frame(xuniques=0:5, obs_freq =y, fit_freq=fit_freq)
  barplot(rbind(obs_freq, fit_freq), 
  args.legend = list(x = "topright"), 
  xlab="No of bomb hits",  
  names.arg = xuniques,  beside=TRUE, 
  col=c("darkblue","red"), 
  legend =c("Observed", "Fitted"), 
  main="Observed and Poisson distribution fitted frequencies 
  for the number of bomb hits in London")