R：模拟正态分布中的点

Question

我正在使用 R 编程语言。

考虑以下情况 - 对于均值为 5 且 sigma=5 的正态分布：

• A) 基于数值积分，观察到大于 7 的点的概率是多少？
• B) 根据模拟（例如 1000 次模拟），观察到大于 7 的点的概率是多少？

理想情况下，A) 和 B) 的答案应该相同：

library(MASS)
mu <- 5
sigma <- 5

# nmerically integrate to find P(X > 7)
 1 - pnorm(7, mean = mu, sd = sigma)
#[1] 0.3445783

# Simulation
set.seed(123)  
simulated_points <- rnorm(1000, mean = mu, sd = sigma)

# Find the percentage of points greater than 5
sum(simulated_points > 7) / length(simulated_points) * 100
#[1] 33.7

我们可以看到，上面的结果是匹配的。

我的问题：假设现在我有一个二维多元正态分布，其中

mu1 = 1, mu2 = 2, sigma1 = 3, sigma2 = 6, sigma12 = sigma=21=5

• A) 基于数值积分，观测到大于 (2,3) 的点的概率是多少？
• B) 根据模拟（例如 1000 次模拟），观察到大于 (2,3) 的点的概率是多少？

我自己尝试这样做：

library(mvtnorm)

mu <- c(1, 2)
Sigma <- matrix(c(9, 5, 5, 36), nrow=2)  # Covariance matrix

   # Numerically integrate

p_x1_gt_2_and_x2_gt_3 <- pmvnorm(lower = -Inf, upper = c(2, 3), mean = mu, sigma = Sigma)

# Simulate 1000 points 
set.seed(123)  
simulated_points <- mvrnorm(100000, mu = mu, Sigma = Sigma)

percentage_gt_2_and_3 <- sum(simulated_points[,1] > 2 & simulated_points[,2] > 3) / nrow(simulated_points) * 100

但是，结果不匹配：

> p_x1_gt_2_and_x2_gt_3
[1] 0.3990654

> percentage_gt_2_and_3
[1] 20.307

有人可以告诉我我做错了什么以及我可以做些什么来解决这个问题吗？

谢谢！

Answer 1

您的 pmvnorm 调用是错误的，这里您正在查看从 -Inf 到 (2,3) 的密度，而您感兴趣的是从 (2,3) 到 +Inf 的密度。

> p_x1_gt_2_and_x2_gt_3 <- pmvnorm(lower = c(2, 3), upper=c(Inf,Inf), mean = mu, sigma = Sigma)
[1] 0.2023229
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"