IT评测·应用市场-qidao123.com技术社区

标题: 基于python的马尔可夫模型初识 [打印本页]

作者: 莫张周刘王 时间: 2024-10-28 04:53
标题: 基于python的马尔可夫模型初识

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 = 50
gambling_goal = 100
gambling_simulation = []
while gambling_money in range(1,gambling_goal):
bet_size = 1
w_or_l = random.randrange(-1, 2, step = 2)
gambling_money += bet_size * w_or_l
gambling_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_goal
prob_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']

复制代码

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 matrix
else:
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]}mc
2 = pd.DataFrame(data = mc
_example2, index = ['hamburger', 'pizza', 'hot-dog'])mc
2

复制代码

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(mc
2.to_numpy(), mc
2.to_numpy())

复制代码

array([[0.32, 0.12, 0.56],
[0.41, 0.18, 0.41],
[0.35, 0.3 , 0.35]])

复制代码

steps = 2def matrix_power(matrix, power):
if power == 0:
return np.identity(len(matrix))
elif power == 1:
return matrix
else:
return np.dot(matrix, matrix_power(matrix, power-1))
matrix_power(mc
2.to_numpy(), steps)

复制代码

array([[0.32, 0.12, 0.56],
[0.41, 0.18, 0.41],
[0.35, 0.3 , 0.35]])

复制代码

matrix_power(mc
2.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(mc
2.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(mc
2.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]

复制代码

免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。

欢迎光临 IT评测·应用市场-qidao123.com技术社区 (https://dis.qidao123.com/)