R:如何计算截断正态分布的均值和协方差
R: how to compute the mean and covariance of a truncated normal distribution
我对查找截断正态随机向量的均值和协方差很感兴趣。假设 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)
截断区域是 Y
的集合,使得 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
的Y
的均值和协方差。
R 中现有的包可以计算截断正态分布的均值和协方差。例如,来自 tmvtnorm
包的 mtmvnorm
:
library(tmvtnorm)
mtmvnorm(mu, sigma, lower = ???, upper = ???)
然而,我拥有的截断集,即满足 AY >= 0
的 Y
集,不能仅用 lower
和 upper
边界来描述。 R 是否有另一种方法来计算截断法线的均值和协方差?
您正确理解(或者可能注意到)这是不是截断多元正态分布。您将 AY>=0
作为 Y
上的线性约束,而不是简单的逐元素 lower/upper 边界。
如果您不是数学爱好者,即追求均值和协方差的显式解,我想一个简单有效的方法是使用 Monte Carlo 模拟。
更具体地说,您可以假设一个足够大的 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
随心所欲
我对查找截断正态随机向量的均值和协方差很感兴趣。假设 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)
截断区域是 Y
的集合,使得 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
的Y
的均值和协方差。
R 中现有的包可以计算截断正态分布的均值和协方差。例如,来自 tmvtnorm
包的 mtmvnorm
:
library(tmvtnorm)
mtmvnorm(mu, sigma, lower = ???, upper = ???)
然而,我拥有的截断集,即满足 AY >= 0
的 Y
集,不能仅用 lower
和 upper
边界来描述。 R 是否有另一种方法来计算截断法线的均值和协方差?
您正确理解(或者可能注意到)这是不是截断多元正态分布。您将 AY>=0
作为 Y
上的线性约束,而不是简单的逐元素 lower/upper 边界。
如果您不是数学爱好者,即追求均值和协方差的显式解,我想一个简单有效的方法是使用 Monte Carlo 模拟。
更具体地说,您可以假设一个足够大的 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
随心所欲