logo标签列表

模拟随机漫步

rrandom-walk
13

Xn可以取-1或1的值，每个值取到的概率都是0.5。并且Sn= Sn-1 + Xn。我该如何计算在时间n观察到的偏差总和Sn = X1 + X2 + : : : + Xn。我想在这里模拟一次随机游走。我已经尝试了以下方法，但不确定是否正确：

rw <- function(n){
    x=numeric(n)
    xdir=c(TRUE, FALSE)
    step=c(1,-1)
    for (i in 2:n)
    if (sample(xdir,1)) {
        x[i]=x[i-1]+sample(step,1)
    } else {
        x[i]=x[i-1]
    }
    list(x=x)
}
请帮忙！
1
您也可以查看R-bloggers上的这篇文章。 - Konrad
4个回答
50

您也可以使用 cumsum 来实现非常简洁和高效的方法。

set.seed(1)

n <- 1000
x <- cumsum(sample(c(-1, 1), n, TRUE))

输入图像描述

3
我认为这是处理它的方式-除了由于向量化而更快（在我的计算机上进行1000000步长的行走时快了44倍），它只需要一行代码就能完成。 - josliber
为什么你会对 c(-1,1) 进行采样，而不是正态分布 rnorm(10000,0,1)？我猜这将更好地等同于白噪声。 - pieterbons
1
@pieterbons 这个问题是关于随机游走的，每一步可以是-1或+1。 - Jake Burkhead
4
本文讨论了计算随机游走的各种基本R方法的时间。本文受到此帖子中的评论和@josilber在Jake Burkhead发表的最快方法帖子中的评论的启发。
下面使用各种方法来计算随机游走。为了实现这一点，每个函数从以下fnc定义的1000个值中提取1或-1。时间测试使用microbenchmark对每种方法进行1000次复制。
fnc <- function(n) sample(c(1L, -1L), n, replace=TRUE)
library(microbenchmark)
microbenchmark(all=cumsum(fnc(1000L)),
      reduce=Reduce("+", fnc(1000L), accumulate=TRUE),
      laplyRpCln=cumsum(unlist(lapply(rep.int(1L, 1000L), fnc))),
      laplyRpAn=cumsum(unlist(lapply(rep.int(1L, 1000L), function(x) fnc(1L)))),
      laplySqAn=cumsum(unlist(lapply(seq_len(1000L), function(x) fnc(1L)))),
      saplyRpCln=cumsum(sapply(rep.int(1L, 1000L), fnc)),
      saplyRpAn=cumsum(sapply(rep.int(1L, 1000L), function(x) fnc(1L))),
      saplySqAn=cumsum(sapply(seq_len(1000L), function(x) fnc(1L))),
      vaplyRpCln=cumsum(vapply(rep.int(1L, 1000L), fnc, FUN.VALUE=0)),
      vaplyRpAn=cumsum(vapply(rep.int(1L, 1000L), function(x) fnc(1L), FUN.VALUE=0)),
      vaplySqAn=cumsum(vapply(seq_len(1000L), function(x) fnc(1L), FUN.VALUE=0)),
      replicate=cumsum(replicate(1000L, fnc(1L))),
      forPre={vals <- numeric(1000L); for(i in seq_along(vals)) vals[i] <- fnc(1L); cumsum(vals)},
      forNoPre={vals <- numeric(0L); for(i in seq_len(1000L)) vals <- c(vals, fnc(1L)); cumsum(vals)},
      times=1000)

在这里，

  • "all" 使用了Jake Burkhead的建议，使用cumsum一次性推出样本。
  • "reduce" 一次性推出样本，但使用Reduce执行求和。
  • laplyRpCln 使用lapply和unlist返回一个向量，并迭代1000个1的实例，通过名称直接调用函数。
  • laplyRpAn 不同之处在于使用匿名函数。
  • laplySqAn 使用匿名函数并使用seq而不是rep创建迭代变量。
  • saplyRpCln、laplyRpAn、laplySqAn与laplyRpCln等相同，只是调用sapply而不是lapply/unlist。
  • vaplyRpCln等与laplyRpCln等相同，只是在lapply/unlist的位置使用vapply。
  • replicate 是对replicate的调用，默认为simplify=TRUE。
  • forPre 使用预分配向量的for循环进行填充。
  • forNoPre 使用for循环创建一个空的numeric(0)向量，然后使用c将其连接起来。

这将返回

Unit: microseconds
       expr      min         lq        mean     median         uq      max neval     cld
        all   25.634    31.0705    85.66495    33.6890    35.3400 49240.30  1000 a      
     reduce  542.073   646.7720   780.13592   696.4775   750.2025 51685.44  1000  b     
 laplyRpCln 4349.384  5026.4015  6433.60754  5409.2485  7209.3405 58494.44  1000   c e  
  laplyRpAn 4600.200  5281.6190  6513.58733  5682.0570  7488.0865 55239.04  1000   c e  
  laplySqAn 4616.986  5251.4685  6514.09770  5634.9065  7488.1560 54263.04  1000   c e  
 saplyRpCln 4362.324  5080.3970  6325.66531  5506.5330  7294.6225 59075.02  1000   cd   
  saplyRpAn 4701.140  5386.1350  6781.95655  5786.6905  7587.8525 55429.02  1000     e  
  saplySqAn 4651.682  5342.5390  6551.35939  5735.0610  7525.4725 55148.32  1000   c e  
 vaplyRpCln 4366.322  5046.0625  6270.66501  5482.8565  7208.0680 63756.83  1000   c    
  vaplyRpAn 4657.256  5347.2190  6724.35226  5818.5225  7580.3695 64513.37  1000    de  
  vaplySqAn 4623.897  5325.6230  6475.97938  5769.8130  7541.3895 14614.67  1000   c e  
  replicate 4722.540  5395.1420  6653.90306  5777.3045  7638.8085 59376.89  1000   c e  
     forPre 5911.107  6823.3040  8172.41411  7226.7820  9038.9550 56119.11  1000      f 
   forNoPre 8431.855 10584.6535 11401.64190 10910.0480 11267.5605 58566.27  1000       g

注意，第一种方法明显是最快的。其次是一次性获取全部样本，然后使用Reduce执行求和。在*apply函数中，“干净”的版本，直接使用函数名称似乎有微小的性能提升，而lapply版本似乎与vapply相当，但考虑到值的范围，这个结论并不完全简单明了。sapply似乎是最慢的，尽管函数调用的方法支配着*apply函数的类型。

两个for循环表现最差，预分配for循环优于随c增长的for循环。

在这里，我正在运行一个补丁版本的3.4.1（大约是2017年8月23日左右打的补丁），运行在openSuse 42.1上。

如果您发现任何错误，请告诉我，我会尽快修复。感谢Ben Bolker促使我更深入地研究最终的函数，我在那里发现了一些错误。

为什么forNoPre比forPre快3倍？这似乎非常可疑。 - Ben Bolker
@BenBolker 谢谢，这促使我再次审视。我忽略了存储每个额外抽取的结果 (c(vals, fnc(1L)) 而不是 vals <- c(vals, fnc(1L)) 在 forNoPre 运行中。还有一个更糟糕的错误，我试图对长度为0的向量进行 seq_along 操作... - lmo
4

这个答案只是为了解释你的代码为什么没有起作用。@jake-burkhead 给出了正确的代码写法。

在这段代码中，你只有一半的时间会迈出一步。这是因为你从 xdir 进行采样以决定是否移动。相反，我建议在循环内使用以下代码：

for(i in 2:n){
  x[i] <- x[i - 1] + sample(step, 1)
}
< p > sample(step, 1) 调用决定步行走动的方向是 1 还是 -1。 < /p> < p > 在生成 x 后，您可以使用 cumsum() 计算部分和。结果将是给定点处部分和的向量。 < /p>
谢谢！但是在我的循环中，我应该在哪里包含“i”？ - user3347124
不用谢。我本来想把 i 改成 step，但是不一致。我已经修正了代码。 - Christopher Louden
在你的代码片段中加入 step = c(-1,1) 的定义会很不错。 - logicbloke
0

这是一种实现方式。

GenerateRandomWalk <- function(k = 250,initial.value = 0) {
  # Add a bionomial at each step
  samples = rbinom(k,1,0.5)
  samples[samples==0] = -1
  initial.value + c(0, cumsum(samples))
}
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接