Question 4c

Parts (i) and (ii)

M = 1000; n = 20; p = 0.2
set.seed(3291) ## Set your own seed 
n = 20
# Storage for all θ̂) values across simulations
theta_values <- numeric(M)


### (ii) ----------------------------
# Generate one sample of size n from Bernoulli(p_true)
# and compute θ̂ = log(p̂ / (1 - p̂)).
# Repeat this process M = 1000 times.
for (i in 1:M) {
  sample <- rbinom(n, size = 1, prob = p)  # generate Bernoulli sample
  p_hat <- mean(sample)                         # sample proportion (p̂)
  
  # Avoid log(0) or log(∞) by restricting p̂ a little away from 0 and 1
  p_hat <- max(0.001, min(0.999, p_hat))
  
  # Compute θ̂ = log(p̂ / (1 - p̂))
  theta_values[i] <- log(p_hat / (1 - p_hat))
}

Part (iii)

# Plot the sampling distribution of θ̂
hist(theta_values, freq = FALSE,
     main = expression(paste("Sampling distribution of ", hat(theta), ", n=20")),
     xlab = expression(hat(theta)))

# Compute bias and variance
theta_true <- log(p / (1 - p))
bias <- mean(theta_values) - theta_true     # bias = E[θ̂] - θ
var_est <- var(theta_values)                # sample variance of θ̂

# Display results
bias
## [1] -0.1700247
var_est
## [1] 0.7815302

Question 4d

M = 1000; n = 100; p = 0.2
set.seed(3291) ## Set your own seed

# Storage for all θ̂) values across simulations
theta_values <- numeric(M)


### (ii) ----------------------------
# Generate one sample of size n from Bernoulli(p_true)
# and compute θ̂ = log(p̂ / (1 - p̂)).
# Repeat this process M = 1000 times.
for (i in 1:M) {
  sample <- rbinom(n, size = 1, prob = p)  # generate Bernoulli sample
  p_hat <- mean(sample)                         # sample proportion (p̂)
  
  # Avoid log(0) or log(∞) by restricting p̂ a little away from 0 and 1
  p_hat <- max(0.001, min(0.999, p_hat))
  
  # Compute θ̂ = log(p̂ / (1 - p̂))
  theta_values[i] <- log(p_hat / (1 - p_hat))
}
# Plot the sampling distribution of θ̂
hist(theta_values, freq = FALSE,
     main = expression(paste("Sampling distribution of ", hat(theta), ", n=100")),
     xlab = expression(hat(theta)))

# Compute bias and variance
theta_true <- log(p / (1 - p))
bias <- mean(theta_values) - theta_true     # bias = E[θ̂] - θ
var_est <- var(theta_values)                # sample variance of θ̂

# Display results
bias
## [1] -0.02560552
var_est
## [1] 0.05784329

Comment:

Comment (4d):
For a sample size of 100 and a true probability of 0.2, the sampling distribution of the log-odds estimator \(\hat{\theta} = \log(\hat{p}/(1 - \hat{p}))\) is approximately normal and centred close to the true log-odds value of \(-1.386\). The bias is very small, and the variance is much smaller compared to when the sample size was 20. This confirms that \(\hat{\theta}\) is an unbiased and consistent estimator as \(n\) increases, and that the normal approximation becomes highly accurate for inference when the sample size is large.

Question 4e

Part (i)

\[ H_0: \theta \ge 0 \quad \text{versus} \quad H_a: \theta < 0 \] since the log-odds function \(\theta = \log\!\left(\frac{p}{1 - p}\right)\) is a monotonic transformation of \(p\).

Part (ii)

set.seed(3291)

# Define parameters
n <- 200          # Sample size
p <- 0.42    # True probability
M <- 1000         # Number of Monte Carlo simulations
z <- qnorm(0.95)  # 95th percentile of the standard normal distribution (z_0.05)

# Define the log-odds (logit) transformation function
g <- function(p) log(p / (1 - p))

# Create an empty logical vector to store whether we reject H0 in each simulation
reject <- logical(M)

# Run the simulation M times
for (m in 1:M) {
  x  <- rbinom(n, 1, p)      # Generate a sample of 0s and 1s (Bernoulli trials)
  ph <- mean(x)                   # Calculate sample proportion (p-hat)
  ph <- pmin(pmax(ph, 1e-8), 1-1e-8)  # Avoid extreme values (0 or 1) for stability
  th <- g(ph)                     # Compute log-odds (theta-hat)

  # Compute the rejection threshold for a one-sided test
  crit <- -z / sqrt(n * ph * (1 - ph))

  # Check if the estimated log-odds is less than the critical value (reject H0)
  reject[m] <- (th < crit)
}

# Compute the proportion of rejections (empirical significance level or power)
prop_reject <- mean(reject)

# Display the result
prop_reject
## [1] 0.728

Part (iii)

As we simulated under \(p = 0.42 < 0.5\), the value \(0.728\) represents the power of the test.
If the simulation had used \(p = 0.5\), it would estimate the Type I error;
if \(p > 0.5\), it would not represent the power of the test;
and if \(p < 0.5\) (as in this case), the proportion corresponds to the power of the test.