library(purrr)
library(stringr)

1.

a.

We already observed all \(k = 20\) successes in \(n = 20\) jobs in a day. Under \[ H_0: p = p_0 = 0.92 , \] the p-value is given by the probability of seeing an outcome at least as extreme as 20 successes in 20 jobs. Note that \[ H_A: p > p_0 = 0.92 , \] to have an desired outcome, we want \[\begin{align*} pv = & P (X \geq k) \\ = & 1 - P (X < k) \\ = & 1 - P (X \leq k - 1) \\ = & \sum_{i = k} ^ {n} \binom{n}{i} p_0 ^ {i} (1 - p_0) ^ {n - i} \\ = & \sum_{i = 20} ^ {20} \binom{20}{i} 0.92 ^ {i} (1 - 0.92) ^ {20 - i} \\ = & 0.92 ^ {20} \\ = & 0.1886933 . \end{align*}\]

n <- 20
k <- 20
p0 <- 0.92

### p-value. Be careful about (k - 1)
pbinom(q = k - 1, size = n, prob = p0, lower.tail = FALSE)
[1] 0.1886933
### p-value. Be careful about (k - 1)
1 - pbinom(q = k - 1, size = n, prob = p0)
[1] 0.1886933
### Special case here
dbinom(k, size = n, prob = p0)
[1] 0.1886933

b.

### Discretize the candidate p0
p0Vec <- seq(0, 1, length.out = 10001)

alp <- 0.05

### Method 1: pbinom supports vector input
pvVec <- 1 - pbinom(q = k - 1, size = n, prob = p0Vec)

### Interval esimate
c(Lower = p0Vec[min(which(pvVec > alp))], Upper = p0Vec[max(which(pvVec > alp))])
 Lower  Upper 
0.8609 1.0000 
### Method 2:

### Compute p-value across p0Vec
pvVec <- map_dbl(p0Vec, function(p0){
  1 - pbinom(q = k - 1, size = n, prob = p0)
})

### Interval esimate
c(Lower = p0Vec[min(which(pvVec > alp))], Upper = p0Vec[max(which(pvVec > alp))])
 Lower  Upper 
0.8609 1.0000

c.

### Remember to specify the alternative correctly
binom.test(x = k, n = n, p = p0, alternative = "greater")

    Exact binomial test

data:  k and n
number of successes = 20, number of trials = 20, p-value = 0.1887
alternative hypothesis: true probability of success is greater than 0.92
95 percent confidence interval:
 0.8608917 1.0000000
sample estimates:
probability of success 
                     1

The ouput agrees with part a and b.

2.

a.

mu1Vec <- c(5, 5)
mu2Vec <- c(5, 5.5)
mudiff <- 0
n <- 50
sd <- sqrt(1)
alp <- 0.05
nSim <- 1000

### For two scenarios
outQ2a <- map2_dbl(mu1Vec, mu2Vec, function(mu1, mu2) {
  ### Repeat the experiment and output logical vector
  I <- map_lgl(seq_len(nSim), function(iSim) {
    ### Generate data sets
    X1 <- rnorm(n, mean = mu1, sd = sd)
    X2 <- rnorm(n, mean = mu2, sd = sd)

    ### 95% two-sided confidence interval
    CI <- t.test(X1, X2, alternative = "two.sided", mu = mudiff, conf.level = 1 - alp)$conf.int

    ### Check if 0 is not covered in the CI
    ### TRUE if not covered
    (mudiff < CI[1]) | (mudiff > CI[2])
  })

  ### Compute the proportion
  sum(I) / nSim
})

names(outQ2a) <- str_c("Scenario ", 1:2)

outQ2a
Scenario 1 Scenario 2 
     0.042      0.681

b. (Optional)

(Not graded) 

### Method 1:
nby <- 5

### For two scenarios
outQ2b <- map2_dbl(mu1Vec, mu2Vec, function(mu1, mu2) {
  ### Repeat the experiment and output logical vector
  I <- map_lgl(seq_len(nSim), function(iSim) {
    n0 <- 0
    X1 <- NULL
    X2 <- NULL

    while (TRUE) {
      ### Current number in each group
      n0 <- n0 + nby

      ### Check if the number exceed the limit
      if (n0 > n) {
        break
      }

      ### Combine data sets
      X1 <- c(X1, rnorm(nby, mean = mu1, sd = sd))
      X2 <- c(X2, rnorm(nby, mean = mu2, sd = sd))

      ### 95% two-sided confidence interval
      CI <- t.test(X1, X2, alternative = "two.sided", mu = mudiff, conf.level = 1 - alp)$conf.int

      ### Check if 0 is not covered in the CI
      ### TRUE if not covered
      Ireject <- (mudiff < CI[1]) | (mudiff > CI[2])
      if (Ireject) {
        break
      }
    }

    Ireject
  })

  ### Compute the proportion
  sum(I) / nSim
})


names(outQ2b) <- str_c("Scenario ", 1:2)

outQ2b
Scenario 1 Scenario 2 
     0.195      0.804 
### For two scenarios
outQ2b <- map2_dbl(mu1Vec, mu2Vec, function(mu1, mu2) {
  ### Repeat the experiment and output logical vector
  I <- map_lgl(seq_len(nSim), function(iSim) {
    ### Generate data sets
    X1 <- rnorm(n, mean = mu1, sd = sd)
    X2 <- rnorm(n, mean = mu2, sd = sd)

    ### Index for subsequence
    Isub <- map(seq(5, n, by = nby), seq_len)

    ### Indicator for reject the null hypothesis
    Ireject <- map_lgl(Isub, function(i) {
      ### 95% two-sided confidence interval for all subsequences
      CI <- t.test(X1[i], X2[i], alternative = "two.sided", mu = mudiff, conf.level = 1 - alp)$conf.int

      ### Check if 0 is not covered in the CI
      ### TRUE if not covered
      (mudiff < CI[1]) | (mudiff > CI[2])
    })
    ### Check if null is reject
    sum(Ireject) != 0
  })

  ### Compute the proportion
  sum(I) / nSim
})

names(outQ2b) <- str_c("Scenario ", 1:2)

outQ2b
Scenario 1 Scenario 2 
     0.212      0.792 
### Method 2:
nby <- 5

### For two scenarios
outQ2b <- map2_dbl(mu1Vec, mu2Vec, function(mu1, mu2) {
  ### Repeat the experiment and output logical vector
  I <- map_lgl(seq_len(nSim), function(iSim) {
    ### Generate data sets
    X1 <- rnorm(n, mean = mu1, sd = sd)
    X2 <- rnorm(n, mean = mu2, sd = sd)

    ### Index for subsequence
    Isub <- map(seq(5, n, by = nby), seq_len)

    ### Indicator for reject the null hypothesis
    Ireject <- map_lgl(Isub, function(i) {
      ### 95% two-sided confidence interval for all subsequences
      CI <- t.test(X1[i], X2[i], alternative = "two.sided", mu = mudiff, conf.level = 1 - alp)$conf.int

      ### Check if 0 is not covered in the CI
      ### TRUE if not covered
      (mudiff < CI[1]) | (mudiff > CI[2])
    })
    ### Check if null is reject
    sum(Ireject) != 0
  })

  ### Compute the proportion
  sum(I) / nSim
})

names(outQ2b) <- str_c("Scenario ", 1:2)

outQ2b
Scenario 1 Scenario 2 
     0.176      0.804

c. (Optional)

(Not graded)

The type I error is given as \[ P (\text{Type I error}) = P (\text{Reject } H_0 \mid H_0 \text{ is TRUE}) , \] which is the proportion we obtained under scenario 1. That means we have type I error 0.042 for standard experiment and 0.176 for Scenario 2. For sequential experiments, we have more chances to reject \(H_0\) compared with standard experiments, as we test for every \(5\) subjects. So we are more likely to conduct type I error, which means the test is over-liberal. Hence, to obtain the desired nominal level, we should better use standard experiment.

3.

nSim <- 100
alp <- 0.05
nVec <- c(10, 50, 100, 200)
pVec <- c(0.01, 0.1, 0.5)

### Experiment with different n
outQ3 <- map(nVec, function(n) {
  ### Experiment with different p
  map(pVec, function(p){
    ### X ~ Bin (n, p)
    X <- rbinom(n = nSim, size = n, prob = p)
    ### T
    T <- sqrt(n) * ((X / n) - p)
    
    ### KS test
    ### Ignore the warnings here
    ks <- ks.test(T, pnorm, mean = 0, sd = sqrt(p * (1 - p)))
    ### Compare p-value with alpha
    ### If TRUE, reject the null hypothesis that two distributions equal
    ks$p.value <= alp
  })
})

### Combine as a matrix
outQ3 <- do.call(rbind, outQ3)
dimnames(outQ3) <- list(str_c("size = ", nVec), str_c("prob = ", pVec))
outQ3
           prob = 0.01 prob = 0.1 prob = 0.5
size = 10  TRUE        TRUE       TRUE      
size = 50  TRUE        TRUE       FALSE     
size = 100 TRUE        TRUE       FALSE     
size = 200 TRUE        FALSE      TRUE

The conclusion would be made accordingly, that is, if we have TRUE, we have sufficient evidence to reject the null and conclude that two distributions are different, otherwise, we do not have sufficient evidence to reject the null and conclude that it is plausible that two distributions are the same.

4. (Optional)

(Not graded) 

nSim <- 100
n <- 8000
side <- 10000
alp <- 0.05

### Repeat the experiment
EVec <- map_dbl(seq_len(nSim), function(iSim) {
  ### Generate data from the square
  X <- matrix(runif(2 * n, min = 0, max = side), ncol = 2)
  colnames(X) <- c("x", "y")

  ### Compute different combinations of distances
  dis <- dist(X)

  ### Compute D, the minimum distance
  D <- min(dis)

  ### Compute E
  E <- 1 - exp(-D^2 / 0.995)

  E
})

### Compare p-value with alpha
### If TRUE, reject the null hypothesis that two distributions equal
ks.test(EVec, punif)$p.value <= alp
[1] FALSE

FALSE here means we do NOT have sufficient evidence to reject the null hypothesis that two distributions equal, which means \(E\) is plausible to be uniform on \([0, 1]\).

5.

(Not graded)

a.

n <- 200
lam1 <- 1/2
lam2 <- 1/30
T <- seq(0, 120, length.out = 10000)

### Simulate the waiting time
tW <- rexp(n, rate = lam1)

### Simulate the time spent
tS <- rexp(n, rate = lam2)

### Compute the arrival time
tA <- cumsum(tW)

### Compute the departure time
tD <- tA + tS

### Over the first 120 minutes
nC <- map_dbl(T, function(t){
  ### Number of customer in the restaurant
  sum((tA < t) & (tD > t))
})

### Plot the number of customers in the restaurant over the first 120 minutes after the restaurant opens
plot(T, nC, type = "l", xlab = "Time (min)", ylab = "Number of customers", main = "The number of customers over the first 120 minutes.")

b.

nSim <- 1000
TQ5bi <- c(10, 20, 30, 120)
nQ5bii <- 10
TQ5bii <- 60
nQ5biii <- 20
TQ5biii <- 60

outQ5b <- map(seq_len(nSim), function(iSim) {
  ### Simulate the waiting time
  tW <- rexp(n, rate = lam1)

  ### Simulate the time spent
  tS <- rexp(n, rate = lam2)

  ### Compute the arrival time
  tA <- cumsum(tW)

  ### Compute the departure time
  tD <- tA + tS

  ### Q5bi
  nQ5bi <- map_dbl(TQ5bi, function(t) {
    ### Number of customer in the restaurant
    sum((tA < t) & (tD > t))
  })

  ### Q5bii
  ### Indicator, TRUE if the condition in the question is satisfied
  IQ5bii <- (sum((tA < TQ5bii) & (tD > TQ5bii))) < nQ5bii

  ### Q5biii
  ### Indicator, TRUE if the condition in the question is satisfied
  IQ5biii <- (sum((tA < TQ5biii) & (tD > TQ5biii))) > nQ5biii

  list(nQ5bi = nQ5bi, IQ5bii = IQ5bii, IQ5biii = IQ5biii)
})

i.

### Compute the average number
colMeans(do.call(rbind, map(outQ5b, `[[`, "nQ5bi")))
[1]  4.286  7.429  9.644 14.848

ii.

### Compute the proportion
sum(map_lgl(outQ5b, `[[`, "IQ5bii")) / nSim
[1] 0.15

iii.

### Compute the proportion
sum(map_lgl(outQ5b, `[[`, "IQ5biii")) / nSim
[1] 0.021