0

假设我想生成 X 和概率 (F) 的前一个对象的 X 和 Y 条件。我的算法：在特定的一天，一个人的状态为 1/0。如果状态为 0，则第二天的状态为rbinom(1,1,0.5)；如果状态为 1，那么第二天的状态将是rbinom(1,1,F[this day])。在这里，我想为 N 个人和 d 天生成此数据。

我的代码：

N <- 10 # number of individual
d <- 10 # days
F <- c(0.9,0.8, 0.7, 0.6, 0.5, 0.4) # Probabilities

    for(j in 1:N){
        x <- NULL
        x[1] <- rbinom(1,1,0.5)
        y <- NULL
        y[1] <- ifelse(x[1]==1, 1, 0)
        for(i in 2:d){ 
            if(x[i-1]==0)
            {
            y[i] <- x[i-1]
            x[i] <- rbinom(1,1,0.05)
            }
            if(x[i-1]==1) 
            {
            y[i] <- y[i-1]+x[i-1]
            x[i] <- rbinom(1,1,F[y[i]]) 
            }
        }

    }

这段代码产生的数据像，

   > x
     [1] 1 1 0 0 0 0 0 0 0 1
    > y
     [1] 1 2 3 0 0 0 0 0 0 0

但我想像这样生成 X 和 Y：

 > x
 [1] 1 1 0 0 0 0 0 0 0 1
> y
 [1] 1 2 0 0 0 0 0 0 0 1

我不知道代码的哪一部分是错误的。有人可以帮忙吗？

4

1 回答 1

0

可以删除外部循环，然后您可以使用该replicate函数在每次迭代时存储输出。我还将该过程扩展到 20 个步骤，以查看您是否可以看到获得两个“事件”的实例。

N <- 10 
d <- 20
F <- c(0.9,0.8, 0.7, 0.6, 0.5, 0.4) 
   replicate(10,{     x <- NULL
        x[1] <- rbinom(1,1,0.5)
        y <- NULL
        y[1] <- ifelse(x[1]==1, 1, 0)
        for(i in 2:d){ 
            if(x[i-1]==0)
            {
            y[i] <- x[i-1]
            x[i] <- rbinom(1,1,0.05)
            }
            if(x[i-1]==1) 
            {
            y[i] <- y[i-1]+x[i-1]
            x[i] <- rbinom(1,1,F[y[i]]) 
            }
        }
      return(rbind(x,y)) } )

经过几次运行后，我确实看到你得到了（我认为）你想要的东西，即 x 和 y 中的两个正值系列。这恰好是 10 次中的第 7 次：

, , 7

  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20]
x    0    0    0    0    0    0    0    0    1     1     1     0     0     0     1     1     1     1     0     0
y    0    0    0    0    0    0    0    0    0     1     2     3     0     0     0     1     2     3     4     0

你从来没有用自然语言描述过你当前的输出有什么不满意的地方，所以我想这仍然只是猜测。我还认为您应该将模拟结果分配给一个对象，以便稍后在 N 变得足够高以至于无法通过肉眼查看时以编程方式进行检查。

res <- replicate( 1000, { code block here } )
于 2015-04-02T16:39:28.083 回答