r - 查找模拟的摘要（即平均值、中位数、众数、se 等），以及该模拟在 R 中的分布

Question

我有一个问题，我生成了一个模拟，其中你基本上在 R 中有一个 100x100 的绘图并从中心 (50,50) 开始，然后朝着一个方向迈出一步，一步一步地尝试到达外面。

我能够生成循环来执行此操作，但现在我想运行这个特定的循环 10,000 次，然后为其生成摘要并查看分布。我只是不确定如何运行循环 10,000 次然后为其生成摘要或看起来如何。我已经包括了到目前为止的内容：

plot(0:100,0:100,type="n")
points(50,50,col="red",pch=16,cex=1.5)
x0<-50
y0<-50
x1<-sample(c(-1,0,1),1)
y1<-sample(c(-1,0,1),1)


for(i in 1:50000){
  x1<-sample(c(-1,0,1),1)
  y1<-sample(c(-1,0,1),1)
  lines(c(x0,x0+1),c(y0,y0+1))
  x0<-x0+x1
  y0<-y0+y1
  if(x0>100|x0<0|y0>100|y0<0)break
}

标签： r

Answer 1

由于我们真正想知道的关于模拟的唯一事情是到达外部需要多少步，我们将从创建一个模拟路径并仅返回步数的函数开始。

random_path_length <- function() {
  x <- 50
  y <- 50
  res <- 0
  while (x >= 0 & x <= 100 & y >= 0 & x <= 100) {
    dx <- sample(c(-1, 0, 1), 1)
    dy <- sample(c(-1, 0, 1), 1)
    x <- x + dx
    y <- y + dy
    res <- res + 1
  }
  res
}

set.seed(1)
random_path_length()
#> [1] 3210

理论上我们可以使用这个函数来模拟 10000 个结果，但问题是它需要很多时间。我建议不是逐个模拟步骤，而是分批模拟步骤以利用矢量化操作。

random_path_length <- function(batch_size = 1000) {
  x <- 50
  y <- 50
  res <- 0
  go <- TRUE
  while (go) {
    # simulate batch_size number of steps
    dx <- sample(c(-1, 0, 1), batch_size, replace = TRUE)
    dy <- sample(c(-1, 0, 1), batch_size, replace = TRUE)
    new_x <- x + cumsum(dx)
    new_y <- y + cumsum(dy)
    # stop path at the point when (if) it reaches the outside
    where_reaches <- which(new_x == 0 | new_x == 100 | new_y == 0 | new_y == 100)
    # this batch didn't reach the outside
    if (length(where_reaches) == 0) {
      res <- res + batch_size
      x <- new_x[[batch_size]]
      y <- new_y[[batch_size]]
    } else {
      where_stop <- where_reaches[[1]]
      res <- res + where_stop
      x <- new_x[[where_stop]]
      y <- new_y[[where_stop]]
      go <- FALSE
    }
  }
  res
}

set.seed(1)
random_path_length()
#> [1] 3023

这里batch_size的参数调节我们一次生成多少步。您可以使用它来找出哪个效果更快。

有了这个函数，我们可以在可接受的时间内模拟 10000 个结果，然后继续计算我们想要的任何统计数据。

set.seed(1)
res <- replicate(10000, random_path_length())
summary(res)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>     229    1110    1784    2210    2858   12731

hist(res)