练习1
某地区有甲、乙、丙三家食品厂生产同一种食品，有一千个用户（或购货点），假定在研究期间无新用户加入也无老用户退出，只有用户的转移，已知 2006 年 5 月份有 500 户是甲厂的顾客；400 户是乙厂的顾客；100 户是丙厂的顾客。6 月份，甲厂有 400 户原来的顾客，上月的顾客有 50 户转乙厂，50 户转丙厂；乙厂有 300 户原来的顾客，上月的顾客有 20 户转甲厂，80 户转丙厂；丙厂有 80 户原来的顾客，上月的顾客有 10 户转甲厂，10 户转乙厂。计算其状态转移概率与平稳分布。
由题意得 6 月份顾客转移表，如下：
| 厂名 | 甲 | 乙 | 丙 | 合计 |
|---|---|---|---|---|
| 甲 | 400 | 50 | 50 | 500 |
| 乙 | 20 | 300 | 80 | 400 |
| 丙 | 10 | 10 | 80 | 100 |
| 合计 | 430 | 360 | 210 | 1000 |
import numpy as np
# 原始数据矩阵
data = np.array([
[400, 50, 50],
[20, 300, 80],
[10, 10, 80]
])
# 每行的合计数
row_totals = np.array([500, 400, 100]).reshape(-1, 1)
# 将数据矩阵每行除以合计数，归一化处理为转移矩阵
P = data / row_totals
print("转移矩阵 P:")
print(np.round(P, 2))
# 计算平稳分布
eigvals, eigvecs = np.linalg.eig(P.T)
eigvec1 = eigvecs[:, np.isclose(eigvals, 1)]
# 将特征向量归一化以使其总和为1
stationary = eigvec1 / eigvec1.sum()
stationary = stationary.real.flatten()
stationary = np.round(stationary, 2)
print("\n平稳分布:")
print(stationary)
#[0.29 0.29 0.43]
练习2
计算机运行状态的预测。设一随机系统状态空间, 总共有4种状态, 记录观测系统所处状态如下:
\[\begin{array}{llllllllll} 4 & 3 & 2 & 1 & 4 & 3 & 1 & 1 & 2 & 3 \\ 2 & 1 & 2 & 3 & 4 & 4 & 3 & 3 & 1 & 1 \\ 1 & 3 & 3 & 2 & 1 & 2 & 2 & 2 & 4 & 4 \\ 2 & 3 & 2 & 3 & 1 & 1 & 2 & 4 & 3 & 1 \end{array} \]若该系统可用马尔科夫模型描述, 且当前系统处于状态4 , 则下一步最可能是状态为何？
求解思路：
- 求解一步转移概率: 从状态 \(i\) 到状态 \(j\) 的概率的估计值为
例如, 从状态1可以转移到状态 \(1 、 2 、 3 、 4\), 那么状态1转移到其他状态的总次数, 就是观测到的数据中, 出现 “ 11 ”" 12 ”"13" "14" 的总次数, 也就是公式中的分母。
从状态 \(i\) 到状态 \(j\) 的次数统计（不需要手算, 代码部分会讲）
| 1 | 2 | 3 | 4 | $$\sum_{j=1}^4 n_{ij}$$ | |
|---|---|---|---|---|---|
| 1 | 4 | 4 | 1 | 1 | 10 |
| 2 | 3 | 2 | 4 | 2 | 11 |
| 3 | 4 | 4 | 2 | 1 | 11 |
| 4 | 0 | 1 | 4 | 2 | 7 |
- 求出一步转移矩阵:
其中第 \(i\) 行第 \(j\) 列的数, 代表着从状态 \(i\) 一步转移到状态 \(j\) 的概率的估计值。题目说当前状态为 4 , 则下一步的状态最有可能是3。
- 则第 \(\mathrm{n}\) 步的状态的概率分布为:
初始概率分布行向量（即初始时每种状态出现的概率）记做 \(P^{(0)}\), 一步转移概率矩阵记做 \(P\)。
注意是把矩阵\(P\)乘 \(n\) 次。
系统初始化后运行到第2步的概率分布: \(P^{(2)}=P^{(0)} P^2=[0.291,0.289,0.271,0.149]\)
系统初始化后运行到第5步的概率分布: \(P^{(5)}=P^{(0)} P^5=[0.294,0.289,0.268,0.149]\)
- 代码求解
求一步转移矩阵:一般需要根据已知数据，求出一步转移矩阵的估计值；
需要遍历已知数据，统计每一组相邻两个状态的数量（找出有多少个 "1 1"，“12",..... )；
如果题目所定的系统状态是文字（例如张三的四种状态），那么在写代码时，可以自行定义不同状态对应的数字，给出一步转移概率和转移矩阵:
import numpy as np
# 原始数据矩阵（去掉最后一列）
data = np.array([
[4, 4, 1, 1],
[3, 2, 4, 2],
[4, 4, 2, 1],
[0, 1, 4, 2]
])
# 计算每行的总和
sums = data.sum(axis=1)
# 计算状态转移矩阵 P
P = data / sums[:, None]
print("状态转移矩阵P:")
print(np.round(P, 2))
# 计算 P 的 2 次方
P2 = np.linalg.matrix_power(P, 2)
print("\nP 的 2 次方:")
print(np.round(P2, 2))
# 计算 P 的 5 次方
P5 = np.linalg.matrix_power(P, 5)
print("\nP 的 5 次方:")
print(np.round(P5, 2))
# 计算平稳分布
eigvals, eigvecs = np.linalg.eig(P.T)
eigvec1 = eigvecs[:, np.isclose(eigvals, 1)]
# 将特征向量归一化以使其总和为1
stationary = eigvec1 / eigvec1.sum()
stationary = stationary.real.flatten()
stationary = np.round(stationary, 2)
print("\n平稳分布:")
print(stationary)
P 的 5 次方:
[[0.29 0.29 0.27 0.15]
[0.29 0.29 0.27 0.15]
[0.29 0.29 0.27 0.15]
[0.29 0.29 0.27 0.15]]
平稳分布:
[0.29 0.29 0.27 0.15]
练习3
考虑中国移动，中国联通和中国电信三家运营商，他们的用户可以互相携号转网，已经当前3家运营商的市场份额，而且也能测试出用户转网的可能性，那么将来3家运营商的市场份额情况如何，即利用当前已知的两项数据，分别是当前的市场份额、用户接下来使用运营商的可能性（即转移概率矩阵），则可预测将来3家运营商的市场份额情况。当前3家运营商分别的市场占比为0.55,0.25和0.2，但是有携号转网政策后，用户可自由的切换运营商，当前从调查数据中可得到，移动用户预期下年还会使用移动的比例是80%，使用电信的比例是15%，联通的比例是5%。电信用户接下来会使用移动的比例是20%，继续使用电信的比例是78%，使用联通的比例是2%。联通用户接下来使用移动的比例是5%，电信为20%，继续使用联通的比例是75%。现希望预期后续比如10年后，3家运营商的市场份额情况如何。
import numpy as np
# 初始市场份额
initial_distribution = np.array([0.55, 0.25, 0.2])
# 状态转移矩阵
P = np.array([
[0.80, 0.15, 0.05],
[0.20, 0.78, 0.02],
[0.05, 0.20, 0.75]
])
# 计算10年后的市场份额分布
distribution_10_years = initial_distribution.dot(np.linalg.matrix_power(P, 10))
# 状态名
states = ["移动", "电信", "联通"]
print("10年后的市场份额分布:")
for state, share in zip(states, np.round(distribution_10_years, 2)):
print(f"{state}: {share}")
# 计算平稳分布
eigvals, eigvecs = np.linalg.eig(P.T)
eigvec1 = eigvecs[:, np.isclose(eigvals, 1)]
# 将特征向量归一化以使其总和为1
stationary = eigvec1 / eigvec1.sum()
stationary = stationary.real.flatten()
stationary = np.round(stationary, 2)
print("\n平稳分布:")
for state, share in zip(states, stationary):
print(f"{state}: {share}")
10年后的市场份额分布:
移动: 0.45
电信: 0.42
联通: 0.13
平稳分布:
移动: 0.45
电信: 0.42
联通: 0.12