问题描述
这是我尝试过的，利用了 mvtnorm 包
Here's what I tried, making use of the mvtnorm package
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 # 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
# 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?我哪里做错了?
Why is the fitted result 0.82? Where did I go wrong?
推荐答案
首先，你不需要模拟任何东西来研究pmvnorm函数:
First, you don't need to simulate anything to study the pmvnorm function:
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)，所以让我们试试看:
Your means vector is approximately (79, 0, 0), so let's try it:
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:与您在参数中指定的非常不同！
The result now is 0.8413447. There's nothing the matter. By specifying only the correlation matrix, you told the software to assume that all variances were unity. In your simulation, the variances were 400, 25, and 25: very different from what you specified in the arguments!
正确的计算使用的是数据的协方差矩阵，而不是其相关矩阵:
pmvnorm(lower=rep(-Inf, 3), upper=c(80, Inf, Inf), mean=means, sigma=cov(df))
结果是0.5178412，与数据相当.
这篇关于如何计算R中的多元正态分布函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！