使用Python的馬爾科夫鏈實例的實踐 一行行的代碼解密馬爾可夫鏈。

當我開始學習物理時，我并不喜歡概率的概念。我對用物理學可以對整個世界進行建模的想法非常振奮，不確定性的想法讓我很生氣:)

事實是，當我們想研究真實的現象時，我們遲早都必須處理一定程度的不確定性。而處理它的唯一方法是獲得對支配我們過程的概率的準確估計。

馬爾科夫鏈是一個很好的方法。馬爾科夫鏈背后的想法是非常簡單的。

未來將發生的一切只取決于現在正在發生的事情。

用數學術語來說，我們說有一個隨機變量X_0, X_1, ..., X_n的序列，可以在某個集合A中取值。然后我們說，如果一個事件的序列是一個馬爾科夫鏈，我們就有。

圖片由我用LaTeX生成 如何用python生成LaTex請參考鏈接：

這聽起來很復雜，但它只不過是上面所表達的概念而已。

另一個假設是，該方程對每一步都有效（不僅是最后一步），而且概率總是相同的（盡管從形式上看，這只對同質馬爾科夫鏈而言是真的）。

現在，可能狀態A的集合通常被表示為樣本空間S，你可以用所謂的過渡概率來描述從S中的一個狀態x到S中的一個狀態y的概率。

但我答應過你，這將是一篇 "手把手 "的文章，所以讓我們開始把這些概念形象化吧!

公平地說，Python不是進行數值模擬的最佳環境。專業研究人員使用的是更復雜的、在某些方面更可靠的語言，如C或Fortran。

盡管如此，本博客的目標是介紹一些非常簡單的概念，使用Python可以使這個學習過程更容易。

import?matplotlib.pyplot?as?plt
import?numpy?as?np
import?pandas?as?pd
import?seaborn?as?sns
plt.style.use('ggplot')
plt.rcParams['font.family']?=?'sans-serif'?
plt.rcParams['font.serif']?=?'Ubuntu'?
plt.rcParams['font.monospace']?=?'Ubuntu?Mono'?
plt.rcParams['font.size']?=?14?
plt.rcParams['axes.labelsize']?=?12?
plt.rcParams['axes.labelweight']?=?'bold'?
plt.rcParams['axes.titlesize']?=?12?
plt.rcParams['xtick.labelsize']?=?12?
plt.rcParams['ytick.labelsize']?=?12?
plt.rcParams['legend.fontsize']?=?12?
plt.rcParams['figure.titlesize']?=?12?
plt.rcParams['image.cmap']?=?'jet'?
plt.rcParams['image.interpolation']?=?'none'?
plt.rcParams['figure.figsize']?=?(12,?10)?
plt.rcParams['axes.grid']=False
plt.rcParams['lines.linewidth']?=?2?
plt.rcParams['lines.markersize']?=?8
colors?=?['xkcd:pale?orange',?'xkcd:sea?blue',?'xkcd:pale?red',?'xkcd:sage?green',?'xkcd:terra?cotta',?'xkcd:dull?purple',?'xkcd:teal',?'xkcd:?goldenrod',?'xkcd:cadet?blue',
'xkcd:scarlet']

因此，讓我們深入了解一下：這是你需要的東西是一堆主流的庫，如pandas、matplotlib、seaborn和numpy。

讓我們從最簡單的場景開始。

隨機漫步

簡單隨機漫步是一個極其簡單的隨機漫步的例子。

第一個狀態是0，然后你以0.5的概率從0跳到1，以0.5的概率從0跳到-1。

圖片由我使用Power Point制作

然后你對x_1, x_2, ..., x_n做同樣的事情。

你認為S_n是時間n的狀態。

有可能證明（實際上非常容易），在時間t+1時處于某種狀態的概率，即一個整數x，只取決于時間t的狀態。

因此，這就是如何生成它。

start?=?0
x?=?[]
n?=?10000
for?i?in?range(n):
????step?=?np.random.choice([-1,1],p=[0.5,0.5])
????start?=?start?+?step
????x.append(start)
????
plt.plot(x)
plt.xlabel('Steps',fontsize=20)
plt.ylabel(r'$S_{n}$',fontsize=20)

而這就是結果：

現在，隨機漫步的想法是模擬如果我們決定從一個點開始，通過投擲一枚完美的硬幣隨機選擇向上或向下，將會發生什么。

這個過程相當簡單，但就其理論應用和特性而言，卻非常有趣。

這個過程的第一個合理的擴展是考慮一個隨機行走，但有一個非完美的硬幣。這意味著，上升的概率與下降的概率不一樣。這被稱為有偏見的隨機漫步。

讓我們考慮以下幾個概率。

?[0.1,0.9]?,?[0.2,0.8],?[0.4,0.6],?
?[0.6,0.4],?[0.8,0.2],[0.9,0.1]

所以我們有6種可能的隨機行走。注意，概率必須是1，因此考慮 "向上 "或 "向下 "的概率即可。

這里是你如何做的。

x?=?[]
p?=?[[0.5,0.5],[0.9,0.1],[0.8,0.2],[0.6,0.4],[0.4,0.6],[0.2,0.8],[0.1,0.9]]
label_p?=?['Simple',r'$p=0.9$',r'$p=0.8$',r'$p=0.6$',r'$p=0.4$',r'$p=0.2$',r'$p=0.1$']
n?=?10000
x?=?[]
for?couple?in?p:
????x_p?=?[]
????start?=?0
????for?i?in?range(n):
????????step?=?np.random.choice([-1,1],p=couple)
????????start?=?start?+?step
????????x_p.append(start)
????x.append(x_p)

而這是我們將其形象化后的結果。

i=0
for?time_series?in?x:
????plt.plot(time_series,?label?=?label_p[i])
????i=i+1
plt.xlabel('Steps',fontsize=20)
plt.ylabel(r'$S_{n}$',fontsize=20)
plt.legend()

賭徒的毀滅鏈

另一種擴展隨機漫步的簡單方法是賭徒的毀滅鏈。 從概念上講，它與隨機漫步非常相似：你從一個狀態x開始，你可以以概率p進入一個狀態y=x+1，或者以概率1-p進入一個狀態y=x-1。

有趣的是，當你到達1或N時，你基本上就被卡住了。你只能永遠停留在這個狀態下，別無他法。

這個函數給定：

起始點（例如3）

第一個可能的值（例如：0）

和最后一個可能的值（例如：5）

步驟數（如10000）

給你最終的狀態。

def?gamblersruinchain(start,first,last,n):
????for?k?in?range(n):
????????if?start==first?or?start==last:
????????????start?=?start
????????else:
????????????step?=?np.random.choice([-1,1],p=[0.5,0.5])
????????????start?=?start?+?step
????return?start

現在，在嘗試這個函數之前，讓我們考慮一個更有趣的情況。

假設我們從狀態3開始。兩步之后，結束在狀態5的概率是多少？

嗯，就是從狀態3到狀態4，然后再從狀態4到狀態5的概率。

LaTeX制作

在我們的例子中，它只是0.25。

如果現在我們問這個方程。

假設我們從狀態3開始。兩步之后，結束在狀態1的概率是多少？

同樣，這是從狀態3到狀態2，再從狀態2到狀態1的概率。

唯一的其他選擇是在兩步之后從狀態3到狀態3。我們可以用一個非常簡單的方法來計算這個概率。由于總的概率必須是1，所以它只是。

圖片由我用LaTeX制作 而如果p=0.5，則又是0.5。

同樣，概率的概念是，如果我們無限次地重復一個實驗，我們應該能夠驗證概率值所暗示的發生情況。

state_list?=?[]
for?i?in?range(100000):
????state_list.append(gamblersruinchain(3,0,5,2))
data_state?=?pd.DataFrame({'Final?State':state_list})
data_occ?=?pd.DataFrame(data_state.value_counts('Final?State')).rename(columns={0:'Count'})
data_occ['Count']?=?data_occ['Count']/100000
sns.barplot(x=data_occ.index,y=data_occ['Count'],palette='plasma')
plt.ylabel('Normalized?Count')

自定義馬爾科夫鏈

前面的模型是眾所周知的，并被用作馬爾科夫鏈的介紹性例子。讓我們試著發揮創意，建立一個全新的、不存在的模型，就像下圖中的模型。

圖片由我使用Power Point制作 我是個糟糕的畫師，但這個模型本身很簡單。

當你看到兩個節點（比方說A和B）之間有一個箭頭時，這意味著你可以從節點A出發，以一定的概率去到節點B，這個概率是用黑色書寫的。

例如，從狀態A到狀態B的概率是0.5。

一個重要的概念是，模型可以用過渡矩陣來概括，它解釋了馬爾科夫鏈中可能發生的一切。這就是我們模型的過渡矩陣。

state_1?=?[0.2,0.5,0.3,0,0]
state_2?=?[0,0.5,0.5,0,0]
state_3?=?[0,0,1,0,0]
state_4?=?[0,0,0,0,1]
state_5?=?[0,0,0,0.5,0.5]
trans_matrix?=?[state_1,state_2,state_3,state_4,state_5]
trans_matrix?=?np.array(trans_matrix)
trans_matrix
array([[0.2,?0.5,?0.3,?0.?,?0.?],
???????[0.?,?0.5,?0.5,?0.?,?0.?],
???????[0.?,?0.?,?1.?,?0.?,?0.?],
???????[0.?,?0.?,?0.?,?0.?,?1.?],
???????[0.?,?0.?,?0.?,?0.5,?0.5]])

如果你仔細觀察這個模型，你可以看到一些非常特別的東西。比方說，你從狀態2跳到狀態3。你能回到狀態2嗎？答案是不能。

同樣的情況也適用于狀態3和狀態1。因此，狀態1、3和2被定義為短暫的狀態。

另一方面，如果你從狀態4開始，總是有可能在某一時刻，你會回到狀態4。同樣的情況也適用于狀態5。這些狀態被稱為反復出現的狀態。

讓我們做一些實驗，以便我們能夠正確理解這個概念。

直觀地講，我們可以看到，從狀態2開始，不回到狀態2的概率隨著步驟數的增加而趨于0。

事實上，從狀態2開始，我們在N步之后發現自己處于狀態2的概率是如下。

圖片由我用LaTeX制作 事實上，如果我們從狀態2到狀態3，就不可能再回到狀態2。讓我們把這個理論函數定義為t(N)，并把它畫出來。

def?t(N):
????step?=?np.arange(1,N+1,1)
????y?=?[]
????for?s?in?step:
????????v?=?0.5**s
????????y.append(v)
????return?y
????
plt.plot(t(10))
plt.ylabel(r'$t(N-1)$',fontsize=20)
plt.xlabel(r'$N-1$',fontsize=20)

現在，讓我們使用馬爾科夫鏈，看看我們是否驗證了同樣的結果。

我們從狀態2開始，在N步之后驗證處于狀態2的概率。在這種情況下，概率只是最終狀態中2的數量與發生次數的比率。為了保持一致，出現的次數需要趨于無窮大。讓我們考慮1000次測試。

這就是我們要使用的函數。

def?prob(N):
????states?=?np.arange(1,6,1)
????steps?=?np.arange(1,N+1,1)
????n=1000
????state_collection?=?[]
????for?k?in?range(n):
????????start?=?2?
????????for?i?in?range(N):
????????????start?=?np.random.choice(states,p=trans_matrix[start-1])
????????if?start==2:
????????????state_collection.append(1)
????????else:
????????????state_collection.append(0)
????state_collection?=?np.array(state_collection)
????return?state_collection.sum()/n

讓我們對各種N使用這個函數，并稱其為p(N)。

def?p(N):
????step?=?np.arange(1,N+1,1)
????y?=?[]
????for?s?in?step:
????????v?=?prob(s)
????????y.append(v)
????return?y
????
p_20?=?p(20)
plt.plot(t(20),label=r'Theory,?$t(N-1)$')
plt.plot(p_20,'x',label=r'Simulation,?$p(N-1)$',color='navy')
plt.ylabel(r'Result',fontsize=20)
plt.xlabel(r'$N-1$',fontsize=20)
plt.legend()

可以看到，我們使用了過渡矩陣來做這個模擬。我們可以使用過渡矩陣來評估我們所考慮的馬爾科夫鏈的所有屬性。

結論

在這本筆記本中，我們已經看到了非常知名的模型，如隨機漫步和賭徒毀滅鏈。然后，我們創建了自己的全新模型，并對其進行了一番研究，發現了一些重要的概念，如過渡矩陣、循環狀態和瞬態。最重要的是，我們已經看到了如何使用Python和非常知名的庫以非常簡單的方式驗證這些概念。

本文由 mdnice 多平臺發布

總結

以上是生活随笔為你收集整理的一行行的代码解密马尔可夫链的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。