while (e_i$X1 < 12 | e_i$X2 < 12) { 错误:参数长度为零

Error in while (e_i$X1 < 12 | e_i$X2 < 12) { : argument is of length zero

在之前的问题()中,我学习了如何进行以下模拟:

第一步:不断生成两个随机数“a”和“b”,直到“a”和“b”都大于12

第 2 步:跟踪在完成第 1 步之前必须生成多少个随机数

步骤 3: 重复步骤 1 和步骤 2 100 次

res <- matrix(0, nrow = 0, ncol = 3)    

for (j in 1:100){
  a <- rnorm(1, 10, 1)
  b <- rnorm(1, 10, 1)
  i <- 1
  while(a < 12 | b < 12) {
    a <- rnorm(1, 10, 1)
    b <- rnorm(1, 10, 1)
    i <- i + 1
  }
  x <- c(a,b,i)
  res <- rbind(res, x)
}

head(res)
      [,1]     [,2] [,3]
x 12.14232 12.08977  399
x 12.27158 12.01319 1695
x 12.57345 12.42135  302
x 12.07494 12.64841  600
x 12.03210 12.07949   82
x 12.34006 12.00365  782

问题:现在,我想对上面的代码稍作修改——而不是单独生成“a”和“b”,我希望它们是“一起”产生(用数学术语来说:“a”和“b”是由两个独立的单变量正态分布产生的,现在我希望它们来自双变量正态分布)。

我尝试自己修改这段代码:

library(MASS)

Sigma = matrix(
    c(1,0.5, 0.5, 1), # the data elements
    nrow=2,              # number of rows
    ncol=2,              # number of columns
    byrow = TRUE)        # fill matrix by rows

res <- matrix(0, nrow = 0, ncol = 3)    

for (j in 1:100){
    e_i = data.frame(mvrnorm(n = 1, c(10,10), Sigma))
    e_i$i <- 1
    while(e_i$X1 < 12 | e_i$X2 < 12) {
        e_i = data.frame(mvrnorm(n = 1, c(10,10), Sigma))
        e_i$i <- i + 1
    }
    x <- c(e_i$X1, e_i$X2  ,i)
    res <- rbind(res, x)
}

res = data.frame(res)

但这会产生以下错误:

Error in while (e_i$X1 < 12 | e_i$X2 < 12) { : argument is of length zero

如果我对您的代码的理解正确,您是在尝试查看在两个值 >=12 之前出现了多少样本,并进行了 100 次试验?这是我会采用的方法:

library(MASS)

for(i in 1:100){
  n <- 1
  while(any((x <- mvrnorm(1, mu=c(10,10), Sigma=diag(0.5, nrow=2)+0.5))<12)) n <- n+1
  if(i==1) res <-            data.frame("a"=x[1], "b"=x[2], n)
  else     res <- rbind(res, data.frame("a"=x[1], "b"=x[2], n))
}

在这里,我在 while() 调用中将 mvrnorm 的结果分配给 x。在同一个调用中,它使用 any() 函数评估其中一个是否小于 12。如果计算结果为 FALSEn(计数器)增加并重复该过程。一旦 TRUE,这些值将附加到您的 data.frame 并返回到 for 循环的开始。


关于您的代码,mvrnorm() 函数在 n=1 时返回一个向量,而不是矩阵,因此两个值都进入 data.frame 中的单个变量:

data.frame(mvrnorm(n = 1, c(10,10), Sigma))

Returns:

  mvrnorm.n...1..c.10..10...Sigma.
1                         9.148089
2                        10.605546

您的 data.frame() 调用中的 matrix() 函数以及对 i 的一些调整将修复您的代码:

library(MASS)

Sigma = matrix(
  c(1,0.5, 0.5, 1), # the data elements
  nrow=2,              # number of rows
  ncol=2,              # number of columns
  byrow = TRUE)        # fill matrix by rows

res <- matrix(0, nrow = 0, ncol = 3)    

for (j in 1:10){
  e_i = data.frame(matrix(mvrnorm(n = 1, c(10,10), Sigma), ncol=2))
  i <- 1
  while(e_i$X1[1] < 12 | e_i$X2[1] < 12) {
    e_i = data.frame(matrix(mvrnorm(n = 1, c(10,10), Sigma), ncol=2))
    i <- i + 1
  }
  x <- c(e_i$X1, e_i$X2  ,i)
  res <- rbind(res, x)
}

res = data.frame(res)