如何从 MVN 模型生成样本？

Question

我正在尝试运行基于本文 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))

Answer 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))

但是现在，您想多次执行此操作以进行 Monte Carlo 模拟。您可以这样做的一种方法是构建一个函数，将 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)

如果你想获得其中的许多，你可以将 that 包装在另一个函数中，然后 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

如何从 MVN 模型生成样本？

How to generate samples from MVN model?

r

montecarlo

multivariate-testing