如何使用R通过蒙特卡洛模拟计算卡方拟合优度检验功效函数

原代码的核心问题 #

你写的错误代码完全偏离了蒙特卡洛功效计算的核心逻辑，问题点如下：

• 没有按循环生成新的随机样本：t检验的模拟中，每次内层循环都会从真实总体抽取新的样本做检验，但你的代码每次循环都用固定的Frequency_true做卡方检验，还生成了完全没有被调用的rchisq随机数，1000次循环跑的是完全相同的检验，结果只会是固定值，没有任何随机波动。
• 抽样分布选错：卡方拟合优度检验面对的是分类计数数据，应该从多项分布抽样生成样本，而非卡方分布。
• 循环逻辑无意义：原假设概率向量Frequency_0是固定的等概率向量，你按它的长度做外层循环，相当于重复计算5次一模一样的结果，没有遍历任何可变参数。
• 频数和概率混淆：你定义的Frequency_true是频数（求和为250），不是概率，不能直接用于抽样。

实现逻辑（完全对齐你已掌握的t检验模拟框架） #

功效计算的本质是：在备择假设成立的前提下，重复抽样、做检验，统计「拒绝原假设」的次数占比，步骤和t检验完全对应：

• 固定全局参数：总样本量、真实总体的分类概率（备择假设分布）、原假设分类概率、显著性水平α、蒙特卡洛重复次数B
• 外层循环：遍历你想要考察的可变参数（比如不同样本量、不同效应量、不同原假设设定，对应t检验里遍历不同的原假设均值）
• 内层循环：
  1. 从真实总体分布（多项分布）抽取1组样本，得到各类别的观测频数
  2. 对这组样本执行卡方拟合优度检验，对比原假设的概率
  3. 记录检验p值是否小于α（即是否拒绝原假设）
• 内层循环结束后，当前参数下的经验功效 = 拒绝原假设的次数 / 总重复次数B

可直接运行的正确代码 #

下面的代码实现了「遍历不同样本量，计算对应卡方拟合优度检验功效」的流程，和你提供的t检验代码结构完全一致：

# 全局参数设置
p_true <- c(50,60,40,47,53) / sum(c(50,60,40,47,53)) # 真实总体的分类概率
p_null <- rep(0.2, 5) # 原假设：5个类别等概率
alpha <- 0.05 # 显著性水平
B <- 1000 # 蒙特卡洛重复次数

# 外层遍历参数：这里以遍历样本量为例，你也可以替换成遍历不同效应量
n_seq <- seq(100, 500, length.out = 10)
Empirical_Power <- rep(NA, length(n_seq))

for(j in 1:length(Empirical_Power)){
  current_n <- n_seq[j]
  Test_Decisions <- rep(NA, B)
  
  for(i in 1:B){
    # 核心步骤：从真实多项分布生成1组样本计数
    sample_counts <- as.vector(rmultinom(n = 1, size = current_n, prob = p_true))
    # 对当前样本做卡方拟合优度检验
    chisq_res <- chisq.test(x = sample_counts, p = p_null)
    # 记录检验决策：TRUE=拒绝原假设，FALSE=不拒绝
    Test_Decisions[i] <- chisq_res$p.value < alpha
  }
  # 计算当前参数下的经验功效
  Empirical_Power[j] <- mean(Test_Decisions)
}

# 可视化功效曲线
plot(n_seq, Empirical_Power, type = "b", 
     xlab = "总样本量", ylab = "经验功效",
     main = "卡方拟合优度检验功效随样本量变化曲线")

提示：如果抽样后某类别的期望频数小于5，卡方检验的渐近结果会有偏差，可以在chisq.test()中加入参数simulate.p.value = TRUE，用蒙特卡洛方法计算精确p值，提升结果可靠性。

内容的提问来源于stack exchange，提问作者Omar Elzoghby