R语言二元蒙特卡洛积分代码报错，求修复（计算二元正态分布概率）

问题背景 #

需要计算服从二元正态分布的变量落在矩形区域0<X<4、0<Y<1内的概率，分布参数为：

• X~N(3,16)（均值3，方差16）
• Y~N(1,1)（均值1，方差1）
• 变量间相关系数为-2/3
  样本量设定为100000，原代码运行后出现报错与警告。

原代码及错误信息 #

原代码 #

for(i in 1:100000){
X[i] = runif(100000)
Y[i] = runif(100000)
Z[i] = runif(100000)
 if((3/(8*(pi)*sqrt(5)))*exp((-9/10)*(((X[i])^2/16)-((17/24)*(X[i]))+(X[i]*Y[i])/3)+((Y[i])^2)-(3*(Y[i]))+41/16) < Z[i])
count = count + 1}

错误与警告信息 #

Error: object 'Z' not found
In addition: Warning messages:
1: In X[i] <- runif(1e+05) :
  number of items to replace is not a multiple of replacement length
2: In Y[i] <- runif(1e+05) :
  number of items to replace is not a multiple of replacement length

错误原因分析 #

1. 变量未初始化：X、Y、Z、count均未预先定义，R中直接向未声明的向量元素赋值会触发报错，且运行效率极低。
2. 循环内样本生成逻辑错误：每次循环调用runif(100000)生成10万个随机数，却只赋值给单个向量元素X[i]，导致"替换长度不匹配"的警告。
3. 蒙特卡洛方法误用：原代码试图用接受-拒绝抽样，但实现逻辑完全错误——对于二元正态分布，直接生成符合要求的样本是更高效准确的方式，无需手动编写复杂的联合密度公式做判断。

修复后的代码及说明 #

正确实现步骤 #

• 定义二元正态分布的均值向量与协方差矩阵
• 直接生成指定数量的二元正态样本
• 统计落在目标区域的样本数，除以总样本量得到概率

修复代码 #

# 若未安装MASS包，先运行 install.packages("MASS")
library(MASS)

# 定义二元正态分布参数
mu <- c(3, 1)  # 均值向量
# 构建协方差矩阵：Cov(X,Y)=ρ*σ_X*σ_Y，其中σ_X=4，σ_Y=1，ρ=-2/3
sigma <- matrix(c(16, (-2/3)*4*1, (-2/3)*4*1, 1), nrow = 2)

# 设置样本量
n <- 100000

# 生成二元正态分布样本
samples <- mvrnorm(n = n, mu = mu, Sigma = sigma)
X <- samples[, 1]
Y <- samples[, 2]

# 统计落在目标区域的样本数量
count <- sum(X > 0 & X < 4 & Y > 0 & Y < 1)

# 计算概率
probability <- count / n
print(probability)

代码说明 #

1. 协方差矩阵计算：二元正态分布的协方差矩阵中，非对角线元素为变量间的协方差，由相关系数×X的标准差×Y的标准差计算得出。
2. 样本生成：MASS包的mvrnorm函数是生成多元正态样本的高效工具，避免手动实现的复杂逻辑与错误。
3. 向量运算替代循环：R中向量运算的效率远高于for循环，统计计数时直接用逻辑向量求和，简洁清晰。

无第三方包替代方案 #

若不想依赖MASS包，可使用mvtnorm包实现：

# 若未安装先运行 install.packages("mvtnorm")
library(mvtnorm)

samples <- rmvnorm(n = n, mean = mu, sigma = sigma)
# 后续统计步骤与上述代码一致

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