醉漢漫步和二維隨機行走路徑（圖片轉自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%的概率。。。

以上結論可參考網頁：mathworld.wolfram.com/P

對於上面的結論，有人調侃：弄丟的狗兒（二維運動）一定能找回家的，而撒出去的鳥兒（三維運動）很可能就回不來了。有人可能會反駁，為啥我家弄丟的狗沒回家啊？那可能是你家的狗活的時間太短，或者狗的運動壓根不是隨機的。。。哈哈哈

下面來看一個時間更長（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學堂」，我將定期更新相關文章。


推薦閱讀：
相关文章