R中多元正态分布的数值积分

Posted 2023-03-24

【中文标题】R中多元正态分布的数值积分

【英文标题】：numerical integration of multivariate normal distribution in R

【发布时间】：2021-11-14 00:25:26

【问题描述】：

我正在尝试找到t1 and t2 的解决方案，条件如下

# x1 and x2 are bivariate normal distributed with 
mu1=mu2=0
cov=matrix(c(1,1/2,1/2,1),2,2)

# trying to find the solution satisfying the equality with  t1=t2
# p(x1>t1) + p(x1<t2,x2>t2) = 0.05

我想知道是否有人可以帮助我编写 R.THX 中的代码！

【参考方案1】：

在 R 中，您感兴趣的概率是

  pnorm(t, lower.tail = FALSE) + 
    pmvnorm(c(-Inf,t), c(t,Inf), sigma = Sigma) 



library(mvtnorm)
Sigma <- rbind(c(1,0.5),c(0.5,1))
f <- function(t)
  pnorm(t, lower.tail = FALSE) + 
    pmvnorm(c(-Inf,t), c(t,Inf), sigma = Sigma) - 0.05


uniroot(f, c(-6,6))
# $root
# [1] 1.916334

【参考方案2】：

p(x1>t) + p(x1<t,x2>t) = 0.05 与p(x1<t,x2<t) = 0.95 相同。因此：

library(mvtnorm)
qmvnorm(0.95, mean = rep(0, 2), sigma = cov)
#$quantile
#[1] 1.916399
#
#$f.quantile
#[1] 2.600961e-06
#
#attr(,"message")
#[1] "Normal Completion"

检查：

pmvnorm(c(1.916399, -Inf), c(Inf, Inf), , mean = rep(0, 2), sigma = cov) +
  pmvnorm(c(-Inf, 1.916399), c(1.916399, Inf), , mean = rep(0, 2), sigma = cov)
#[1] 0.04999262
#attr(,"error")
#[1] 2e-16
#attr(,"msg")
#[1] "Normal Completion"

【讨论】：

一个用途是做一个错误。我发现相反。

我显然错了。 P(x1 > -1.916334) 是 > 1/2。

啊不抱歉，我也发现了。