rm(list=ls())

# # # # # # # # # # # # # # # # # # # # # # # # # # # 
# # # Consequences of rounding # # # # # # # # # # 
# # # # # # # # # # # # # # # # # # # # # # # # # # # 


#      Proschan (1963) reported on the times between failures of the
#      air-conditioning equipment in 10 Boeing 720 aircraft. The
#      aircondit data frame contains the intervals for the ninth
#      aircraft while aircondit7 contains those for the seventh
#      aircraft.
# 
#      Both data frames have just one column.  Note that the data have
#      been sorted into increasing order.
# 
#   Both datasets are also available in the library "boot"

aircond <- c(3, 5, 7, 18, 43, 85, 91, 98, 100, 130, 230, 487);
aircond7 <- c(3, 5, 5, 13, 14, 15, 22, 22, 23, 30, 36, 39, 44, 46, 50, 72, 79, 88, 97, 102, 139, 188, 197, 210);

values <- unique(aircond7);
(ni <- table(aircond7));
n <- sum(ni);

hist(aircond7)

a <- values - 0.5;
b <- values + 0.5;

# Initial estimator 
lambdahat0 <- 1/mean(aircond7);

lambdahat <- lambdahat0;

n.iter <- 9;

for(i in 1:n.iter){
# #   E-step 
 Xk <- (a*exp(-a*lambdahat) - b*exp(-b*lambdahat))/(exp(-a*lambdahat) - exp(-b*lambdahat)) + 1/lambdahat;
#  M-step 
 lambdahat <- n/sum(ni*Xk);
 print(c(i, lambdahat), digits=15)
}

# Compare with the estimate of lambda ignoring rounding 
lambdahat0;


# # # 
# Suppose now a more severe interval censoring (intervals are 0-10, 10-20, 20-30, ...)
(aircond7.new <- floor(aircond7/10)*10 + 5)
values <- unique(aircond7.new);
(ni <- table(aircond7.new))
n <- sum(ni);

a <- values - 5;
b <- values + 5;
# unique(c(a,b))


# Initial estimator 
lambdahat0 <- 1/mean(aircond7.new);

lambdahat <- lambdahat0;

n.iter <- 9;

for(i in 1:n.iter){
# #   E-step 
 Xk <- (a*exp(-a*lambdahat) - b*exp(-b*lambdahat))/(exp(-a*lambdahat) - exp(-b*lambdahat)) + 1/lambdahat;
#  M-step 
 lambdahat <- n/sum(ni*Xk);
 print(c(i, lambdahat), digits=15)
}

# Compare with the estimate of lambda ignoring rounding 
lambdahat0;


# # # 
# Suppose now even a more severe interval censoring (intervals are 0-50, 50-100, 100-150, ...)
(aircond7.new <- floor(aircond7/50)*50 + 25)
values <- unique(aircond7.new);
(ni <- table(aircond7.new))
n <- sum(ni);

a <- values - 25;
b <- values + 25;
# unique(c(a,b))

# Initial estimator 
lambdahat0 <- 1/mean(aircond7.new);

lambdahat <- lambdahat0;

n.iter <- 9;

for(i in 1:n.iter){
# #   E-step 
 Xk <- (a*exp(-a*lambdahat) - b*exp(-b*lambdahat))/(exp(-a*lambdahat) - exp(-b*lambdahat)) + 1/lambdahat;
#  M-step 
 lambdahat <- n/sum(ni*Xk);
 print(c(i, lambdahat), digits=15)
}

# Compare with the estimate of lambda ignoring rounding 
lambdahat0;