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

Question

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

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

当两个分布的标准差等于1（协方差矩阵在主对角线上只有一个）时，结果是相同的。但是，如果将协方差矩阵中的两个条目更改为两个或三个......结果不完全相同。

（我希望）我已经正确阅读了帮助文档，并且还有这个文档（https://cran.r-project.org/web/packages/mvtnorm/vignettes/MVT_Rnews.pdf）。

这里在stackoverflow（How to calculate multivariate normal distribution function in R）我发现这是因为也许我的协方差矩阵是错误的定义。

但是到现在为止我无法找到答案...

所以我的问题：为什么我的代码返回不同的结果，当标准差不等于一个？

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

非常感谢提前！

现在我的代码：

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来判断。