思维导图：

19.1

用蒙特卡罗积分法求：
∫−∞∞x2exp⁡(−x22)dx\int_{-\infty}^{\infty} x^{2} \exp \left(-\frac{x^{2}}{2}\right) d x ∫−∞∞x2exp(−2x2)dx

首先将被积函数分解为分布函数与待求期望的函数的乘积：
∫−∞∞x2exp⁡(−x22)dx=2π∫−∞∞x212πexp⁡(−x22)dx=2πE[x2]\begin{aligned}&\int_{-\infty}^{\infty} x^{2} \exp \left(-\frac{x^{2}}{2}\right) d x \\ &=\sqrt{2 \pi} \int_{-\infty}^{\infty} x^{2} \frac{1}{\sqrt{2 \pi}} \exp \left(-\frac{x^{2}}{2}\right) d x \\ &=\sqrt{2 \pi} E\left[x^{2}\right]\end{aligned} ∫−∞∞x2exp(−2x2)dx=2π∫−∞∞x22π1exp(−2x2)dx=2πE[x2]
即，在正态分布下，计算x2x^{2}x2的期望值。

# 编写代码求解
import numpy as np# num_samples = 100
def Monte_Carlo_solution(num_samples): samples = np.random.normal(loc=0, scale=1, size=num_samples)result = np.sqrt(2*np.pi) * (samples**2).mean()return resultfor num_samples in [10, 100, 300, 500, 600, 800, 1000, 5000, 10000, 100000]:result = Monte_Carlo_solution(num_samples=num_samples)print(f'the result is: {result:.2f} with {num_samples} samples')

the result is: 2.35 with 10 samples
the result is: 2.85 with 100 samples
the result is: 2.30 with 300 samples
the result is: 2.54 with 500 samples
the result is: 2.81 with 600 samples
the result is: 2.46 with 800 samples
the result is: 2.46 with 1000 samples
the result is: 2.52 with 5000 samples
the result is: 2.48 with 10000 samples
the result is: 2.50 with 100000 samples

19.2

证明如果马尔可夫链是不可约的，且有一个状态是非周期的，那么其他所有的状态也是非周期的，即这个马尔可夫链是非周期的。

思想证明（反证法）：
假设存在一个状态是周期的，假设周期为T，我们假设初始时刻状态即为这个有周期的状态，接着将马尔可夫链分割，1～T，T+1～2T，2T+1～3T……，因为考虑的是时间齐次的马尔可夫链，所以，每一段中每时刻的状态都只由前一个状态确定，初始时刻加上经历时间（转移核/转移概率矩阵作用次数）就决定了该时刻的状态，从而，这几段都是从周期状态开始，经历相同时间，因此是完全等价的。再考虑不可约，如果存在某个时刻t，其转移到某个别的状态的概率一定不为0，这个t一定落在上面分段中的一段上，而每段等价，说明每段都会存在这样一个时刻，而且由于等价，其在每段内的相对位置（时刻）也是相同的，那么这些等价的时刻也会构成一个周期，那么这两另外的态也是周期的，同理，所有的态都不可约，因此都会有这么一个周期存在，从而整个链都是周期的，与假设矛盾。因此，对于不可约的马尔可夫链，如果一个状态是非周期的，其他的状态也只能都是非周期的，从而整个马尔可夫链是非周期的。

19.3

验证具有以下转移概率矩阵的马尔可夫链是可约的，但是非周期的。
P=[1/21/2001/201/2001/200001/21]P=\left[\begin{array}{cccc}1 / 2 & 1 / 2 & 0 & 0 \\ 1 / 2 & 0 & 1 / 2 & 0 \\ 0 & 1 / 2 & 0 & 0 \\ 0 & 0 & 1 / 2 & 1\end{array}\right] P=⎣⎢⎢⎡1/21/2001/201/2001/201/20001⎦⎥⎥⎤

# 直接代码验证，比较直观
trans_prob = np.array([[0.5, 0.5, 0, 0],[0.5, 0, 0.5, 0], [0, 0.5, 0, 0],[0, 0, 0.5, 1]])
n = 30
Markov_sequence = []
start_index = np.random.randint(0, 4)
state = np.zeros(4).reshape(-1, 1)
state[start_index] = 1for i in range(n):Markov_sequence.append(state.argmax())state = trans_prob @ statestate = state + np.random.rand(4).reshape(-1, 1) * 1e-5  # 使得概率相等时，能够随机选取一个状态state_idx = state.argmax()state = np.zeros(4).reshape(-1, 1)state[state_idx] = 1print(Markov_sequence)

[1, 2, 1, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3]

可以发现，一旦进入第4个状态，那么就一直在这个状态了，上面例子中第1状态永远无法得到，所以是可约的，如果开始就是第4状态，那么，1、2、3状态都无法转到;
上面例子也可以看到，没有周期性，是非周期的。

19.4

验证具有以下转移概率矩阵的马尔可夫链是不可约的，但是周期性的。
P=[01/200101/2001/201001/20]P=\left[\begin{array}{cccc}0 & 1 / 2 & 0 & 0 \\ 1 & 0 & 1 / 2 & 0 \\ 0 & 1 / 2 & 0 & 1 \\ 0 & 0 & 1 / 2 & 0\end{array}\right] P=⎣⎢⎢⎡01001/201/2001/201/20010⎦⎥⎥⎤

# 直接代码验证，比较直观
trans_prob = np.array([[0, 0.5, 0, 0],[1, 0, 0.5, 0], [0, 0.5, 0, 1],[0, 0, 0.5, 0]])
n = 500
Markov_sequence = []
start_index = np.random.randint(0, 4)
state = np.zeros(4).reshape(-1, 1)
state[start_index] = 1for i in range(n):Markov_sequence.append(state.argmax())state = trans_prob @ statestate = state + np.random.rand(4).reshape(-1, 1) * 1e-5  # 使得概率相等时，能够随机选取一个状态state_idx = state.argmax()state = np.zeros(4).reshape(-1, 1)state[state_idx] = 1print(Markov_sequence)

[0, 1, 0, 1, 0, 1, 2, 1, 0, 1, 0, 1, 2, 3, 2, 3, 2, 3, 2, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 2, 1, 0, 1, 0, 1, 0, 1, 2, 1, 2, 1, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 1, 2, 3, 2, 1, 2, 1, 2, 1, 0, 1, 2, 3, 2, 1, 0, 1, 2, 3, 2, 1, 2, 3, 2, 1, 2, 1, 2, 1, 2, 3, 2, 1, 0, 1, 0, 1, 2, 1, 0, 1, 0, 1, 2, 3, 2, 3, 2, 3, 2, 1, 0, 1, 0, 1, 2, 1, 0, 1, 2, 3, 2, 1, 2, 3, 2, 1, 2, 1, 0, 1, 2, 1, 2, 3, 2, 1, 0, 1, 0, 1, 2, 3, 2, 3, 2, 1, 0, 1, 2, 1, 0, 1, 2, 1, 2, 3, 2, 1, 2, 3, 2, 1, 0, 1, 2, 1, 0, 1, 0, 1, 0, 1, 0, 1, 2, 3, 2, 3, 2, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 2, 1, 0, 1, 0, 1, 0, 1, 2, 1, 2, 1, 2, 3, 2, 3, 2, 3, 2, 3, 2, 1, 0, 1, 2, 1, 0, 1, 0, 1, 2, 1, 0, 1, 0, 1, 0, 1, 2, 1, 2, 3, 2, 3, 2, 1, 2, 3, 2, 3, 2, 1, 2, 1, 2, 1, 0, 1, 2, 3, 2, 1, 0, 1, 2, 1, 0, 1, 2, 3, 2, 1, 0, 1, 2, 3, 2, 3, 2, 1, 2, 1, 2, 3, 2, 3, 2, 3, 2, 3, 2, 1, 0, 1, 2, 3, 2, 3, 2, 1, 2, 1, 0, 1, 0, 1, 2, 1, 0, 1, 2, 3, 2, 3, 2, 1, 2, 3, 2, 3, 2, 1, 0, 1, 0, 1, 2, 1, 0, 1, 0, 1, 2, 3, 2, 1, 2, 1, 2, 3, 2, 3, 2, 1, 0, 1, 2, 1, 0, 1, 0, 1, 2, 1, 0, 1, 0, 1, 0, 1, 0, 1, 2, 1, 2, 3, 2, 1, 2, 3, 2, 3, 2, 3, 2, 3, 2, 1, 2, 1, 2, 3, 2, 1, 2, 1, 0, 1, 2, 3, 2, 1, 2, 3, 2, 1, 2, 1, 2, 3, 2, 3, 2, 3, 2, 1, 2, 3, 2, 3, 2, 1, 2, 1, 0, 1, 0, 1, 2, 3, 2, 3, 2, 3, 2, 1, 2, 3, 2, 3, 2, 3, 2, 1, 2, 1, 2, 1, 2, 1, 0, 1, 2, 1, 2, 3, 2, 3, 2, 1, 2, 1, 2, 3, 2, 3, 2, 1, 0, 1, 2, 3, 2, 3, 2, 3, 2, 1, 0, 1, 0, 1, 0, 1, 2, 1, 0, 1, 2, 3, 2, 1, 2, 3, 2, 1, 2, 3, 2, 3, 2, 1, 2, 3, 2, 1, 0, 1, 2, 3, 2, 1, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3]

可以发现0, 1, 2, 3所代表的第1，2，3，4状态都可以出现，因此是不可约的。接下来可以下周期，注意到：

(np.nonzero(aaa == 1)[0] - 1) % 2

array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0, 0, 0, 0, 0, 0, 0])

其中小括号内部根据周期性的定义，将所要研究的状态作为0时刻出发的态，我们把所有从第二状态出发（0时刻），到达第二状态的时刻找出来，发现其公约数中有2（并非非周期只有1），因此是周期的，根据19.2的证明过程，不可约的马尔可夫链中，有一个态是周期的，那么所有态都是周期的，且周期是相同的，下面验证一下。

((np.nonzero(aaa == 0)[0]) % 2).any()

False

((np.nonzero(aaa == 2)[0] - 6) % 2).any()

False

((np.nonzero(aaa == 3)[0] - 13) % 2).any()

False

得证。

19.5

证明可逆马尔可夫链一定是不可约的。

简单思想证明：对于可逆马尔可夫链，其性质决定了，如果从平稳分布起始，如果一个分布在未来能够到达，那么过去也能到达，如果可约那么过去和未来将都无法到达，只能是孤立的，而孤立的就不在马尔可夫链上了。

19.6

从一般的Metropolis-Hastings算法推导出单分量Metropolis-Hastings算法。

好像不是很难理解吧，需要推导吗？

19.7

假设进行伯努利实验，后验概率为P(θ∣y)P(\theta |y)P(θ∣y)，其中变量y∈{0,1}y \in \{0,1\}y∈{0,1}表示实验可能的结果，变量θ\thetaθ表示结果为1的概率。再假设先验概率P(θ)P(\theta)P(θ)遵循Beta分布B(α,β)B(\alpha, \beta)B(α,β)，其中α=1,β=1\alpha=1,\beta=1α=1,β=1：似然函数P(y∣θ)P(y |\theta)P(y∣θ)遵循二项分布Bin(n,k,θ)Bin(n,k,\theta)Bin(n,k,θ)，其中n=10,k=4n=10,k=4n=10,k=4，即实验进行10次其中结果为1的次数为4。试用Metropolis-Hastings算法求后验概率分布P(θ∣y)∝P(θ)P(y∣θ)P(\theta |y)\propto P(\theta)P(y |\theta)P(θ∣y)∝P(θ)P(y∣θ)的均值和方差。（提示：可以采用Metropolis选择，即假设建议分布是对称的。）

import numpy as np
# import math
# 采用Metropolis选择，并使用书中给出的与两个随机变量绝对值相关的建议分布的特例，即一个以原状态为中心的高斯分布# 假设初始的概率为0.5
theta0 = 0.5
# 选定迭代次数
n = 20000
# 执行循环
theta = theta0
theta_list = [theta0, ]
for i in range(n):theta_old = thetatheta = np.random.normal(loc=theta_old, scale=0.2, size=1)[0]# 限制取值范围if theta < 0:theta = 0elif theta > 1:theta = 1acceptor = min(1,(theta**4 * (1 - theta)**6) / (theta_old**4 * (1 - theta_old)**6))  # 注意先验分布在给定条件下为均匀分布random_num = np.random.rand(1)if random_num > acceptor:theta = theta_oldtheta_list.append(theta)
# 选定认为的达到平稳分布的时刻m
m = 2000
# 计算待求的函数
expectation = np.array(theta_list[m:]).mean()
variation = np.array(theta_list[m:]).var()
# 打印结果
print(f'iteration {n} with {m} to be stable, exectation:{expectation}, variation:{variation}')

iteration 20000 with 2000 to be stable, exectation:0.4146535483495147, variation:0.0188274672691176

可以重复多次试验m，n的选择，从而找到收敛的结果，上面超参数的情况下，基本上期望值在0.41～0.42，方差基本为0.019

19.8

设谋试验可能有五种结果，其出现的概率分别为：
θ4+18,θ4,η4,η4+38,12(1−θ−η)\frac{\theta}{4}+\frac{1}{8}, \quad \frac{\theta}{4}, \quad \frac{\eta}{4}, \quad \frac{\eta}{4}+\frac{3}{8}, \quad \frac{1}{2}(1-\theta-\eta) 4θ+81,4θ,4η,4η+83,21(1−θ−η)
模型含有两个参数θ\thetaθ和η\etaη，都介于0和1之间。现有22次试验结果的观测值为：
y=(y1,y2,y3,y4,y5)=(14,1,1,1,5)y=\left(y_{1}, y_{2}, y_{3}, y_{4}, y_{5}\right)=(14,1,1,1,5) y=(y1,y2,y3,y4,y5)=(14,1,1,1,5)
其中yiy_{i}yi表示22次试验中第iii个结果出现的次数，i=1,2,...,5i=1,2,...,5i=1,2,...,5。试用吉布斯抽样估计参数θ\thetaθ和η\etaη的均值和方差。

这个问题相当两个分量θ1,θ2\theta_{1},\theta_{2}θ1,θ2（这里是θ,η\theta,\etaθ,η）的多变量的取样，即看作有两个分量的多变量，从而应用“单分量”吉布斯抽样算法。抽取时，按照p(θ∣η,y)p(\theta | \eta, y)p(θ∣η,y)对θ\thetaθ分量抽样，按照p(η∣θ,y)p(\eta | \theta, y)p(η∣θ,y)对η\etaη分量抽样。当η\etaη已知时，根据已经有的22次的试验结果，θ\thetaθ有三种可取的值，这里认为其条件分布为这三个可取值的均匀分布，η\etaη的情况与此类似。

# 定义一个函数限制参数的范围
def normalization(x):x = max(x, 0)x = min(x, 1)return x

# 初始化
theta, eta = 0.5, 0.5
# 选定迭代次数
n = 50000
# 执行迭代循环
sample_list = [(theta, eta), ]
for i in range(n):theta = (1/26, (51 - 56 * eta)/66, (2 - 2 * eta)/7)[np.random.randint(0, 3)]theta = normalization(theta)eta = (0, (2 - 2 * theta)/7, (-11 - 4 * theta)/14)[np.random.randint(0, 3)]eta = normalization(eta)sample_list.append((theta, eta))
# 选定认为的达到平稳分布的时刻m
m = 5000
# 计算待求的函数
theta_list = [theta for theta in zip(*sample_list[m:])][0]
eta_list = [eta for eta in zip(*sample_list[m:])][1]
theta_expectation = np.array(theta_list).mean()
theta_variation = np.array(theta_list).var()
eta_expectation = np.array(eta_list).mean()
eta_variation = np.array(eta_list).var()
# 打印结果
print(f'iteration {n} with {m} to be stable')
print(f'theta_exectation:{theta_expectation}, theta_variation:{theta_variation}')
print(f'eta_exectation:{eta_expectation}, eta_variation:{eta_variation}')

iteration 50000 with 5000 to be stable
theta_exectation:0.3400718094477886, theta_variation:0.08291996043361917
eta_exectation:0.06193353766656458, eta_variation:0.0100438612941617

仍然可以尝试多次确定m，n，最终θ\thetaθ的期望大约为0.34，η\etaη期望大约为0.062。方差分别是0.08，0.01。将期望作为其取值，那么得到的五个概率依此为:

a = theta_variation
b = eta_expectation
a/4 + 1/8, a/4, b/4, b/4 + 3/8, 0.5*(1-a-b)

(0.1457299901084048,0.020729990108404792,0.015483384416641145,0.39048338441664115,0.4275732509499081)

概率之和一定是1，因为原来的已知条件的加和与两个参数无关，这个题目说明了，利用有限次（不能完成大数极限）的试验，可以借助马尔可夫蒙特卡罗的吉布斯抽样算法估计参数。

《统计学习方法（第2版）》李航第19章马尔可夫蒙特卡罗法 MCMC 思维导图笔记及课后全部习题答案（步骤详细，包含Metropolis算法，吉布斯算法代码实现）第十九章相关推荐

《统计学习方法（第2版）》李航第20章潜在狄利克雷分配 LDA Dirichlet 思维导图笔记及课后全部习题答案（步骤详细，包含吉布斯抽样算法）狄利克雷分布期望推导
思维导图: 20.1 推导狄利克雷分布数学期望公式. 首先写出Dirichlet分布的概率密度函数: ρ(θ)=Γ(α0)Γ(α1)⋯P(αn)∏i=1nθiαi−1\rho(\theta)=\fra ...
python做马尔科夫模型预测法_李航《统计学习方法》第十章——用Python实现隐马尔科夫模型...
相关文章: 李航<统计学习方法>第二章--用Python实现感知器模型(MNIST数据集) 李航<统计学习方法>第三章--用Python实现KNN算法(MNIST数据集) 李航 ...
思维导图下载注册安全_2019安全工程师《建筑实务》第二章第一节考点：物料提升机思维导图...
2019安全工程师<安全生产专业实务-建筑施工安全>第二章第一节考点:物料提升机思维导图,本节的大部分知识点前两节塔式起重机和施工升降机相似,大家可以对比之前考点的思维导图来理解记忆,本知 ...
李航书上隐马尔科夫模型案例的实验结果复现
主要是转载自: https://www.cnblogs.com/pinard/p/7001397.html 我自己加了些注释以及修复了一个在python2.7下面运行的bug 1) 评估观察序列概率. ...
《Python编程从入门到实践》记录之第2章变量和简单数据类型总结（思维导图）
<Python编程从入门到实践>第2章变量和简单数据类型知识总结:
《机器学习》周志华 --第3章线性模型思维导图+笔记+习题
基本形式问题描述:给定由d个属性描述的示例x=(x1;x2;x3-xd),其中xi是x的第i个属性上的取值,线性模型试图学得一个通过属性的线性组合来进行预测函数, 函数形式:f(x) = w1x1+ ...
自动控制原理第1章——自动控制系统的基本概念（思维导图）
第1章--自动控制系统的基本概念闭环控制系统结构图: 内容参考: 王建辉.顾树生老师所主编的自动控制原理(第2版) 程鹏老师所主编的自动控制原理(第2版) 如有错误或者不足之处,欢迎大家留言指正!
计算机网络第五章（谢希仁）--运输层思维导图
《机器学习》周志华第10章降维与度量学习思维导图+笔记+习题
K-Means与LVQ都试图以类簇中心作为原型指导聚类,其中K-Means通过EM算法不断迭代直至收敛,LVQ使用真实类标辅助聚类:高斯混合聚类采用高斯分布来描述类簇原型:密度聚类则是将一个核心对象所 ...

《统计学习方法（第2版）》李航第19章马尔可夫蒙特卡罗法 MCMC 思维导图笔记及课后全部习题答案（步骤详细，包含Metropolis算法，吉布斯算法代码实现）第十九章

19.1

19.2

19.3

19.4

19.5

19.6

19.7

19.8

《统计学习方法（第2版）》李航第19章马尔可夫蒙特卡罗法 MCMC 思维导图笔记及课后全部习题答案（步骤详细，包含Metropolis算法，吉布斯算法代码实现）第十九章相关推荐

最新文章

热门文章

《统计学习方法（第2版）》李航 第19章 马尔可夫蒙特卡罗法 MCMC 思维导图笔记 及 课后全部习题答案（步骤详细， 包含Metropolis算法，吉布斯算法代码实现）第十九章

19.1

19.2

19.3

19.4

19.5

19.6

19.7

19.8

《统计学习方法（第2版）》李航 第19章 马尔可夫蒙特卡罗法 MCMC 思维导图笔记 及 课后全部习题答案（步骤详细， 包含Metropolis算法，吉布斯算法代码实现）第十九章相关推荐

最新文章

热门文章

《统计学习方法（第2版）》李航第19章马尔可夫蒙特卡罗法 MCMC 思维导图笔记及课后全部习题答案（步骤详细，包含Metropolis算法，吉布斯算法代码实现）第十九章

《统计学习方法（第2版）》李航第19章马尔可夫蒙特卡罗法 MCMC 思维导图笔记及课后全部习题答案（步骤详细，包含Metropolis算法，吉布斯算法代码实现）第十九章相关推荐