返回不同结果的二元正态分布的密度（pdf）的两个计算公式

Question

使用代码计算二元正态分布的密度。在这里，我使用了两个公式，它们应该 return 相同的结果。

第一个公式使用了mvtnorm包的dmvnorm，第二个公式使用了维基百科的公式(https://en.wikipedia.org/wiki/Multivariate_normal_distribution)。

当两个分布的标准差都为1时（协方差矩阵的主对角线上只有1），结果是一样的。但是，当您将协方差矩阵中的两个条目更改为二分之一或三分之一时……结果并不完全相同。

（我希望）我已经正确阅读了帮助以及本文档 (https://cran.r-project.org/web/packages/mvtnorm/vignettes/MVT_Rnews.pdf)。

我在 Whosebug () 上找到了这个，因为也许我的协方差矩阵定义错误。

但是直到现在我都找不到答案……

所以我的问题是： 为什么我的代码 return 在标准偏差不等于 1 时得到不同的结果？

我希望我提供了足够的信息...但是如果缺少某些内容，请发表评论。我会编辑我的问题。

非常感谢！

现在我的代码：

   library(mvtnorm)  # for loading the package if necessary

    mu=c(0,0)
    rho=0
    sigma=c(1,1)  # the standard deviation which should be changed to two or one third or… to see the different results
    S=matrix(c(sigma[1],0,0,sigma[2]),ncol=2,byrow=TRUE)

    x=rmvnorm(n=100,mean=mu,sigma=S)
    dim(x)  # for control
    x[1:5,]  # for visualization

    # defining a function
    Comparison=function(Points=x,mean=mu,sigma=S,quantity=4) {
    for (i in 1:quantity) {
           print(paste0("The ",i," random point"))
           print(Points[i,])
           print("The following two results should be the same")
           print("Result from the function 'dmvnorm' out of package 'mvtnorm'")
           print(dmvnorm(Points[i,],mean=mu,sigma=sigma,log=FALSE))
           print("Result from equation out of wikipedia")
           print(1/(2*pi*S[1,1]*S[2,2]*(1-rho^2)^(1/2))*exp((-1)/(2*(1-rho^2))*(Points[i,1]^2/S[1,1]^2+Points[i,2]^2/S[2,2]^2-(2*rho*Points[i,1]*Points[i,2])/(S[1,1]*S[2,2]))))
           print("----")
           print("----")
    } # end for-loop     
    } # end function

    # execute the function and compare the results
    Comparison(Points=x,mean=mu,sigma=S,quantity=4)

Answer 1

请记住，S 是方差-协方差矩阵。您从维基百科使用的公式使用标准偏差而不是方差。因此，您需要将对角线项的平方根代入公式。这也是为什么当你选择 1 作为对角线条目时它起作用的原因（方差和 SD 都是 1）。

查看下面修改后的代码：

 library(mvtnorm)  # for loading the package if necessary

 mu=c(0,0)
 rho=0
 sigma=c(2,1)  # the standard deviation which should be changed to two or one      third or… to see the different results
 S=matrix(c(sigma[1],0,0,sigma[2]),ncol=2,byrow=TRUE)

 x=rmvnorm(n=100,mean=mu,sigma=S)
 dim(x)  # for control
 x[1:5,]  # for visualization

 # defining a function
 Comparison=function(Points=x,mean=mu,sigma=S,quantity=4) {
   for (i in 1:quantity) {
     print(paste0("The ",i," random point"))
     print(Points[i,])
     print("The following two results should be the same")
     print("Result from the function 'dmvnorm' out of package 'mvtnorm'")
     print(dmvnorm(Points[i,],mean=mu,sigma=sigma,log=FALSE))
     print("Result from equation out of wikipedia")
     SS <- sqrt(S)
     print(1/(2*pi*SS[1,1]*SS[2,2]*(1-rho^2)^(1/2))*exp((-1)/(2*(1-rho^2))*(Points[i,1]^2/SS[1,1]^2+Points[i,2]^2/SS[2,2]^2-(2*rho*Points[i,1]*Points[i,2])/(SS[1,1]*SS[2,2]))))
     print("----")
    print("----")
  } # end for-loop     
} # end function

# execute the function and compare the results
Comparison(Points=x,mean=mu,sigma=S,quantity=4)

所以你在定义 sigma 时的评论是不正确的。在您的代码中，sigma 是方差，如果您根据构造 S.

的方式来判断，则不是标准差

返回不同结果的二元正态分布的密度（pdf）的两个计算公式

Two calculation formulas of density (pdf) of a bivariate normal distribution returning different results

r

normal-distribution

probability-density