醉汉漫步和二维随机行走路径（图片转自http://blog.sciencenet.cn/blog-677221-1071588.html）
假设张三在家喝醉酒后，然后出来溜达，每一秒钟会随机朝东南西北中的一个方向走一步，如此继续下去。。。那么有人就会问：张三走了1000步后会在哪呢？5000步呢？是不是随着时间增长，张三离家越来越远呢？
上述问题被称为“醉汉漫步”，数学家把它抽象成一个模型，叫随机行走（random walk）。由于醉汉只能在地面上游荡，因此是一种二维随机行走。
这是一个典型的马尔可夫过程。因为醉汉在t+1时刻的状态（即位置），仅仅由他在t时刻的状态以及他随机选择的方向所决定，与过去（t之前）所走过的路径无关。
这个过程看似非常随机，很难解决。但是如果用电脑程序来模拟，通过绘图很容易把醉汉的行走轨迹以及醉汉离家的距离展现出来。
下面我们用R来模拟这个过程：
通过最开始的140步来细看这个过程，看下面的动画：图中黑点为原点（即醉汉的家）；蓝点和蓝线是醉汉经过的线路；红点是醉汉当前时刻的位置，图下方用红字显示了当前位置的坐标。
经过观察发现，醉汉在最开始的140步中有好几次回到了原点（家）。
绘制上面动画的R代码如下：
############ 编写二维随机行走函数
random_walk_2d <- function(x0=0,
y0=0,
N=100,
random_seed=1234){
##### 参数意义:
# x0：醉汉初始x坐标；
# y0：醉汉初始y坐标；
# N：醉汉行走的步数；
# random_seed：随机数种子，改变它可产生不一样的随机序列。
##### 返回值:
# 一个数据框(3列)，包含每一步的序号(0, 1, 2, 3......N)以及每一步的x和y坐标。
set.seed(random_seed) ##设置随机数种子
rand <- sample.int(4, N, replace=TRUE) ##产生可重复的随机序列，例如：2, 3, 1, 2, 4......
x <- x_temp <- x0
y <- y_temp <- y0
for(i in 1:N){
if(rand[i]==1L) ##如果为1, 向前走一步，x坐标加1
x_temp <- x_temp+1L
else if(rand[i]==2L) ##如果为2, 向后走一步，x坐标减1
x_temp <- x_temp-1L
else if(rand[i]==3L) ##如果为3, 向上走一步，y坐标加1
y_temp <- y_temp+1L
else ##如果为4, 向下走一步，y坐标减1
y_temp <- y_temp-1L
x <- c(x, x_temp)
y <- c(y, y_temp)
}
result <- data.frame(step=0:N, x=x, y=y) ##结果数据为数据框
return(result) ##返回结果
}
############ 绘制最开始140步的动画
library(animation) ##加载绘制动画所需的包
result <- random_walk_2d(N=50000, random_seed=12)
n <- 140
result_plot <- result[1:n, ] ##只选取前140步
x <- result_plot$x
y <- result_plot$y
max_xy <- max(abs(range(c(x,y)))) ##计算x和y坐标的最大范围
ani.options(interval=0.3, ani.width=800, ani.height=800, autobrowse=FALSE) ##设置动画参数
saveGIF(
for(i in 1:n)
{
plot(x[1:i], y[1:i], type='o', asp=1, xaxs="i",yaxs="i",
xlab='x',ylab='y',cex.lab=2.5, cex.axis=1.5,cex.main=3,
xlim=c(-max_xy, max_xy), ylim=c(-max_xy, max_xy),
col='blue', pch=1, cex=2, lwd=2.5,
main=paste0('N = ', i-1L))
points(0, 0, pch=19, cex=2) ##显示坐标原点
points(x[i], y[i], pch=19, col='red', cex=2) ##设置当前点为红实点
text(0, -max_xy, paste0('(', x[i], ', ', y[i], ')'), pos=3, col='red', cex=2.5) ##显示当前位置的坐标
abline(h=seq(-max_xy, max_xy), col="gray50", lty=3) ##绘制水平线
abline(v=seq(-max_xy, max_xy), col="gray50", lty=3) ##绘制竖直线
}, movie.name='2d_random_walk_1.gif')
对于“醉汉回家”这个问题，美籍匈牙利数学家波利亚做了详细研究。
波利亚
对于随机行走，波利亚1921年得到的结论总结如下：
-
在一维空间中，一个随机选择下一步往哪儿走的醉汉有无穷多的机会回到家中
-
在二维网格道路上，醉汉在经过足够长的时间后也一定能够回到家中
在波利亚的研究基础之上，后来的数学家们对多维情况进行了探索，结论如下图：
多维情况下，回到原点的概率
在三维的道路网格上，醉汉即使经过足够长的时间后也只有34%的概率能够回到他自己的家中；四维情况只有大约19%的概率。。。
以上结论可参考网页：http://mathworld.wolfram.com/PolyasRandomWalkConstants.html
对于上面的结论，有人调侃：弄丢的狗儿（二维运动）一定能找回家的，而撒出去的鸟儿（三维运动）很可能就回不来了。有人可能会反驳，为啥我家弄丢的狗没回家啊？那可能是你家的狗活的时间太短，或者狗的运动压根不是随机的。。。哈哈哈
下面来看一个时间更长（5万步）的随机行走模拟：
可以看到，醉汉经过的轨迹非常随机。
绘制上面动画的代码如下：
n <- 50001
result_plot <- result[1:n, ]
x <- result_plot$x
y <- result_plot$y
max_xy <- max(abs(range(c(x,y))))
ani.options(interval=0.3, ani.width=800, ani.height=800, autobrowse=FALSE)
saveGIF(
for(i in seq.int(1, n, 500)) ##每隔500步绘制一张
{
plot(x[1:i], y[1:i], type='l', asp=1, xaxs="i",yaxs="i",
xlab='x',ylab='y',cex.lab=2.5, cex.axis=1.5,cex.main=3,
xlim=c(-max_xy, max_xy), ylim=c(-max_xy, max_xy),
col='blue', pch=1, cex=2,lwd=1.2,
main=paste0('N = ', i-1L))
points(0, 0, pch=19, cex=2)
points(x[i], y[i], pch=19, col='red', cex=2)
}, movie.name='2d_random_walk_2.gif')
本博客中的核心代码是最上面的random_walk_2d()函数，通过改变随机数种子（random_seed参数）可以模拟出大量不重复的二维随机行走过程。然后根据这些模拟结果可以去检验已有相关理论是否正确？这里就不介绍了，感兴趣的可以去研究研究。
如今，随机行走模型是一种解决随机问题的重要方法，它与人类生活息息相关，例如醉汉行走的轨迹、扩散过程、股票的涨跌等都可以用它来模拟。它已经被广泛应用到数学、物理、生物、医学以及经济学等领域。
感谢您的阅读！想了解更多有关R语言技巧，请关注我的微信公众号“R语言和Python学堂”，我将定期更新相关文章。
网友评论