mirror of
https://github.com/microsoft/ai-edu.git
synced 2026-09-21 04:26:39 +08:00
@@ -32,16 +32,23 @@ def try_once(n_doors: int):
|
||||
|
||||
return no_change_but_win, win_after_change
|
||||
|
||||
if __name__ == "__main__":
|
||||
|
||||
def try_n_doors(n_doors):
|
||||
total = 100000
|
||||
n_doors = 8
|
||||
n_win_0 = 0
|
||||
n_win_1 = 0
|
||||
for i in range(total):
|
||||
no_change_but_win, win_after_change = try_once(n_doors)
|
||||
n_win_0 += win_after_change
|
||||
n_win_1 += no_change_but_win
|
||||
|
||||
print(str.format("{0}扇门:", n_doors))
|
||||
print(n_win_0, n_win_1)
|
||||
print(str.format("更换选择而中奖的概率={0} \n不更换而中奖的概率={1}",
|
||||
print(str.format("更换选择而中奖的概率={0} \n不换选择而中奖的概率={1}",
|
||||
n_win_0/total, n_win_1/total))
|
||||
|
||||
if __name__ == "__main__":
|
||||
n_doors = 3
|
||||
try_n_doors(n_doors)
|
||||
|
||||
n_doors = 8
|
||||
try_n_doors(n_doors)
|
||||
|
||||
Binary file not shown.
|
Before Width: | Height: | Size: 64 KiB |
|
Before Width: | Height: | Size: 12 KiB After Width: | Height: | Size: 12 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 66 KiB |
@@ -6,7 +6,7 @@
|
||||
|
||||
参赛者面前有三扇关闭着的门,其中一扇的后面是一辆汽车,选中后面有车的那扇门就可以赢得该汽车,而另外两扇门后面则各藏有一只山羊。当参赛者选定了一扇门,但未去开启它的时候,主持人会开启剩下两扇门中的一扇,露出其中一只山羊。主持人其后会问参赛者要不要换另一扇仍然关上的门。问题是:换另一扇门是否会增加参赛者赢得汽车的机率?
|
||||
|
||||
<img src="./images/ThreeDoors1.png" width="500">
|
||||
<img src="./img/ThreeDoors1.png" width="500">
|
||||
|
||||
图1
|
||||
|
||||
@@ -102,7 +102,7 @@
|
||||
|
||||
另外,有些读者看表格可能有困难,所以我们把这些情况变成概率数字放在图 2 中。
|
||||
|
||||
<img src="./images/ThreeDoors2.png" width="600">
|
||||
<img src="./img/ThreeDoors2.png" width="600">
|
||||
|
||||
图 2
|
||||
|
||||
@@ -263,7 +263,7 @@ if __name__ == "__main__":
|
||||
```
|
||||
66259 33741
|
||||
更换选择而中奖的概率=0.66259
|
||||
不更换而中奖的概率=0.33741
|
||||
不换选择而中奖的概率=0.33741
|
||||
```
|
||||
|
||||
更换选择中奖的概率约等于 $\frac{2}{3}$,不更换选择中奖的概率约等于 $\frac{1}{3}$,与穷举法和理论推导结论一致。
|
||||
@@ -273,7 +273,7 @@ if __name__ == "__main__":
|
||||
```
|
||||
14554 12480
|
||||
更换选择而中奖的概率=0.14554
|
||||
不更换而中奖的概率=0.1248
|
||||
不换选择而中奖的概率=0.1248
|
||||
```
|
||||
|
||||
比如当 n_doors=8 时,按公式 2,结果是 $0.1248 \approx \frac{1}{8}$;按公式 3,结果是 $0.14554 \approx \frac{8-1}{8 \times (8-2)}=\frac{7}{48}$。
|
||||
@@ -288,4 +288,4 @@ if __name__ == "__main__":
|
||||
1. 强化学习中,基于模型的理论部分,也是以概率论为基础的,可以先热复习一下热身。
|
||||
2. 图 2 中,其实是由“策略-动作-状态”组成的,这与强化学习的理论基本一致。
|
||||
3. 使用代码模拟,也是一种“聪明的笨办法”,利用计算机快速模拟实际环境的交互,这也是强化学习的重要方法。
|
||||
4. 用代码理解公式及理论知识,是程序员的一种重要技能。
|
||||
4. 对于没有受过训练的可以直接阅读并理解公式推导的读者来说,用代码理解公式及理论知识,是一种重要手段。
|
||||
|
||||
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
@@ -0,0 +1,51 @@
|
||||
import numpy as np
|
||||
import os
|
||||
import json
|
||||
|
||||
# 作为目标的转移概率矩阵
|
||||
P = np.array([
|
||||
[0.1, 0.3, 0.0, 0.6],
|
||||
[0.8, 0.0, 0.2, 0.0],
|
||||
[0.0, 0.9, 0.1, 0.0],
|
||||
[0.0, 0.3, 0.3, 0.4]
|
||||
])
|
||||
|
||||
# 采样
|
||||
def sample(n_samples, n_states, start_state):
|
||||
states = [i for i in range(n_states)]
|
||||
# 状态转移序列
|
||||
X = []
|
||||
# 开始采样
|
||||
X.append(start_state)
|
||||
current = start_state
|
||||
for i in range(n_samples):
|
||||
next = np.random.choice(states, p=P[current])
|
||||
X.append(next)
|
||||
current = next
|
||||
#endfor
|
||||
return X
|
||||
|
||||
def save_file(X, file_name):
|
||||
# 把0123变成ABCD
|
||||
Y = [chr(x+65) for x in X]
|
||||
#print(Y)
|
||||
# 保存Y到文件
|
||||
json_list = json.dumps(Y)
|
||||
file = open(file_name, "w")
|
||||
file.write(json_list)
|
||||
file.close()
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
# 采样数量
|
||||
n_samples = 10000
|
||||
# 状态空间
|
||||
n_states = 4
|
||||
# 起始状态(从0开始)
|
||||
start_state = 1
|
||||
X = sample(n_samples, n_states, start_state)
|
||||
#print(X)
|
||||
# 保存文件
|
||||
root = os.path.split(os.path.realpath(__file__))[0]
|
||||
file_name = os.path.join(root, "CarData.txt")
|
||||
save_file(X, file_name)
|
||||
@@ -0,0 +1,35 @@
|
||||
import numpy as np
|
||||
import json
|
||||
import os
|
||||
|
||||
# 计算从a_i转移到a_j的概率
|
||||
def calculate_P(X, a_i, a_j):
|
||||
n_i = 0
|
||||
n_j = 0
|
||||
for x in range(len(X)-1):
|
||||
if a_i == X[x]:
|
||||
n_i += 1
|
||||
if a_j == X[x+1]:
|
||||
n_j += 1
|
||||
print(n_i, n_j, n_j/n_i)
|
||||
|
||||
def open_file(file_name):
|
||||
file = open(file_name, "r")
|
||||
lines = file.read()
|
||||
file.close()
|
||||
data_list = json.loads(lines)
|
||||
return data_list
|
||||
|
||||
if __name__ == "__main__":
|
||||
# 状态空间
|
||||
n_states = 4
|
||||
# 读取文件
|
||||
root = os.path.split(os.path.realpath(__file__))[0]
|
||||
file_name = os.path.join(root, "CarData.txt")
|
||||
data_list = open_file(file_name)
|
||||
# 把 ABCD 变成 0123
|
||||
X = [ord(x)-65 for x in data_list]
|
||||
# 计算转移矩阵
|
||||
calculate_P(X, 0, 1) # 1代表B店,0代表A店
|
||||
|
||||
|
||||
@@ -0,0 +1,37 @@
|
||||
import numpy as np
|
||||
import json
|
||||
import os
|
||||
|
||||
def calculate_Matrix(n_states, X):
|
||||
P_counter = np.zeros((n_states, n_states))
|
||||
for i in range(len(X)-1):
|
||||
a_i = X[i]
|
||||
a_j = X[i+1]
|
||||
P_counter[a_i, a_j] += 1
|
||||
#endfor
|
||||
# 计算各列之和
|
||||
sum = np.sum(P_counter, axis=1, keepdims=True)
|
||||
print("各个状态出现的次数:\n",sum)
|
||||
P = P_counter / sum
|
||||
return P
|
||||
|
||||
def open_file(file_name):
|
||||
file = open(file_name, "r")
|
||||
lines = file.read()
|
||||
file.close()
|
||||
data_list = json.loads(lines)
|
||||
return data_list
|
||||
|
||||
if __name__ == "__main__":
|
||||
# 状态空间
|
||||
n_states = 4
|
||||
# 读取文件
|
||||
root = os.path.split(os.path.realpath(__file__))[0]
|
||||
file_name = os.path.join(root, "CarData.txt")
|
||||
data_list = open_file(file_name)
|
||||
# 把 ABCD 变成 0123
|
||||
X = [ord(x)-65 for x in data_list]
|
||||
# 计算转移矩阵
|
||||
P = calculate_Matrix(n_states, X)
|
||||
print("概率转移矩阵:")
|
||||
print(np.around(P, 1))
|
||||
@@ -0,0 +1,19 @@
|
||||
import numpy as np
|
||||
|
||||
P = np.array([
|
||||
[0.1, 0.3, 0.0, 0.6],
|
||||
[0.8, 0.0, 0.2, 0.0],
|
||||
[0.0, 0.9, 0.1, 0.0],
|
||||
[0.0, 0.3, 0.3, 0.4]
|
||||
])
|
||||
|
||||
def calculate_day(X, P, day):
|
||||
X_curr = X.copy()
|
||||
for i in range(day):
|
||||
print(str.format("day {0}: {1} ", i, X_curr))
|
||||
X_next = np.dot(X_curr, P)
|
||||
X_curr = X_next.copy()
|
||||
|
||||
if __name__=="__main__":
|
||||
X = np.array([0,1,0,0])
|
||||
calculate_day(X, P, 6)
|
||||
@@ -0,0 +1,23 @@
|
||||
import numpy as np
|
||||
|
||||
P = np.array([
|
||||
[0.1, 0.3, 0.0, 0.6],
|
||||
[0.8, 0.0, 0.2, 0.0],
|
||||
[0.0, 0.9, 0.1, 0.0],
|
||||
[0.0, 0.3, 0.3, 0.4]
|
||||
])
|
||||
|
||||
# 计算K步转移概率矩阵
|
||||
def K_step_matrix(P, K):
|
||||
Pk=P.copy()
|
||||
for i in range(K-1):
|
||||
Pk=np.dot(P,Pk)
|
||||
#print(Pk)
|
||||
return Pk
|
||||
|
||||
if __name__=="__main__":
|
||||
X = np.array([0,1,0,0])
|
||||
P5 = K_step_matrix(P, 5)
|
||||
print("5步转移矩阵:\n", P5)
|
||||
X5 = np.dot(X, P5)
|
||||
print("第 5 天的情况:", X5)
|
||||
@@ -0,0 +1,22 @@
|
||||
import numpy as np
|
||||
|
||||
P = np.array([
|
||||
[0.1, 0.3, 0.0, 0.6],
|
||||
[0.8, 0.0, 0.2, 0.0],
|
||||
[0.0, 0.9, 0.1, 0.0],
|
||||
[0.0, 0.3, 0.3, 0.4]
|
||||
])
|
||||
|
||||
def Check_Convergence(P):
|
||||
P_curr = P.copy()
|
||||
for i in range(100000):
|
||||
P_next=np.dot(P,P_curr)
|
||||
print("迭代次数 =",i+1)
|
||||
print(P_next)
|
||||
if np.allclose(P_curr, P_next):
|
||||
break
|
||||
P_curr = P_next
|
||||
return P_next
|
||||
|
||||
if __name__=="__main__":
|
||||
Pn = Check_Convergence(P)
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 43 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 79 KiB |
@@ -0,0 +1,439 @@
|
||||
|
||||
|
||||
## 租车问题
|
||||
|
||||
### 1 提出问题
|
||||
|
||||
在北京有一个叫做“四通八达”的租车连锁公司,在海淀区、朝阳区、通州区、丰台区开设了 4 个门店,编号分别为 A,B,C,D。客户可以从任意一个店租车,用完后可以归还到任意一个店。为了简化问题,我们假设客户都是早晨租车开走,晚上归还。
|
||||
|
||||
为了招徕顾客,公司特地买了一辆高档宝驴车,以普通价格出租。公司经理想知道如果有一天这辆车从B号店出租了,2 天后的早晨会最有可能在哪个店出现?5 天后又会如何?
|
||||
|
||||
公司统计了去年全年的租车、还车记录,抽样以后发现车辆在四个店的流动情况是这样的:
|
||||
|
||||
```
|
||||
["B", "A", "A", "D", "D", "C", "B", "A",......]
|
||||
```
|
||||
|
||||
数据序列的含义是,有一辆车:
|
||||
|
||||
- 第一天从B店被租走,晚上还到A店;
|
||||
- 第二天从A店被租走,晚上还到A店;
|
||||
- 第三天从A店被租走,晚上还到D店;
|
||||
- 第四天从D店被租走,晚上还到D店;
|
||||
- 第五天从D店被租走,晚上还到C店;
|
||||
- ......
|
||||
|
||||
序列中一共是10000条记录,可以认为是100辆车在100天内的借还地点数据的集合。这其中有什么规律呢?
|
||||
|
||||
其实这是一个标准的**转移概率**的问题。
|
||||
|
||||
### 2 转移概率
|
||||
|
||||
转移概率是马尔可夫链中的重要概念,若马氏链分为 m 个状态组成(在本例中 $m=4$),历史数据统计是由这 m 个状态所组成的序列。从任意一个状态出发,经过任意一次转移,必然出现状态 $1,2,……,m$ 中的一个,这种状态之间的转移称为转移概率。
|
||||
|
||||
在本问题中,有 4 个门店,所以一共就有 4 个状态。
|
||||
|
||||
如何根据历史数据计算出转移概率呢?其实就是条件概率:
|
||||
$$
|
||||
P(A|B)=\frac{P(A,B)}{P(B)} \tag{1}
|
||||
$$
|
||||
|
||||
具体到本问题中,公式如下:
|
||||
|
||||
$$
|
||||
P_{i,j}=\Pr\{X_{t+1}=a_j|X_t = a_i\}=\frac{\Pr\{X_{t+1}=a_j,X_t=a_i\}}{\Pr\{X_t=a_i\}}=\frac{a_j在a_i后出现的次数}{a_i出现的总次数} \tag{2}
|
||||
$$
|
||||
|
||||
公式 2 的含义是,在某一时刻 $t$,处于 $a_i$ 状态时,下一个时刻 $t+1$ 时,转移到 $a_j$ 状态的概率。
|
||||
|
||||
比较公式 2 和 1,形式上是完全相同的,只是里面的变量更具体化一些。
|
||||
|
||||
用上面的具体数据举例说明,如果我们想统计从 A 店租车,在 B 店还车的概率,应该这样做:
|
||||
|
||||
1. 设置两个计数器,一个 $nA = 0$,一个 $nB = 0$;
|
||||
1. 在上述序列数据中先找到所有的 A,出现一次就记录 $nA = nA+1$,假设一共出现 1000 次;
|
||||
2. 再找到紧接着 A 后面出现的 B,如果有,则计数器 $nB=nB+1$,假设一共 110 次;
|
||||
- 如果是 ['A','B','C'] 的顺序,则 $nA$ 计数,$nB$ 也计数。
|
||||
- 如果是 ['A','C','B'] 的顺序,中间有个'C',则 $nA$ 计数,$nB$ 不计数。
|
||||
3. 所有数据都统计完后计算:$P_{A,B} = \frac{nB}{nA}=\frac{110}{1000}=0.11$
|
||||
|
||||
```Python
|
||||
# 计算从a_i转移到a_j的概率
|
||||
def calculate_P(X, a_i, a_j):
|
||||
n_i = 0
|
||||
n_j = 0
|
||||
for x in range(len(X)-1):
|
||||
if a_i == X[x]:
|
||||
n_i += 1
|
||||
if a_j == X[x+1]:
|
||||
n_j += 1
|
||||
print(n_i, n_j, n_j/n_i)
|
||||
#end def
|
||||
calculate_P(X, 0, 1) # 0代表A店, 1代表B店
|
||||
```
|
||||
```
|
||||
[Out]
|
||||
2703 819 0.3029966703662597
|
||||
```
|
||||
上述代码统计出 A 店一共出现 2703 次,B 紧接在 A 后出现 819 次,所以 $P_{A,B}=\frac{819}{2703}\approx 0.3$。
|
||||
|
||||
同理可以统计出所有的转移概率,以 A 店为例:
|
||||
|
||||
- 从A号店租车后
|
||||
- 还到A号店的概率是0.1
|
||||
- 还到B号店的概率是0.3
|
||||
- 还到C号店的概率是0.0
|
||||
- 还到D号店的概率是0.6
|
||||
|
||||
当然,也统计了从 B,C,D 店租车还车的记录,绘制在图 1 中,避免赘述。
|
||||
|
||||
<img src="./img/Car1.png" width="500">
|
||||
|
||||
图 1
|
||||
|
||||
### 3 解决经理的问题
|
||||
|
||||
重复一下需求:公司经理想知道如果有一天(命名为第 0 天)这辆车从 B 号店出租了,2 天后的早晨会最有可能在哪个店出现?5 天后又会如何?
|
||||
|
||||
一般情况下,读者会根据图 1 从 B 店出发,根据概率顺藤摸瓜地计算出第 1 天的情况,再计算出第 2 天的情况。
|
||||
|
||||
- 第 1 天早晨出现在各门店的概率是:$[0.8,\ 0.0,\ 0.2,\ 0.0]$。
|
||||
- 第 2 天早晨出现在各门店的概率是......有点儿复杂,我们绘制出图 2 来帮助整理思路。
|
||||
|
||||
<img src="./img/Car2.png" width="600">
|
||||
|
||||
图 2
|
||||
|
||||
从图2一眼就可以看出来,第 2 天早晨该车出现在各门店的概率就是两个连续的概率之乘积,比如:
|
||||
- 第 1 天
|
||||
- 出现在A店(橙色)的概率是 0.8;
|
||||
- 出现在C店(红色)的概率是 0.2;
|
||||
- 但是不可能出现在B,D店;
|
||||
|
||||
- 第 2 天
|
||||
- 由于A店有0.3的概率还到B店,所以出现在B店的概率是$0.8 \times 0.3=0.24$;
|
||||
- 由于C店有0.9的概率还到B店,所以出现在B店的概率是$0.2 \times 0.9=0.18$;
|
||||
|
||||
所以,该车第3天早晨出现在B店的概率是 $0.24+0.18=0.42$。出现在其它店的数字也可以同理得到。
|
||||
|
||||
其它店的所有情况列在表 1 中,便于统计,数字的颜色和图2 是一一对应的,表示是哪个门店,方便读者对照理解。
|
||||
|
||||
表1
|
||||
|
||||
|从$\rightarrow$到|A店|B店|C店|D店|第1天|
|
||||
|:-:|-|-|-|-|-|
|
||||
|**A店**|$\color{orange}{0.8\times0.1=0.08}$|$\color{orange}{0.8\times0.3=0.24}$|$\color{orange}{0.8\times0.0=0.00}$|$\color{orange}{0.8\times0.6=0.48}$|$\color{orange}{0.8}$|
|
||||
|**B店**|$\color{green}{0.0\times0.8=0.0}$|$\color{green}{0.0\times0.0=0.0}$|$\color{green}{0.0\times0.2=0.0}$|$\color{green}{0.0\times0.0=0.0}$|$\color{green}{0.0}$|
|
||||
|**C店**|$\color{red}{0.2\times0.0=0.00}$|$\color{red}{0.2\times0.9=0.18}$|$\color{red}{0.2\times0.1=0.02}$|$\color{red}{0.2\times0.0=0.00}$|$\color{red}{0.2}$|
|
||||
|**D店**|$\color{blue}{0.0\times0.0=0.0}$|$\color{blue}{0.0\times0.3=0.0}$|$\color{blue}{0.0\times0.3=0.0}$|$\color{blue}{0.0\times0.4=0.0}$|$\color{blue}{0.0}$|
|
||||
|**第2天**|$\color{orange}{0.08}$|$\color{green}{0.24+0.18=0.42}$|$\color{red}{0.02}$|$\color{blue}{0.48}$|$1.0$|
|
||||
|
||||
数据解读:
|
||||
|
||||
- 表 1 除去表头,中间部分的 4x4 区域,和图2的数据是一致的。
|
||||
- 最后一列是中间4列的和,所以正好是第1天的出现概率。
|
||||
- 最后一行是中间4行的和,所以是第2天的出现概率。
|
||||
|
||||
OK! 2 天后的问题解决了,那么 5 天后呢?这么计算太麻烦了,我们引入转移概率矩阵的概念来帮助解决问题。
|
||||
|
||||
|
||||
### 4 转移概率矩阵
|
||||
|
||||
把图1变成表2,方便使用计算机来解决问题。
|
||||
|
||||
表2
|
||||
|
||||
|P: 从$\rightarrow$到|A|B|C|D||输出总和|
|
||||
|:-:|-|-|-|-|-|-|
|
||||
|**A**|0.1|0.3|0.0|0.6||1.0|
|
||||
|**B**|0.8|0.0|0.2|0.0||1.0|
|
||||
|**C**|0.0|0.9|0.1|0.0||1.0|
|
||||
|**D**|0.0|0.3|0.3|0.4||1.0|
|
||||
||||||||
|
||||
|输入总和|0.9|1.5|0.6|1.0|4.0|
|
||||
|
||||
- 是一个 $n \times n$ 的方阵(中间的 4x4 区域)。
|
||||
- 矩阵各元素都是非负的,即:$0 \le P_{i,j} \le 1$。
|
||||
- 各行元素(输出概率)之和为 1,即:$\sum^n_{j=1} P_{i,j}=1$。
|
||||
- 对于输入概率总和没有限制,但是最终的总和(右下角)在行列上肯定应该一致,都是4.0(因为有个 4 个门店),即:$\sum_{i=1}^n \sum_{j=1}^n P_{i,j} = n$。
|
||||
|
||||
如果从原始数据计算的话,可以用更简单的代码一次搞定所有门店的数据,不需要像前面的代码那样一个个地计算:
|
||||
|
||||
```Python
|
||||
def calculate_Matrix(n_states, X):
|
||||
P_counter = np.zeros((n_states, n_states))
|
||||
for i in range(len(X)-1):
|
||||
a_i = X[i]
|
||||
a_j = X[i+1]
|
||||
P_counter[a_i, a_j] += 1
|
||||
#endfor
|
||||
# 计算各列之和
|
||||
sum = np.sum(P_counter, axis=1, keepdims=True)
|
||||
print("各个状态出现的次数:\n",sum)
|
||||
P = P_counter / sum
|
||||
return P
|
||||
#enddef
|
||||
|
||||
# 把原始数据中的字符 ABCD 变成 0123
|
||||
n_states= 4
|
||||
X = [ord(x)-65 for x in data_list]
|
||||
# 计算转移矩阵
|
||||
P = calculate_Matrix(n_states, X)
|
||||
print("转移概率矩阵:")
|
||||
print(np.around(P, 1)) # 保留1位小数
|
||||
```
|
||||
```
|
||||
[Out]
|
||||
各个状态出现的次数:
|
||||
[[2703.]
|
||||
[3040.]
|
||||
[1568.]
|
||||
[2689.]]
|
||||
转移概率矩阵:
|
||||
[[0.1 0.3 0. 0.6]
|
||||
[0.8 0. 0.2 0. ]
|
||||
[0. 0.9 0.1 0. ]
|
||||
[0. 0.3 0.3 0.4]]
|
||||
```
|
||||
|
||||
|
||||
|
||||
用矩阵来表示
|
||||
|
||||
$$
|
||||
P =
|
||||
\begin{pmatrix}
|
||||
P_{11} & P_{12} & P_{13} & P_{14}
|
||||
\\
|
||||
P_{21} & P_{22} & P_{23} & P_{24}
|
||||
\\
|
||||
P_{31} & P_{32} & P_{33} & P_{34}
|
||||
\\
|
||||
P_{41} & P_{42} & P_{43} & P_{44}
|
||||
\end{pmatrix}=
|
||||
\begin{pmatrix}
|
||||
0.1 & 0.3 & 0.0 & 0.6
|
||||
\\
|
||||
0.8 & 0.0 & 0.2 & 0.0
|
||||
\\
|
||||
0.0 & 0.9 & 0.1 & 0.0
|
||||
\\
|
||||
0.0 & 0.3 & 0.3 & 0.4
|
||||
\end{pmatrix}
|
||||
\tag{3}
|
||||
$$
|
||||
|
||||
而对于一个通用的问题,矩阵形式是:
|
||||
|
||||
$$
|
||||
P =
|
||||
\begin{pmatrix}
|
||||
P_{11} & \cdots & P_{1n}
|
||||
\\
|
||||
\vdots & \ddots & \vdots
|
||||
\\
|
||||
P_{n1} & \cdots & P_{nn}
|
||||
\end{pmatrix}
|
||||
\tag{4}
|
||||
$$
|
||||
|
||||
有了矩阵,我们可以方便地使用矩阵运算来代替复杂的循环逻辑
|
||||
|
||||
首先定义第 0 天的初始向量为:
|
||||
$$
|
||||
X_0=(0,1,0,0) \tag{5}
|
||||
$$
|
||||
|
||||
表示豪华宝驴车目前肯定在 B 门店。
|
||||
|
||||
那么第 1 天在哪个门店的问题可以转换为:
|
||||
|
||||
$$
|
||||
X_1 = X_0 P=(0, 1 , 0 , 0)
|
||||
\begin{pmatrix}
|
||||
0.1 & 0.3 & 0.0 & 0.6
|
||||
\\
|
||||
0.8 & 0.0 & 0.2 & 0.0
|
||||
\\
|
||||
0.0 & 0.9 & 0.1 & 0.0
|
||||
\\
|
||||
0.0 & 0.3 & 0.3 & 0.4
|
||||
\end{pmatrix}=
|
||||
(0.8,\ 0.0,\ 0.2,\ 0.0)
|
||||
\tag{6}
|
||||
$$
|
||||
|
||||
那么第2天在哪个门店的问题可以递归计算出:
|
||||
|
||||
$$
|
||||
X_2 = X_1 P=(0.8,\ 0.0,\ 0.2,\ 0.0)
|
||||
\begin{pmatrix}
|
||||
0.1 & 0.3 & 0.0 & 0.6
|
||||
\\
|
||||
0.8 & 0.0 & 0.2 & 0.0
|
||||
\\
|
||||
0.0 & 0.9 & 0.1 & 0.0
|
||||
\\
|
||||
0.0 & 0.3 & 0.3 & 0.4
|
||||
\end{pmatrix}=
|
||||
(0.08, \ 0.42, \ 0.02, \ 0.48) \tag{7}
|
||||
$$
|
||||
|
||||
式7的结果和表1完全一致。
|
||||
|
||||
推广到一般情况,计算第 $n+1$ 天的概率时,需要使用第 $n$ 天的结果
|
||||
|
||||
$$
|
||||
X_{n+1}=X_n P \tag{8}
|
||||
$$
|
||||
|
||||
|
||||
我们可以很方便地写出代码
|
||||
|
||||
```Python
|
||||
P = np.array([
|
||||
[0.1, 0.3, 0.0, 0.6],
|
||||
[0.8, 0.0, 0.2, 0.0],
|
||||
[0.0, 0.9, 0.1, 0.0],
|
||||
[0.0, 0.3, 0.3, 0.4]
|
||||
])
|
||||
|
||||
def calculate_day(P, day):
|
||||
X = np.array([0,1,0,0])
|
||||
X_curr = X.copy()
|
||||
for i in range(day):
|
||||
print(str.format("day {0}: {1} ", i+1, X_curr))
|
||||
X_next = np.dot(X_curr, P)
|
||||
X_curr = X_next
|
||||
#enddef
|
||||
|
||||
calculate_day(P, 5) # 5表示计算5天的数据
|
||||
```
|
||||
```
|
||||
[Out]
|
||||
day 0: [0 1 0 0]
|
||||
day 1: [0.8 0. 0.2 0. ]
|
||||
day 2: [0.08 0.42 0.02 0.48]
|
||||
day 3: [0.344 0.186 0.23 0.24 ]
|
||||
day 4: [0.1832 0.3822 0.1322 0.3024]
|
||||
day 5: [0.32408 0.26466 0.18038 0.23088]
|
||||
```
|
||||
|
||||
OK!经理的问题解决了,在第 5 天的时候,该豪华宝驴车在 4 个门店出现的概率依次是:
|
||||
|
||||
$[0.32408,\ 0.26466,\ 0.18038,\ 0.23088]$
|
||||
|
||||
### 5 K-步转移概率矩阵
|
||||
|
||||
从上面的例子我们可以看到,只有一步的转移概率通常是不能满足实际需要的,通常我们需要知道 K 步后的情况如何。
|
||||
|
||||
比如,5步后的概率应该是:
|
||||
|
||||
$$
|
||||
X_5 = X_4 P=(X_3P)P=((X_2P)P)P=(((X_1P)P)P)P=((((X_0P)P)P)P)P=X_0 P^5 \tag{9}
|
||||
$$
|
||||
|
||||
式 4 可以定义为一步转移概率矩阵,那么从式9的推导可以得到式 10,定义为 K 步转移概率:
|
||||
|
||||
$$
|
||||
P^K =
|
||||
\begin{pmatrix}
|
||||
P_{11} & \cdots & P_{1n}
|
||||
\\
|
||||
\vdots & \ddots & \vdots
|
||||
\\
|
||||
P_{n1} & \cdots & P_{nn}
|
||||
\end{pmatrix}^K
|
||||
\tag{10}
|
||||
$$
|
||||
|
||||
代码如下
|
||||
|
||||
```Python
|
||||
# 计算K步转移概率矩阵
|
||||
def K_step_matrix(P, K):
|
||||
Pk=P.copy()
|
||||
for i in range(K-1):
|
||||
Pk=np.dot(Pk,P)
|
||||
return Pk
|
||||
```
|
||||
读者可能会注意到循环次数是 K-1,而不是 K。假设我们计算 K=2 步转移矩阵,那么一步矩阵P和自己做一次矩阵相乘就可以了,而不是 2 次。所以 K 步矩阵做 K-1 次矩阵相乘。另外一个细节是 np.dot(P,Pk) 和 np.dot(Pk,P) 都会得到相同的结果。
|
||||
|
||||
得到 K 步转移矩阵后,可以用如下代码简单地直接计算第 5 天的情况:
|
||||
|
||||
```Python
|
||||
if __name__=="__main__":
|
||||
X = np.array([0,1,0,0])
|
||||
P5 = K_step_matrix(P, 5)
|
||||
print(P5)
|
||||
X5 = np.dot(X, P5)
|
||||
print(X5)
|
||||
```
|
||||
```
|
||||
[Out]
|
||||
5步转移矩阵:
|
||||
[[0.25585 0.31989 0.14028 0.28398]
|
||||
[0.32408 0.26466 0.18038 0.23088]
|
||||
[0.19728 0.36459 0.14005 0.29808]
|
||||
[0.26448 0.29469 0.15843 0.2824 ]]
|
||||
第 5 天的情况: [0.32408 0.26466 0.18038 0.23088]
|
||||
```
|
||||
|
||||
与上面的迭代方法计算结果一致。
|
||||
|
||||
|
||||
|
||||
### 6 更多步的情况
|
||||
|
||||
从 K 步转移可以向更远的地方思考,如果 K 为无穷大时,这个转移概率矩阵是什么情况呢?
|
||||
|
||||
我们用代码做一个试验
|
||||
|
||||
```Python
|
||||
def N_step_matrix(P):
|
||||
P_curr = P.copy()
|
||||
for i in range(100000):
|
||||
P_next=np.dot(P,P_curr)
|
||||
print("迭代次数=",i+1)
|
||||
print(P_next)
|
||||
if np.allclose(P_curr, P_next):
|
||||
break
|
||||
P_curr = P_next
|
||||
return P_next
|
||||
#enddef
|
||||
N_step_matrix(P)
|
||||
```
|
||||
注意,在代码中我们使用 np.allclose(P_curr, P_next) 来判断上一次迭代的结果和本次迭代的结果的差值,如果小于1e-6则认为已经收敛,停止循环。运行过程如下所示:
|
||||
```
|
||||
[Out]
|
||||
迭代次数 = 1
|
||||
[[0.25 0.21 0.24 0.3 ]
|
||||
[0.08 0.42 0.02 0.48]
|
||||
[0.72 0.09 0.19 0. ]
|
||||
[0.24 0.39 0.21 0.16]]
|
||||
迭代次数 = 2
|
||||
[[0.193 0.381 0.156 0.27 ]
|
||||
[0.344 0.186 0.23 0.24 ]
|
||||
[0.144 0.387 0.037 0.432]
|
||||
[0.336 0.309 0.147 0.208]]
|
||||
|
||||
......
|
||||
|
||||
迭代次数 = 27
|
||||
[[0.26966338 0.30337037 0.1573036 0.26966265]
|
||||
[0.26966195 0.30337167 0.15730289 0.26966349]
|
||||
[0.26966412 0.3033697 0.15730396 0.26966222]
|
||||
[0.26966285 0.30337085 0.15730334 0.26966296]]
|
||||
迭代次数 = 28
|
||||
[[0.26966264 0.30337105 0.15730323 0.26966309]
|
||||
[0.26966353 0.30337024 0.15730367 0.26966257]
|
||||
[0.26966217 0.30337147 0.157303 0.26966336]
|
||||
[0.26966296 0.30337075 0.15730339 0.2696629 ]]
|
||||
```
|
||||
|
||||
我们惊奇地发现,迭代了28次后就已经收敛到很小的误差了。可以认为在经过多次的 “租、还、租、还” 循环后,某辆车在四个门店的出现概率是:$[0.27,\ 0.30,\ 0.16,\ 0.27]$,精确到两位小数。
|
||||
|
||||
从实际情况来看:
|
||||
- 可能是因为通州区太远了,所以顾客都偏好与就近还车。
|
||||
- 海淀区和朝阳区的停车场要大一些,以便可以停更多的车。
|
||||
- 丰台区无所谓,因为属于地广人稀的地段,去那里还车的大概都是去机场顾客。
|
||||
Reference in New Issue
Block a user