R：如何计算截断正态分布的均值和协方差

Posted 2023-03-24

【中文标题】R：如何计算截断正态分布的均值和协方差

【英文标题】：R: how to compute the mean and covariance of a truncated normal distribution

【发布时间】：2021-08-21 20:19:36

【问题描述】：

我有兴趣找到截断的正态随机向量的均值和协方差。假设Y 是一个包含[Y1 Y2 Y3] 的向量。 Y 遵循具有以下均值和协方差的多元正态分布：

mu <- c(0.5, 0.5, 0.5)
sigma <- matrix(c(  1,  0.6, 0.3,
                    0.6,    1, 0.2,
                    0.3,  0.2,   2), 3, 3)

截断区域是Ys 的集合，使得AY >= 0。例如，

A <- matrix(c(1, -2, -0.5, 1.5, -2, 0, 3, -1, -1, 4, 0, -2), byrow = TRUE, nrow = 4)
> A
     [,1] [,2] [,3]
[1,]  1.0   -2 -0.5
[2,]  1.5   -2  0.0
[3,]  3.0   -1 -1.0
[4,]  4.0    0 -2.0

对于Y的以下抽奖，不满足AY >= 0：

set.seed(3)
Y <- rmvnorm(n = 1, mean = mu, sigma = sigma)
> all(A %*% as.matrix(t(Y)) >= 0)
[1] FALSE

但是对于Y 的其他绘制，它们将满足AY >= 0，我想找到满足AY >= 0 的那些Ys 的均值和协方差。

R 中存在计算截断正态分布的均值和协方差的现有包。例如，mtmvnorm 来自 tmvtnorm 包：

library(tmvtnorm)
mtmvnorm(mu, sigma, lower = ???, upper = ???)

但是，我拥有的截断集，即满足AY >= 0 的Ys 集，不能仅用lower 和upper 边界来描述。 R 是否有另一种方法来计算截断法线的均值和协方差？

【问题讨论】：

或许你可以在CrossValidated问这个问题？

@MartinGal 我一开始是这样做的。但是人们投票结束了这个问题..

【参考方案1】：

您正确理解（或可能注意到）这是不是截断多元正态分布。您将 AY>=0 作为对 Y 的线性约束，而不是简单的逐元素下限/上限。

---

如果您不是数学专家，即追求均值和协方差的显式解决方案，我想一种直接有效的方法是使用 蒙特卡洛模拟。

更具体地说，您可以假设足够大的N 来生成足够大的样本集Y，然后过滤掉满足约束AY>=0 的样本。反过来，您可以计算所选样本的均值和协方差。尝试如下

N <- 1e7
Y <- rmvnorm(n = N, mean = mu, sigma = sigma)
Y_h <- subset(Y, colSums(tcrossprod(A, Y) >= 0) == nrow(A))
mu_h <- colMeans(Y_h)
sigma_h <- cov(Y_h)

你会看到

> mu_h
[1]  0.8614791 -0.1365222 -0.3456582
> sigma_h
          [,1]       [,2]       [,3]
[1,] 0.5669915 0.29392671 0.37487421
[2,] 0.2939267 0.36318397 0.07193513
[3,] 0.3748742 0.07193513 1.37194669

---

另一种方式遵循类似的想法，但我们可以假设所选样本的集合大小，即N样本Y都应该使AY>=0站立。然后我们可以使用while 循环来做到这一点

N <- 1e6
Y_h <- list()
nl <- 0
while (nl < N) 
  Y <- rmvnorm(n = N, mean = mu, sigma = sigma)
  v <- subset(Y, colSums(tcrossprod(A, Y) >= 0) == nrow(A))
  nl <- nl + nrow(v)
  Y_h[[length(Y_h) + 1]] <- v

Y_h <- head(do.call(rbind, Y_h), N)
mu_h <- colMeans(Y_h)
sigma_h <- cov(Y_h)

你会看到

> mu_h
[1]  0.8604944 -0.1364895 -0.3463887
> sigma_h
          [,1]       [,2]       [,3]
[1,] 0.5683498 0.29492573 0.37524248
[2,] 0.2949257 0.36352022 0.07252898
[3,] 0.3752425 0.07252898 1.37427521

注意：第二个选项的优势在于，它可以为您提供足够多的选定Y_h。

【讨论】：

谢谢。有没有办法加快第二种解决方案？我发布了关于并行化 while 循环的问题（它确实在尝试加快第二个解决方案）。我已经用future 尝试了解决方案，但我想知道是否还有其他方法可以加快速度。

@Adrian 我看到了你的问题。我不认为future 可以加快你的代码。你可以试试我更新的答案***.com/a/67863714/12158757

@Adrian 我也更新了这个问题的答案，所以它不再慢了。