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

In an earlier question (R: Logical Conditions Not Being Respected), I learned how to make the following simulation :

Step 1: Keep generating two random numbers "a" and "b" until both "a" and "b" are greater than 12

Step 2: Track how many random numbers had to be generated until it took for Step 1 to be completed

Step 3: Repeat Step 1 and Step 2 100 times

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

Question: Now, I am trying to make a slight modification to the above code - Instead of "a" and "b" being produced separately, I want them to be produced "together" (in math terms: "a" and "b" were being produced from two independent univariate normal distributions, now I want them to come from a bivariate normal distribution).

I tried to modify this code myself:

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)

But this is producing the following error:

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

CodePudding user response：

If I understand your code correctly you are trying to see how many samples occur before both values are >=12 and doing that for 100 trials? This is the approach I would take:

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

Here I am assigning the results of a mvrnorm to x within the while() call. In that same call, it evaluates whether either are less than 12 using the any() function. If that evaluates to FALSE, n (the counter) is increased and the process repeated. Once TRUE, the values are appended to your data.frame and it goes back to the start of the for-loop.

Regarding your code, the mvrnorm() function is returning a vector, not a matrix, when n=1 so both values go into a single variable in the 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

The matrix() function within your data.frame() calls, along with some tweaks to your use of i, will fix your code:

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)