r - R中的征费步行模拟

Question

我正在尝试生成一系列数字来模拟 R 中的征费步行。目前我正在使用以下代码：

alpha=2
n=1000
x=rep(0,n)
y=rep(0,n)

for (i in 2:n){
   theta=runif(1)*2*pi
   f=runif(1)^(-1/alpha)
   x[i]=x[i-1]+f*cos(theta)
   y[i]=y[i-1]+f*sin(theta)
}

代码按预期工作，我能够根据我的要求生成数字。下图显示了这样的 Levy Walk：

以下直方图证实了生成的数字（即 f）实际上属于幂律：

我的问题如下：生成的步长（即 f）非常大。我可以修改代码以使步长仅在某个范围内 [fmin, fmax] 吗？

PS我故意没有对代码进行矢量化。

Answer 1

尝试使用这个：

f=runif(1, fmax^(-alpha), fmin^(-alpha))^(-1/alpha)

请注意，您需要0 < fmin < fmax.

顺便说一句，您可以像这样对代码进行矢量化：

theta <- runif(n-1)*2*pi
f <- runif(n-1, fmax^(-alpha), fmin^(-alpha))^(-1/alpha)
x <- c(0, cumsum(f*cos(theta)))
y <- c(0, cumsum(f*sin(theta)))

Answer 2

为了精确起见，您在这里模拟的是 Lévy 飞行。要使其成为 Lévy walk，您应该允许粒子从每次飞行的开始到结束“行走”（for例如，使用 ）。如果你用你的代码绘制你的模拟结果，plot(x, y, type = "o")你会看到飞行中没有位置（没有步行）。

Answer 3

library(ggplot2)
library(gridExtra)

alpha= 5
n= 1000
x= rep(0,n)
y= rep(0,n)

fmin= 1
fmax= n

for (i in 2:n){
   theta= runif(n-1)*2*pi
   f= runif(n-1, fmax^(-alpha), fmin^(-alpha))^(-1/alpha)
   x= c(0, cumsum(f*cos(theta)))
   y= c(0, cumsum(f*sin(theta)))
}
ggplot(data.frame(x=x, y=y), aes(x, y))+geom_point()+geom_path()