这是Jupyter Notebook转换成的markdown。
部分内容参考了这里的答案。
import numpy as np
from sympy import Matrix
第一题
假设传感器并不会在使用过程中损坏，只可能一开始就是坏的或者好的。令随机变量\(X\)，其值为1表示传感器是好的，为0表示传感器是坏的
\[\begin{aligned} p(X=0\mid z_{1:n} < 1)&=\frac{p(z_{1:n} < 1\mid X=0)p(X=0)}{p(z_{1:n} < 1)}\\ &=\frac{1\times 0.01}{p(z_{1:n} < 1\mid X=0)p(X=0) + p(z_{1:n} < 1\mid X=1)p(X=1)}\\ &=\frac{0.01}{0.01 + \frac{1}{3^N} \times 0.99} \end{aligned} \]第二题
第1问
查查表很容易得到：
\(0.2\times0.4\times0.2=0.016\)
第2、3问
要随机生成天气，只要取出每天概率分布中概率最大的天气就好，现在主要求一求随着天数增加，分布会不会稳定：
tran_mat = np.array(
[[0.8, 0.4, 0.2],
[0.2, 0.4, 0.6],
[0, 0.2, 0.2]]
)
weathers = np.diag([1, 1, 1], 0)
wea_str = ["sunny", "cloudy", "rainy"]
for idx in range(3):
print("Let 1st day's weather to be ", wea_str[idx], "!", sep='')
print("day", "sunny", "cloudy", "rainy", "predict", sep="\t")
weather = np.reshape(weathers[idx], (3, 1))
for day in range(12): # 输出10天的结果
weather = np.matmul(tran_mat, weather)
if (day >= 8):
print(day+2, end='\t')
for w in weather:
print("%.3f" % w, end='\t')
print(wea_str[np.argmax(weather)])
print("================")
Let 1st day's weather to be sunny!
day sunny cloudy rainy predict
10 0.643 0.285 0.071 sunny
11 0.643 0.286 0.071 sunny
12 0.643 0.286 0.071 sunny
13 0.643 0.286 0.071 sunny
================
Let 1st day's weather to be cloudy!
day sunny cloudy rainy predict
10 0.642 0.286 0.072 sunny
11 0.642 0.286 0.072 sunny
12 0.643 0.286 0.072 sunny
13 0.643 0.286 0.071 sunny
================
Let 1st day's weather to be rainy!
day sunny cloudy rainy predict
10 0.642 0.286 0.072 sunny
11 0.642 0.286 0.072 sunny
12 0.643 0.286 0.072 sunny
13 0.643 0.286 0.071 sunny
================
第4问
这是一个马尔科夫过程（我并不是不是很了解这个问题），要求闭式解主要是要求状态转移矩阵的幂，需要用到线性代数中关于相似对角化的知识，此处就直接用Python的包求解啦：
from sympy import Matrix
# 还是原来的状态转移矩阵，只不过换一个数据类型
sympy_tran = Matrix(3, 3,
[0.8, 0.4, 0.2,
0.2, 0.4, 0.6,
0, 0.2, 0.2])
# 求相似对角化矩阵和
(diagonalize_mat, diag_mat) = sympy_tran.diagonalize()
diag_mat
\(\displaystyle \left[\begin{matrix}1.0 & 0 & 0\\0 & 0.482842712474619 & 0\\0 & 0 & -0.082842712474619\end{matrix}\right]\)
假设第\(t\)天的状态为\(A_t\)，状态转移矩阵\(M\)可以相似对角化得到\(D=P^{-1}MP\)。观察上边的计算结果可得：
\[\lim_{N\to\infty}D^N= \begin{bmatrix} 1&0&0\\0&0&0\\0&0&0 \end{bmatrix} \]sta_tran = diagonalize_mat * Matrix(3, 3, [1, 0, 0, 0, 0, 0, 0, 0, 0])\
* diagonalize_mat.inv()
sta_tran
\(\displaystyle \left[\begin{matrix}0.642857142857143 & 0.642857142857143 & 0.642857142857143\\0.285714285714286 & 0.285714285714286 & 0.285714285714286\\0.0714285714285714 & 0.0714285714285714 & 0.0714285714285714\end{matrix}\right]\)
故有：
\[\begin{aligned} A_t &= (PDP^{-1})^NA_0\\ &=PD^NP^{-1}A_0\\ &=\begin{bmatrix} 0.643&0.643&0.643\\0.286&0.286&0.286\\0.071&0.071&0.071& \end{bmatrix}A_0 \end{aligned} \]又因为A_0各元素和为1，所以无论\(A_0\)如何取值，最后的稳态分布都是\(\begin{bmatrix}0.643\\0.286\\0.071\end{bmatrix}\)
第5问
sta_distrib = np.array([sta_tran[0, 0], sta_tran[1, 0], sta_tran[2, 0]], dtype=float)
entropy = -np.sum(np.log2(sta_distrib) * sta_distrib)
entropy
1.1981174211304033
故熵是1.198。
第6问
设某一天的天气是A，其后一天的天气是B，则已知\(p(B\mid A)\)，求\(p(A\mid B)\)。
\[\begin{aligned} p(A\mid B) &= \frac{p(B\mid A) p(A)}{p(B)}\\ \end{aligned} \]其中\(p(A)\)、\(p(B)\)即是稳态分布，有
def reverse_transition(a2b, p_a, p_b):
'''求状态转移矩阵的反过来的转移矩阵'''
b2a = np.zeros((3, 3))
for a in range(3):
for b in range(3):
b2a[a, b] = a2b[b, a] * p_a[a] / p_b[b]
return b2a
rev_tran_mat = reverse_transition(tran_mat, sta_distrib, sta_distrib)
print("to->ye", "sunny", "cloudy", "rainy", sep="\t")
for idx in range(3):
print(wea_str[idx]+"\t", "%.4f\t%.4f\t%.4f"
%(rev_tran_mat.T[idx, 0],
rev_tran_mat.T[idx, 1],
rev_tran_mat.T[idx, 2]))
to->ye sunny cloudy rainy
sunny 0.8000 0.1778 0.0222
cloudy 0.4500 0.4000 0.1500
rainy 0.0000 0.8000 0.2000
第6问
马尔可夫过程要求对当下的状态能够充分估计未来的状态，所以在只知道当前天气的条件下添加天气变量将使得以上推导不再有效，解决方法是在当前状态中再加上天气。
第三题
注意这题是建立在第2题的基础上的
记号约定：
- 晴天、多云和雨天分别记为1~3
- \(X_n=s\)（或\(X_n=c\)、\(X_n=r\)）简写为\(x_n\)
已知今天天气到明天天气的状态转移矩阵\(P_1\)和天气到传感器测量结果的转移矩阵\(P_2\)：
\[P_1= \begin{bmatrix} 0.8&0.4&0.2\\0.2&0.4&0.6\\0&0.2&0.2 \end{bmatrix},\ P_2= \begin{bmatrix} 0.6&0.3&0\\0.4&0.7&0\\0&0&1 \end{bmatrix} \]第1问
等价于在\(n=5,X_1=1,z_{2:n}=\{2,2,3,1\},i=1\)时求下式的值：
\[p(X_n=i\mid z_{2:n},x_1),\quad n\ge2 \]将该式化简可得
\[\begin{aligned} p(X_n=i\mid z_{2:n},x_1) &=\eta p(z_n\mid z_{2:n-1}, X_n=i, x_1)p(X_n=i\mid z_{2:n-1},x_1)\\ &=\eta p(z_n\mid X_n=i)\sum\limits_{j=1}^3p(X_n=i\mid z_{2:n-1}, X_{n-1}=j, x_1)p(X_{n-1}=j\mid z_{2:n-1}, x_1)\\ &=\eta p(z_n\mid X_n=i)\sum\limits_{j=1}^3p(X_n=i\mid X_{n-1}=j)p(X_{n-1}=j\mid z_{2:n-1}, x_1)\\ &=\eta p(z_n\mid X_n=i)\sum\limits_{j=1}^3p(X_n=i\mid X_{n-1}=j)p(X_{n-1}=j\mid z_{2:n-1}, x_1) \end{aligned} \]其中\(p(X_{n-1}=j\mid z_{2:n-1}, x_1)\)与\(p(X_n=i\mid z_{2:n},x_1)\)结构一致，可见可以按\(n=2\sim5\)的顺序求的该题的答案。
首先求迭代的初始条件
SUNNY = 0
CLOUDY = 1
RAINY = 2
# 当天天气到传感器的转移矩阵
wea2sen = np.array(
[[0.6, 0.3, 0],
[0.4, 0.7, 0],
[0, 0, 1]]
)
# 当天天气到后一天的天气的转移矩阵
wea2tom = np.array(
[[0.8, 0.4, 0.2],
[0.2, 0.4, 0.6],
[0, 0.2, 0.2]]
)
# 传感器天气序列，从第2天开始
sensor_data = [CLOUDY, CLOUDY, RAINY, SUNNY]
# 初始天气，第1天的真实天气
wea_init = SUNNY
# 本题条件下第2、3、4、5天的真实天气的概率分布
p_sync = np.zeros((3, 4))
# 第二天的真实天气概率分布
for wea in range(3):
p_sync[wea, 0] = wea2sen[sensor_data[0], wea] * wea2tom[wea, wea_init]
# 归一化，从而避免求p(z_n|z_{n-1})
p_sync[:, 0] /= np.sum(p_sync[:, 0])
p_sync
array([[0.69565217, 0. , 0. , 0. ],
[0.30434783, 0. , 0. , 0. ],
[0. , 0. , 0. , 0. ]])
最后迭代即可得到第五天的概率分布，满足的公式如下（上边推的此处再写一遍）：
\[p(X_n=i\mid z_{2:n},x_1)=\eta p(z_n\mid X_n=i)\sum\limits_{j=1}^3p(X_n=i\mid X_{n-1}=j)p(X_{n-1}=j\mid z_{2:n-1}, x_1) \]for day in range(3, 6):
for wea in range(3):
for sum_wea in range(3):
p_sync[wea, day-2] += wea2tom[wea, sum_wea] * p_sync[sum_wea, day-3]
p_sync[wea, day-2] *= wea2sen[sensor_data[day-2], wea]
p_sync[:, day-2] /= np.sum(p_sync[:, day-2])
p_sync
array([[0.69565217, 0.59770115, 0. , 0.4 ],
[0.30434783, 0.40229885, 0. , 0.6 ],
[0. , 0. , 1. , 0. ]])
第2问
方式1
此时每一天只能知道截止到当天的传感器信息，这和第1问求的结果形式一样，直接套用一遍：
# 传感器天气序列，从第2天开始
sensor_data_2 = [SUNNY, SUNNY, RAINY]
# 初始天气，第1天的真实天气
wea_init_2 = SUNNY
# 本题条件下第2、3、4天的真实天气的概率分布
p_sync_2 = np.zeros((3, 3))
# 第二天的真实天气概率分布
for wea in range(3):
p_sync_2[wea, 0] = wea2sen[sensor_data[0], wea] * wea2tom[wea, wea_init_2]
# 归一化，从而避免求p(z_n|z_{n-1})
p_sync_2[:, 0] /= np.sum(p_sync_2[:, 0])
# 迭代求各天的分布
for day in range(3, 5):
for wea in range(3):
for sum_wea in range(3):
p_sync_2[wea, day-2] += wea2tom[wea, sum_wea] * p_sync_2[sum_wea, day-3]
p_sync_2[wea, day-2] *= wea2sen[sensor_data[day-2], wea]
p_sync_2[:, day-2] /= np.sum(p_sync_2[:, day-2])
p_sync_2
array([[0.69565217, 0.59770115, 0. ],
[0.30434783, 0.40229885, 0. ],
[0. , 0. , 1. ]])
方式2
首先还是弄清楚要求什么，并将其化简：
\[\begin{aligned} p(X_k=i\mid z_{2:n},x_1) &=\sum\limits_{j=1}^3p(X_k=i\mid z_{2:n}, X_{k-1}=j,x_1)p(X_{k-1}=j\mid z_{2:n}, x_1)\\ &=\sum\limits_{j=1}^3p(X_k=i\mid z_{k:n}, X_{k-1}=j)p(X_{k-1}=j\mid z_{2:n}, x_1) \end{aligned} \]可以发现\(p(X_{k-1}=j\mid z_{2:n}, x_1)\)与待求的目标形式是相同的，所以思路遇上题一样，求出初始条件然后迭代。而初始条件$$p(X_2=i\mid z_{2:n},x_1)$$又和上式中的另一项\(p(X_k=i\mid z_{k:n}, X_{k-1}=j)\)是同构的，因此关键在于求得\(p(X_k=i\mid z_{k:n}, X_{k-1}=j)\)。将其化简：
\[\begin{aligned} p(X_k=i\mid z_{k:n}, X_{k-1}=j) &=\eta p(z_{k:n}\mid X_k=i, X_{k-1}=j)p(X_k=i\mid X_{k-1}=j)\\ &=\eta p(z_{k:n}\mid X_k=i)p(X_k=i\mid X_{k-1}=j) \end{aligned} \]故首先要求\(p(z_{k:n}\mid X_k=i)\)。
def get_sen_seq_prob(sen_seq, init_wea):
'''已知第一天的天气，求一组从第二天开始的天气序列发生的概率'''
prob = 1
wea = np.reshape(np.zeros(3), (3, 1))
wea[init_wea] = 1
for day in range(0, len(sen_seq)):
sen = np.matmul(wea2sen, wea)
prob *= sen[sen_seq[day], 0]
wea = np.matmul(wea2tom, wea)
return prob
p_sync_3 = np.zeros((3, 3))
# 得到第2天的分布
for wea in range(3):
p_sync_3[wea, 0] = get_sen_seq_prob(sensor_data_2, wea) * wea2tom[wea, wea_init]
p_sync_3[:, 0] /= np.sum(p_sync_3[:, 0])
# 得到随后几天的分布
for day in range(3, 5):
# 求$p(X_k=i\mid z_{k:n}, X_{k-1}=j)$
p_temp = np.zeros((3, 3))
for wea_sum in range(3):
for wea in range(3):
p_temp[wea, wea_sum] = get_sen_seq_prob(sensor_data_2[day-2:], wea) \
* wea2tom[wea, wea_sum]
# 有的序列无论初始天气如何，都不可能达成，这时不能进行归一化
if (p_temp[:, wea_sum].any() > 1e-8):
p_temp[:, wea_sum] /= np.sum(p_temp[:, wea_sum])
# 求结果
for wea in range(3):
for wea_sum in range(3):
p_sync_3[wea, day-2] += p_temp[wea, wea_sum] * p_sync_3[wea_sum, day-3]
p_sync_3
array([[0.8, 0. , 0. ],
[0.2, 1. , 0. ],
[0. , 0. , 1. ]])
第3问
上一问已经求出每一天的概率分布了，第3和4天都是一定的，只有第2天可能是晴天或多云，所以结果显然是晴天、多云、雨天，概率是0.8。
第四题
第1问
\[p(x) = \frac{1}{\sqrt{2\pi}\sigma_1}e^{-\frac{(x-\mu_1)^2}{2\sigma_1^2}}\\ p(z\mid x) = \frac{1}{\sqrt{2\pi}\sigma_2}e^{-\frac{(z-x)^2}{2\sigma_2^2}} \]第2问
首先求\(p(z)\)：
\[\begin{aligned} p(z) &=\int_{-\infty}^{\infty}p(z\mid x)p(x)\mathrm{d}x\\ &=\frac{1}{{2\pi\sigma_1\sigma_2}}\int_{-\infty}^{\infty}\exp\left(-\frac{(x-\mu_1)^2}{2\sigma_1^2}-\frac{(z-x)^2}{2\sigma_2^2}\right)\mathrm{d}x\\ &=\frac{1}{{6\pi\sigma_2^2}}\int_{-\infty}^{\infty}\exp\left(\frac{-(x-\mu_1)^2+9(z-x)^2}{18\sigma_2^2}\right)\mathrm{d}x\\ &=\frac{1}{{6\pi\sigma_2^2}}\int_{-\infty}^{\infty}\exp\left(\frac{-10(x-\frac{9z+\mu_1}{\sqrt{10}})^2+\frac{9}{10}(z-\mu_1)^2}{18\sigma_2^2}\right)\mathrm{d}x\\ &=\frac{1}{\sqrt{2\pi}\sigma_3}e^{-\frac{(z-\mu_1)^2}{2\sigma_3^2}}\int_{-\infty}^{\infty}\frac{1}{\sqrt{2\pi}(3\sigma_2/\sqrt{10})}\exp\left(-\frac{(x-\frac{9z+\mu_1}{\sqrt{10}})^2}{2(3\sigma_2/\sqrt{10})^2}\right)\mathrm{d}x\\ &=\frac{1}{\sqrt{2\pi}\sigma_3}e^{-\frac{(z-\mu_1)^2}{2\sigma_3^2}} \end{aligned} \]其中\(\sigma_3=\sqrt{10}\sigma_2\)。进一步化简即可得\(p(x\mid z)\)：
\[\begin{aligned} p(x\mid z) &=\frac{p(z\mid x)p(x)}{p(z)}\\ &=\frac{1}{\sqrt{2\pi}(\frac{3}{\sqrt{10}}\sigma_2)}\exp\left[\frac{(x-\frac{9z+\mu_1}{10})^2}{2\cdot(\frac{3}{\sqrt{10}}\sigma_2)^2}\right] \end{aligned} \]\(z\)已知的条件下，则就是个正态分布。
第3问
0，因为对于连续的概率分布来说观测值等于某个具体的值的概率都是0。
第五题
使用\(p(a)=p(a\mid b)p(b)\)化简2.16：
\[\begin{aligned} p(x,y\mid z) &=p(x\mid y, z)p(y\mid z)\\ &=p(y\mid x, z)p(x\mid z)\\ &=p(x\mid z)p(y\mid z) \end{aligned} \]即可得2.17和2.18。