基于python的马尔可夫模型初识
- **1.什么是随机过程?**
- **1.1模拟赌徒的毁灭Gambler's Ruin**
- **2.马尔可夫链(Markov Chains)**
- **2.1马尔可夫链模拟**
- **2.2马尔可夫转移概率图**
- **2.3无记忆性:给定现在,未来独立于过去**
- **2.4 n n n 步转移矩阵**
- 3.示例
import numpy as np
import pandas as pd
import random
import matplotlib.pyplot as plt
#matplotlib parameters
plt.rcParams["figure.figsize"] = (12, 8)
plt.rcParams.update({'font.size': 14})
1.什么是随机过程?
随机过程是随机变量的集合: { X t , t ∈ I } , \{X_t,t \in I \}, {Xt,t∈I},其中: X t X_t Xt是时间 t t t的随机变量集, I I I是过程的索引集。
我们将在本教程中重点介绍的离散时间随机过程是随机变量序列。
1.1模拟赌徒的毁灭Gambler’s Ruin
我们可以用著名的赌徒破产的例子,这是一个随机过程。
在我们的例子中,一个赌徒从$50开始。为了简单起见,他们只能以$1的增量下注。他们每次下注只能赢或输$1。他们会继续赌博,直到输光所有的钱(带着$0离开)或者赢了$100。我们可以用以下符号将其形式化:
$\begin{equation}
X_{g} =
\begin{cases}
$+1, & \text{with a probability of 50%}\
$-1, & \text{with a probability of 50%}\
\end{cases}
\end{equation}$
式中: X g X_g Xg为以美元计的赌博结果.
Source: Introduction to Stochastic Processes with R by Robert P. Dobrow
def gamblers_ruin():gambling_money = 50gambling_goal = 100gambling_simulation = []while gambling_money in range(1,gambling_goal):bet_size = 1w_or_l = random.randrange(-1, 2, step = 2)gambling_money += bet_size * w_or_lgambling_simulation.append(gambling_money)return gambling_simulation
plt.plot(gamblers_ruin())
plt.yticks(np.arange(-20,120,10))
plt.axhline(y=0, color='r', linestyle='-')
plt.axhline(y=100, color='black', linestyle='-')
plt.xlabel('Number of bets')
plt.ylabel('Winnings')
plt.title('Gambling Simulation')
def prob_of_ruin(gambling_goal, initial_gambling_money):return (gambling_goal - initial_gambling_money)/gambling_goalprob_of_ruin(100,50)
sim_list = []while len(sim_list) < 500:sim_list.append(gamblers_ruin()[-1])np.mean(sim_list)
Source: Introduction to Stochastic Processes with R by Robert P. Dobrow
2.马尔可夫链(Markov Chains)
马尔可夫链是一种随机过程。
马尔可夫链是随机变量( X t X_t Xt)的集合,其中未来状态 ( j j j)仅取决于当前状态 ( i i i)。马尔可夫链可以是离散的或连续的。
对于马尔可夫链转移矩阵(表示为 P P P):
- 每一行和为1,即: ∑ j P i j = 1 \sum\limits_{j}P_{ij} = 1 j∑Pij=1.
- 概率必须是非负的,即: P i j ≥ 0 ∀ i , j P_{ij} \geq 0 \ \ \forall \ \ i,j Pij≥0 ∀ i,j
Source: Markov Chain from Wolfram MathWorld
2.1马尔可夫链模拟
马尔可夫链可以由初始分布和转移矩阵模拟。
例子中,初始状态是纽约市New York。从最初的状态,我们可以旅行到:巴黎Paris,开罗Cairo,首尔Seoul,甚至纽约市New York City.
转移矩阵包含从一个状态移到另一个状态的一步转移概率。
mc_example = {'NYC': [.25,0,.75,1],'Paris': [.25,.25,0,0],'Cairo': [.25,.25,.25,0],'Seoul': [.25,.5,0,0]}mc = pd.DataFrame(data = mc_example, index = ['NYC', 'Paris', 'Cairo', 'Seoul'])
2.2马尔可夫转移概率图
我们可以用以下符号将马尔可夫链在初始起点(纽约市)的运动数理化:
P ( X 1 = New York ∣ X 0 = New York ) = P ( X 1 = Paris ∣ X 0 = New York ) = P ( X 1 = Cairo ∣ X 0 = New York ) = P ( X 1 = Seoul ∣ X 0 = New York ) = 25 % P(X_1 = \text{New York}|X_0 = \text{New York})=P(X_1 = \text{Paris}|X_0 = \text{New York})=P(X_1 = \text{Cairo}|X_0 = \text{New York}) = P(X_1 = \text{Seoul}|X_0 = \text{New York}) = 25\% P(X1=New York∣X0=New York)=P(X1=Paris∣X0=New York)=P(X1=Cairo∣X0=New York)=P(X1=Seoul∣X0=New York)=25%
travel_sim = []
travel_sim.append(mc.iloc[0].index[0])
city = np.random.choice(mc.iloc[0].index, p = mc.iloc[0])
travel_sim.append(city)while len(travel_sim) < 25:city = np.random.choice(mc.iloc[mc.index.get_loc(city)].index, p = mc.iloc[mc.index.get_loc(city)])travel_sim.append(city)
travel_sim
['NYC','Cairo','NYC','Seoul','NYC','Paris','Paris','Paris','Seoul','NYC','Cairo','Cairo','NYC','NYC','NYC','Cairo','NYC','NYC','Paris','Cairo','NYC','NYC','Paris','Cairo','Cairo']
mc
NYC Paris Cairo Seoul
NYC 0.25 0.25 0.25 0.25
Paris 0.00 0.25 0.25 0.50
Cairo 0.75 0.00 0.25 0.00
Seoul 1.00 0.00 0.00 0.00
2.3无记忆性:给定现在,未来独立于过去
马尔可夫链的一个重要特征是它们是无记忆的。看例子可知,唯一重要的状态是当前状态。 如果我们的从纽约开始 ( X 0 = N Y C X_0 = NYC X0=NYC) 去了巴黎 ( X 1 = P a r i s X_1 = Paris X1=Paris), 然后到下一个城市 ( X 2 X_2 X2)只取决于从巴黎出发的可能性. 最初在纽约开始出发这一事实并不影响到 X 2 . X_2. X2.
2.4 n n n 步转移矩阵
mc.to_numpy()
array([[0.25, 0.25, 0.25, 0.25],[0. , 0.25, 0.25, 0.5 ],[0.75, 0. , 0.25, 0. ],[1. , 0. , 0. , 0. ]])
def matrix_power(matrix, power):if power == 0:return np.identity(len(matrix))elif power == 1:return matrixelse:return np.dot(matrix, matrix_power(matrix, power-1))
matrix_power(mc.to_numpy(), 2)
array([[0.5 , 0.125 , 0.1875, 0.1875],[0.6875, 0.0625, 0.125 , 0.125 ],[0.375 , 0.1875, 0.25 , 0.1875],[0.25 , 0.25 , 0.25 , 0.25 ]])
np.dot(np.dot(mc.to_numpy(), mc.to_numpy()))
array([[0.5 , 0.125 , 0.1875, 0.1875],[0.6875, 0.0625, 0.125 , 0.125 ],[0.375 , 0.1875, 0.25 , 0.1875],[0.25 , 0.25 , 0.25 , 0.25 ]])
matrix_power(mc.to_numpy(), 2).sum(axis=1)
array([1., 1., 1., 1.])
上述两步转移矩阵告诉我们的是,如果我们观察纽约的状态 i i i,那么在两步之后,有50%的概率回到纽约市,有12.5%的概率会前往巴黎,有18.75%的概率会到达开罗或首尔。
for i in range(1,15,1):print(f'n Step Transition Matrix at the nth power {i}\n', matrix_power(mc.to_numpy(), i),'\n')
n Step Transition Matrix at the nth power 1[[0.25 0.25 0.25 0.25][0. 0.25 0.25 0.5 ][0.75 0. 0.25 0. ][1. 0. 0. 0. ]] n Step Transition Matrix at the nth power 2[[0.5 0.125 0.1875 0.1875][0.6875 0.0625 0.125 0.125 ][0.375 0.1875 0.25 0.1875][0.25 0.25 0.25 0.25 ]] n Step Transition Matrix at the nth power 3[[0.453125 0.15625 0.203125 0.1875 ][0.390625 0.1875 0.21875 0.203125][0.46875 0.140625 0.203125 0.1875 ][0.5 0.125 0.1875 0.1875 ]] n Step Transition Matrix at the nth power 4[[0.453125 0.15234375 0.203125 0.19140625][0.46484375 0.14453125 0.19921875 0.19140625][0.45703125 0.15234375 0.203125 0.1875 ][0.453125 0.15625 0.203125 0.1875 ]] n Step Transition Matrix at the nth power 5[[0.45703125 0.15136719 0.20214844 0.18945312][0.45703125 0.15234375 0.20214844 0.18847656][0.45410156 0.15234375 0.203125 0.19042969][0.453125 0.15234375 0.203125 0.19140625]] n Step Transition Matrix at the nth power 6[[0.45532227 0.15209961 0.20263672 0.18994141][0.4543457 0.15234375 0.20288086 0.19042969][0.45629883 0.15161133 0.20239258 0.18969727][0.45703125 0.15136719 0.20214844 0.18945312]] n Step Transition Matrix at the nth power 7[[0.45574951 0.15185547 0.20251465 0.18988037][0.45617676 0.15167236 0.20239258 0.1897583 ][0.45556641 0.15197754 0.20257568 0.18988037][0.45532227 0.15209961 0.20263672 0.18994141]] n Step Transition Matrix at the nth power 8[[0.45570374 0.15190125 0.20252991 0.18986511][0.45559692 0.15196228 0.20256042 0.18988037][0.45570374 0.15188599 0.20252991 0.18988037][0.45574951 0.15185547 0.20251465 0.18988037]] n Step Transition Matrix at the nth power 9[[0.45568848 0.15190125 0.20253372 0.18987656][0.45569992 0.1518898 0.20252991 0.18988037][0.45570374 0.15189743 0.20252991 0.18986893][0.45570374 0.15190125 0.20252991 0.18986511]] n Step Transition Matrix at the nth power 10[[0.45569897 0.15189743 0.20253086 0.18987274][0.45570278 0.15189743 0.20252991 0.18986988][0.45569229 0.15190029 0.20253277 0.18987465][0.45568848 0.15190125 0.20253372 0.18987656]] n Step Transition Matrix at the nth power 11[[0.45569563 0.1518991 0.20253181 0.18987346][0.45569301 0.15190005 0.20253253 0.18987441][0.4556973 0.15189815 0.20253134 0.18987322][0.45569897 0.15189743 0.20253086 0.18987274]] n Step Transition Matrix at the nth power 12[[0.45569623 0.15189868 0.20253164 0.18987346][0.45569706 0.15189826 0.2025314 0.18987328][0.45569605 0.15189886 0.2025317 0.1898734 ][0.45569563 0.1518991 0.20253181 0.18987346]] n Step Transition Matrix at the nth power 13[[0.45569624 0.15189873 0.20253164 0.1898734 ][0.45569609 0.15189883 0.20253168 0.1898734 ][0.45569618 0.15189873 0.20253165 0.18987344][0.45569623 0.15189868 0.20253164 0.18987346]] n Step Transition Matrix at the nth power 14[[0.45569618 0.15189874 0.20253165 0.18987342][0.45569618 0.15189873 0.20253165 0.18987344][0.45569623 0.15189873 0.20253164 0.18987341][0.45569624 0.15189873 0.20253164 0.1898734 ]]
我们可以看到,随着步数的增加,概率会趋于一种被称为稳态分布(也称为稳态向量)的状态。
假设这次从首尔开始我们的旅程。在这种情况下,我们可以将初始分布表示为: α = [ 0 , 0 , 0 , 1 ] \alpha = [0,0,0,1] α=[0,0,0,1]。现在我们来找出从现在起两次行程后最终回到首尔的概率。
initial_dist = np.asarray([0,0,0,1])mc_p2 = matrix_power(mc.to_numpy(),2)np.dot(initial_dist,mc_p2)
array([0.25, 0.25, 0.25, 0.25])
3.示例
mc_example2 = {'hamburger': [.2,.3,.5],'pizza': [.6,0,0],'hot-dog': [.2,.7,.5]}mc2 = pd.DataFrame(data = mc_example2, index = ['hamburger', 'pizza', 'hot-dog'])
mc2
hamburger pizza hot-dog
hamburger 0.2 0.6 0.2
pizza 0.3 0.0 0.7
hot-dog 0.5 0.0 0.5
np.dot(mc2.to_numpy(), mc2.to_numpy())
array([[0.32, 0.12, 0.56],[0.41, 0.18, 0.41],[0.35, 0.3 , 0.35]])
steps = 2
def matrix_power(matrix, power):if power == 0:return np.identity(len(matrix))elif power == 1:return matrixelse:return np.dot(matrix, matrix_power(matrix, power-1))
matrix_power(mc2.to_numpy(), steps)
array([[0.32, 0.12, 0.56],[0.41, 0.18, 0.41],[0.35, 0.3 , 0.35]])
matrix_power(mc2.to_numpy(), steps).sum(axis=1)
array([1., 1., 1.])
for i in range(1,11,1):print(f'n Step Transition Matrix at the nth step {i}\n', matrix_power(mc2.to_numpy(), i),'\n')
n Step Transition Matrix at the nth step 1[[0.2 0.6 0.2][0.3 0. 0.7][0.5 0. 0.5]] n Step Transition Matrix at the nth step 2[[0.32 0.12 0.56][0.41 0.18 0.41][0.35 0.3 0.35]] n Step Transition Matrix at the nth step 3[[0.38 0.192 0.428][0.341 0.246 0.413][0.335 0.21 0.455]] n Step Transition Matrix at the nth step 4[[0.3476 0.228 0.4244][0.3485 0.2046 0.4469][0.3575 0.201 0.4415]] n Step Transition Matrix at the nth step 5[[0.35012 0.20856 0.44132][0.35453 0.2091 0.43637][0.35255 0.2145 0.43295]] n Step Transition Matrix at the nth step 6[[0.353252 0.210072 0.436676][0.351821 0.212718 0.435461][0.351335 0.21153 0.437135]] n Step Transition Matrix at the nth step 7[[0.35201 0.2119512 0.4360388][0.3519101 0.2110926 0.4369973][0.3522935 0.210801 0.4369055]] n Step Transition Matrix at the nth step 8[[0.35200676 0.211206 0.43678724][0.35220845 0.21114606 0.43664549][0.35215175 0.2113761 0.43647215]] n Step Transition Matrix at the nth step 9[[0.35215677 0.21120406 0.43663917][0.35210825 0.21132507 0.43656668][0.35207925 0.21129105 0.43662969]] n Step Transition Matrix at the nth step 10[[0.35211216 0.21129406 0.43659378][0.35210251 0.21126495 0.43663254][0.35211801 0.21124755 0.43663443]]
initial_dist = np.asarray([0,1,0])for i in range(1,11,1):print(f'start pizza then n Step Transition Vector at the nth step {i}\n', np.dot(initial_dist,matrix_power(mc2.to_numpy(), i)),'\n')
start pizza then n Step Transition Vector at the nth step 1[0.3 0. 0.7] start pizza then n Step Transition Vector at the nth step 2[0.41 0.18 0.41] start pizza then n Step Transition Vector at the nth step 3[0.341 0.246 0.413] start pizza then n Step Transition Vector at the nth step 4[0.3485 0.2046 0.4469] start pizza then n Step Transition Vector at the nth step 5[0.35453 0.2091 0.43637] start pizza then n Step Transition Vector at the nth step 6[0.351821 0.212718 0.435461] start pizza then n Step Transition Vector at the nth step 7[0.3519101 0.2110926 0.4369973] start pizza then n Step Transition Vector at the nth step 8[0.35220845 0.21114606 0.43664549] start pizza then n Step Transition Vector at the nth step 9[0.35210825 0.21132507 0.43656668] start pizza then n Step Transition Vector at the nth step 10[0.35210251 0.21126495 0.43663254]