如何计算 R 中的多元正态分布函数
How to calculate multivariate normal distribution function in R
这是我尝试过的，利用 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？我哪里错了？
首先，你不需要模拟任何东西来研究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，和数据很吻合
这是我尝试过的，利用 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？我哪里错了？
首先，你不需要模拟任何东西来研究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，和数据很吻合