HW10comp

Ex 21: Please generate random numbers from a Binomial distribution with n=20 (number of trials)and p=0.15 (probability of success), and compare the coverage and length of 95% asymptotic confidence intervals and 95% exact confidence intervals. Number of replications = 1000, and seed numbers from 1 to 1000, respectively.

Newton Raphson function

newtonraphson <- function(ftn, x0, tol = 1e-9, max.iter = 100) {
        x <- x0 # x0: the initial value
        fx <- ftn(x)
        iter <- 0
        while ((abs(fx[1]) > tol) & (iter < max.iter)) {
                x <- x - fx[1]/fx[2]
                fx <- ftn(x)
                iter <- iter + 1
                # cat("At iteration", iter, "value of x is:", x, "\n")
        }
        if (abs(fx[1]) > tol) {
                # cat("Algorithm failed to converge\n")
                return(NULL)
        } else { # abs(fx[1]) <= tol
                # cat("Algorithm converged\n")
                return(x)
        }
}

Basic settings

n <- 20;p<- 0.15;no.rep <- 1000

coverage_asym <- rep(NA,no.rep)
length_asym <- rep(NA,no.rep)
coverage_exact <- rep(NA,no.rep)
length_exact <- rep(NA,no.rep)

for (i in 1:no.rep) {
        set.seed(i)
        binom_data <- rbinom(1, size = n, prob = p)
        phat <- binom_data/n
        asy_l <- phat-(qnorm(0.975)*sqrt((phat*(1-phat)/20)))
        asy_u <- phat+(qnorm(0.975)*sqrt((phat*(1-phat)/20)))
        if (asy_l <= 0) {
                asy_l <- 0
        }
        coverage_asym[i] <- ((asy_l <= p) & (asy_u >= p))
        length_asym[i] <- (asy_u - asy_l)
        
        # ci_asym <- binom.confint(binom_data, n, methods = "asymptotic")
        # coverage_asym[i] <- (ci_asym[1] <= p) & (ci_asym[2] >= p)
        # length_asym[i] <- ci_asym[,6] - ci_asym[,5]
        

        ftn1 <- function(p){
                fp1 <- (-0.975)
                dfp1 <- 0
                for(i in 0:(binom_data-1)){
                        fp1 <- fp1 + (choose(20,i)*(p^i)*((1-p)^(20-i)))
                        dfp1 <- dfp1 + choose(20,i)*(i*(p^(i-1))*((1-p)^(20-i))-(p^i)*(20-i)*((1-p)^(19-i)))
                }
                return(c(fp1,dfp1))
        }
        
        ftn2 <- function(p){
                fp2 <- (-0.025)
                dfp2 <- 0
                for(i in 0:binom_data){
                        fp2 <- fp2 + (choose(20,i)*(p^i)*((1-p)^(20-i)))
                        dfp2 <- dfp2 + choose(20,i)*(i*(p^(i-1))*((1-p)^(20-i))-(p^i)*(20-i)*((1-p)^(19-i)))
                }
                return(c(fp2,dfp2))
        }
        
        
        exact_l <- newtonraphson(ftn1,0.05,tol = 1e-9,1000)
        exact_u <- newtonraphson(ftn2,0.3,tol = 1e-9,1000)
        coverage_exact[i] <- (exact_l <= p) & (exact_u >= p)
        length_exact[i] <- exact_u - exact_l
        # ci_exact <- binom.confint(binom_data, n, methods = "exact")
        # coverage_exact[i] <- (ci_exact[1] <= p) & (ci_exact[2] >= p)
        # length_exact[i] <- ci_exact[,6] - ci_exact[,5]
}

coverage_asym <- mean(coverage_asym)
coverage_exact <- mean(coverage_exact)

length_asym <- mean(length_asym)
length_exact <- mean(length_exact)

cat("coverage of asymptotic vs exact CI:", coverage_asym,"vs",coverage_exact, "\n")

## coverage of asymptotic vs exact CI: 0.803 vs 0.975

cat("CI length of asymptotic vs exact method:", length_asym,"vs",length_exact, "\n")

## CI length of asymptotic vs exact method: 0.2765599 vs 0.3136801