R模拟中计算ED值并统计ED>7.815的数量及占比求助

我已有一段能生成随机样本并重复大量模拟的R代码，运行正常，可产出计算新变量ED所需的数据。ED的计算公式为：ED = N²*(theta_new)%*%sigma_inverse%*%trans。我需要在每次模拟中获取并打印ED值，随后统计所有模拟中ED>7.815的总数及其占比，但编写的代码运行异常，代码如下：

# Initialize the counter
  n_greater_total <- 0
  
  # Calculate the ED value for each simulation
  ED <- c()
  for (i in 1:n_simulations) {
    
    # Calculate the ED value for the full sample
    ED_current <- N^2* (theta_new)%*%sigma_inverse %*% trans
    ED <- c(ED, ED_current)
    
    # Print the ED value
    print(paste("ED for simulation", i, "is", ED_current))
    
    # Check if the ED value is greater than 7.815
    if (ED_current >= 7.815) {
      
      # If yes, increment the counter
      n_greater_total <- n_greater_total + 1
    }
  }
  
  # Calculate the total number of ED values greater than 7.815
  total_greater <- sum(ED > 7.815)
  
  # Calculate the percentage of ED values greater than 7.815
  percentage <- (total_greater / n_simulations) * 100
  
  # Print the results
  print("Total number of values greater than 7.815:")
  print(total_greater)
  
  print("Percentage of ED values  greater than 7.815:")
  print(percentage)

核心问题 #

1. 循环内未更新关键变量：theta_new、sigma_inverse、trans未在循环内重新计算，导致所有模拟复用同一组值，不符合模拟逻辑。
2. 矩阵运算结果类型问题：%*%返回矩阵对象，直接用于打印和判断会引发类型错误，需转为数值。
3. 动态扩展向量效率低：ED <- c(ED, ED_current)每次循环都重新分配内存，模拟次数多时性能差。

修正后的完整代码 #

# 初始化模拟参数（替换为你的实际参数）
n_simulations <- 100  # 模拟次数
N <- 50               # 公式中的N值

# 预分配ED向量，提升运行效率
ED <- numeric(n_simulations)
n_greater_total <- 0

for (i in 1:n_simulations) {
  # --------------------------
  # 此处放入你的样本生成与变量计算代码
  # 示例（替换为你实际的逻辑）：
  # theta_new <- rnorm(2)
  # sigma_inverse <- matrix(runif(4), nrow=2)
  # trans <- theta_new
  # --------------------------
  
  # 计算ED并转换为数值类型
  ED_current <- as.numeric(N^2 * theta_new %*% sigma_inverse %*% trans)
  ED[i] <- ED_current
  
  # 打印当前模拟的ED值
  cat(sprintf("第%d次模拟的ED值：%.4f\n", i, ED_current))
  
  # 判断是否超过阈值
  if (ED_current > 7.815) {
    n_greater_total <- n_greater_total + 1
  }
}

# 统计并输出结果
total_greater <- sum(ED > 7.815)
percentage <- (total_greater / n_simulations) * 100

cat("\nED值大于7.815的总数：", total_greater, "\n")
cat("ED值大于7.815的占比：", sprintf("%.2f%%", percentage), "\n")

关键修正说明 #

• 将样本生成、theta_new等变量的计算逻辑移入循环，保证每次模拟使用新的样本数据。
• 用as.numeric()将矩阵运算结果转为数值，避免类型错误。
• 预分配ED向量，大幅提升大规模模拟时的运行效率。
• 用cat()替代print()，输出格式更简洁易读。

内容的提问来源于stack exchange，提问作者Amal Ghania