如何从 MVN 模型生成样本？答案

【问题标题】：How to generate samples from MVN model?如何从 MVN 模型生成样本？
【发布时间】：2020-02-25 17:33:18
【问题描述】：

我正在尝试根据本文here 通过示例 5.1 在 R 上运行一些代码。我想模拟以下内容：

我的 R 背景不是很好，所以我有以下代码，如何从中生成直方图和样本？

xseq<-seq(0, 100, 1)  
n<-100
Z<- pnorm(xseq,0,1)
U<- pbern(xseq, 0.4, lower.tail = TRUE, log.p = FALSE)
Beta <- (-1)^U*(4*log(n)/(sqrt(n)) + abs(Z))

【问题讨论】：

您可以使用 MASS 包中的 mvrnorm 模拟多元正态分布。 stat.ethz.ch/R-manual/R-devel/library/MASS/html/mvrnorm.html
@jandraor 谢谢，但我怎样才能改变上面的系数？
为了模拟单变量随机正态样本，您应该使用 rnorm，而不是 pnorm。在伯努利分布的情况下，您可以使用 n = 1 的 rbinom 对其进行模拟。查看 mvrnorm 文档以生成示例。您需要指定均值和协方差矩阵。
您能告诉我们有关 $X_i$ 集合中变量的所需协方差矩阵的任何信息吗？（例如，它们应该是 i.i.d. 标准正态变量吗？）
@AaronMontgomery，生成的 i.i.d.按标准多元正态，是的

标签： r montecarlo multivariate-testing

【解决方案1】：

一些有用的工具演示：

rnorm(1)                 # generates one standard normal variable
rnorm(10)                # generates 10 standard normal variables
rnorm(1, 5, 6)           # generates 1 normal variable with mu = 5, sigma = 6
                         # not needed for this problem, but perhaps worth saying anyway

rbinom(5, 1, 0.4)        # generates 5 Bernoulli variables that are 1 w/ prob. 0.4

所以，要生成一个 beta 实例：

n <- 100                 # using the value you gave; I have no idea what n means here
u <- rbinom(1, 1, 0.4)   # make one Bernoulli variable
z <- rnorm(1)            # make one standard normal variable
beta <- (-1)^u * (4 * log(n) / sqrt(n) + abs(z))

但现在，您希望多次执行此操作以进行蒙特卡洛模拟。您可以这样做的一种方法是构建一个函数，将beta 作为其输出，并使用replicate() 函数，如下所示：

n <- 100                    # putting this here because I assume it doesn't change
genbeta <- function(){      # output of this function will be one copy of beta
  u <- rbinom(1, 1, 0.4)
  z <- rnorm(1)
  return((-1)^u * (4 * log(n) / sqrt(n) + abs(z)))
}

# note that we don't need to store beta anywhere directly; 
# rather, it is just the return()ed value of the function we defined

betadraws <- replicate(5000, genbeta())
hist(betadraws)

这将产生 5000 个 beta 变量副本并将它们放入直方图中的效果。

还有其他方法可以做到这一点——例如，可以只制作一个随机变量的大矩阵并直接使用它——但我认为这将是最清晰的开始方法。

编辑：我意识到我完全忽略了第二个等式，这可能是你不想要的。

我们现在已经制作了一个beta 值的向量，您可以在上面的replicate() 函数的第一个参数中控制向量的长度。在下面的后续示例中，我将其保留为 5000。

要获取 Y 向量的随机样本，您可以使用类似：

x <- replicate(5000, rnorm(17))     
  # makes a 17 x 5000 matrix of independent standard normal variables
epsilon <- rnorm(17)
  # vector of 17 standard normals
y <- x %*% betadraws + epsilon
  # y is now a 17 x 1 matrix (morally equivalent to a vector of length 17)

如果您想获得其中的许多，您可以将那个包装在另一个函数中并 replicate() 它。

或者，如果您不想要 Y 向量，而只想要一个 Y_i 分量：

x <- rnorm(5000)    
  # x is a vector of 5000 iid standard normal variables 
epsilon <- rnorm(1)
  # epsilon_i is a single standard normal variable
y <- t(x) %*% betadraws + epsilon
  # t() is the transpose function; y is now a 1 x 1 matrix

【讨论】：