用r语言编译： 假设样本X1, ..., Xn 以及Y1, ...Yl 分别是来自总体 X~N(µ1, σ12), Y~N(µ2, σ22) 的简单随机抽样。考虑以下假设检验问题（显著性水平α = 0.05） H0 : µ1 = µ2, σ12=σ22；H1 : µ1 ≠ µ2 或者 σ12 ≠ σ22 假设样本量n = 10, l = 15，完成以下问题： (1) 生成一组X1, ..., Xn,Y1, ...Yl 的数据 (2) 若不考虑假设检验问题，求参数µ1 , µ2, σ12 , σ22的最大似然估计 (3) 在原假设下,设µ =µ1=µ2,σ2 =σ12 =σ22,求µ, σ2的最大似然估计 (4) 写出该检验的对数似然比统计量LRT (5) 使用蒙特卡罗随机模拟重复(1)~(4)，重复次数为m = 10000，画出对数似然比统计量LRT的抽样分布 (6) 若样本量n = 100, l = 150，重复（5）中的步骤，此时对数似然比统计量LRT的抽样分布是什么？和(5)中的抽样分布有什么不同？

首先，我们可以使用R语言代码生成X1, ..., Xn, Y1, ..., Yl的数据：

set.seed(123)
n <- 10
l <- 15
mu1 <- 1
mu2 <- 2
sigma1 <- 1
sigma2 <- 2
x <- rnorm(n, mu1, sigma1)
y <- rnorm(l, mu2, sigma2)

接下来，我们可以求参数µ1, µ2, σ12, σ22的最大似然估计：

log_likelihood <- function(params, x, y) {
  mu1 <- params[1]
  mu2 <- params[2]
  sigma1_sq <- params[3]
  sigma2_sq <- params[4]
  n <- length(x)
  l <- length(y)
  -n/2*log(2*pi*sigma1_sq) - sum((x - mu1)^2)/(2*sigma1_sq) - l/2*log(2*pi*sigma2_sq) - sum((y - mu2)^2)/(2*sigma2_sq)
}

optim_res <- optim(c(mu1=0, mu2=0, sigma1_sq=1, sigma2_sq=1), log_likelihood, x=x, y=y, method="L-BFGS-B", lower=c(-Inf, -Inf, 0, 0))
mu1_hat <- optim_res$par[1]
mu2_hat <- optim_res$par[2]
sigma1_sq_hat <- optim_res$par[3]
sigma2_sq_hat <- optim_res$par[4]

在原假设下，我们可以使用µ = µ1 = µ2, σ^2 = σ1^2 = σ2^2来求µ, σ^2的最大似然估计：

log_likelihood_reduced <- function(params, x, y) {
  mu <- params[1]
  sigma_sq <- params[2]
  n <- length(x)
  l <- length(y)
  -n/2*log(2*pi*sigma_sq) - sum((x - mu)^2)/(2*sigma_sq) - l/2*log(2*pi*sigma_sq) - sum((y - mu)^2)/(2*sigma_sq)
}

optim_res_reduced <- optim(c(mu=0, sigma_sq=1), log_likelihood_reduced, x=x, y=y, method="L-BFGS-B", lower=c(-Inf, 0))
mu_hat_reduced <- optim_res_reduced$par[1]
sigma_sq_hat_reduced <- optim_res_reduced$par[2]

接下来，我们可以写出该检验的对数似然比统计量LRT：

lrt_stat <- -2*(log_likelihood_reduced(c(mu_hat_reduced, sigma_sq_hat_reduced), x, y) - log_likelihood(c(mu1_hat, mu2_hat, sigma1_sq_hat, sigma2_sq_hat), x, y))

然后，我们可以使用蒙特卡罗方法重复(1)-(4)，重复次数为m = 10000，画出对数似然比统计量LRT的抽样分布：

m <- 10000
lrt_stat_sim <- numeric(m)
for (i in 1:m) {
  x_sim <- rnorm(n, mu1, sigma1)
  y_sim <- rnorm(l, mu2, sigma2)
  optim_res_sim <- optim(c(mu1=0, mu2=0, sigma1_sq=1, sigma2_sq=1), log_likelihood, x=x_sim, y=y_sim, method="L-BFGS-B", lower=c(-Inf, -Inf, 0, 0))
  mu1_hat_sim <- optim_res_sim$par[1]
  mu2_hat_sim <- optim_res_sim$par[2]
  sigma1_sq_hat_sim <- optim_res_sim$par[3]
  sigma2_sq_hat_sim <- optim_res_sim$par[4]
  lrt_stat_sim[i] <- -2*(log_likelihood_reduced(c(mu_hat_reduced, sigma_sq_hat_reduced), x_sim, y_sim) - log_likelihood(c(mu1_hat_sim, mu2_hat_sim, sigma1_sq_hat_sim, sigma2_sq_hat_sim), x_sim, y_sim))
}

ggplot(data.frame(lrt=lrt_stat_sim), aes(x=lrt)) + geom_density() + theme_classic()

如果我们将样本量增大到n=100, l=150，重复(5)中的步骤，我们可以看到对数似然比统计量LRT的抽样分布更加接近于标准正态分布，这是因为大样本会导致LRT的抽样分布更加稳定。因此，我们可以看到与(5)中的抽样分布相比，(6)中的抽样分布更加钟形对称。