如何在R中计算多元正态分布函数

Question

以下是我尝试的内容，使用了mvtnorm软件包

样本数据集

library(mvtnorm)

set.seed(2357)
df <- data.frame(
  x = rnorm(1000, mean=80, sd=20),
  y = rnorm(1000, mean=0, sd=5),
  z = rnorm(1000, mean=0, sd=5)
)

head(df)
      x      y       z
1 70.38  1.307  0.2005
2 59.76  5.781 -3.5095
3 54.14 -1.313 -1.9022
4 79.91  7.754 -6.2076
5 87.07  1.389  1.1065
6 75.89  1.684  6.2979

适应多元正态分布并检查P(x <= 80) ~ 0.5

# Get the dimension means and correlation matrix
means <- c(x=mean(df$x), y=mean(df$y), z=mean(df$z))
corr <- cor(df)

# Check P(x <= 80)
sum(df$x <= 80)/nrow(df)  # 0.498
pmvnorm(lower=-Inf, upper=c(80, Inf, Inf), mean=means, corr=corr)  # 0.8232

为什么拟合结果是0.82？我哪里出错了？

- Ben

Answer 1

首先，您不需要模拟任何内容来研究pmvnorm函数：

pmvnorm(lower=rep(-Inf, 3), upper=c(80, Inf, Inf), mean=c(80,0,0), corr=diag(rep(1,3)))

结果是 0.5，正如您所预期的那样。

您的均值向量大约为 (79, 0, 0)，因此让我们尝试一下：

pmvnorm(lower=rep(-Inf, 3), upper=c(80, Inf, Inf), mean=c(79,0,0), corr=diag(rep(1,3)))

现在的结果是0.8413447。没有问题。通过仅指定相关矩阵，您告诉软件假定所有方差都是统一的。在您的模拟中，方差分别为400、25和25：与您在参数中指定的非常不同！正确的计算使用数据的协方差矩阵，而不是其相关矩阵：

pmvnorm(lower=rep(-Inf, 3), upper=c(80, Inf, Inf), mean=means, sigma=cov(df))

结果为0.5178412，与数据相符。

- whuber

---

@Frank，感谢您的有益编辑。我会记住这个惯例在以后的帖子中！ - whuber