martingale、markov chain、Monte Carlo、MCMC

文章结构如下：

1: MCMC

1.1 MCMC是什么
1.2 为什么需要MCMC
2： 蒙特卡罗

2.1 引入
2.2 均匀分布，Box-Muller 变换
2.3 拒绝接受采样（Acceptance-Rejection Sampling）
2.4 接受拒绝采样的直观解释
2.5 接受拒绝采样方法有效性证明
2.6 接受拒绝采样方法python实现
2.7 蒙特卡罗方法小结
3： 马尔科夫链

3.1 马尔科夫链概述
3.2 马尔科夫链模型状态转移矩阵的性质
3.3 基于马尔科夫链采样
3.4 马尔科夫链采样小结
3.5 补充：PageRank：Google的民主表决式网页排名技术
4： 马尔可夫链蒙特卡罗算法

4.1 马尔科夫链的细致平稳条件(Detailed Balance Condition)
4.2 MCMC采样
5： M-H采样

5.1 M-H采样算法
5.2 M-H采样python实现
5.3 M-H采样小结
6：Gibbs采样

6.1 重新寻找合适的细致平稳条件
6.2 二维Gibbs采样
6.3 多维Gibbs采样
6.4 二维Gibbs采样python实现
6.5 Gibbs采样小结
7： 参考文献

一、基本知识

1. 条件概率

条件概率是指在某件事情已经发生的前提下，另一件事情在此基础上发生的概率，举例来说P(A丨B)表示B发生的基础上，A也发生的概率，基本公式为：

2. 条件期望

在上述概率下的期望我们称之为条件期望。
计算：

3. 随机过程

在概率论概念中，随机过程是随机变量的集合。若一随机系统的样本点是随机函数，则称此函数为样本函数，这一随机系统全部样本函数的集合是一个随机过程。实际应用中，样本函数的一般定义在时间域或者空间域。

顾名思义，它其实就是个过程，比如今天下雨，那么明天下不下雨呢？后天下不下雨呢？从今天下雨到明天不下雨再到后天下雨，这就是个过程。那么怎么预测N天后到底下不下雨呢？这其实是可以利用公式进行计算的，随机过程就是这样一个工具，把整个过程进行量化处理，用公式就可以推导出来N天后的天气状况，下雨的概率是多少，不下雨的概率是多少。说白了，随机过程就是一些统计模型，利用这些统计模型可以对自然界的一些事物进行预测和处理，比如天气预报，比如股票，比如市场分析，比如人工智能。

随机过程是一族时间函数的集合，随机过程的每个样本函数（一个随机信号）是一个确定的时间函数x(t)（t=0~T），随机过程在一个确定的时刻t1是一个随机变量X(t1)。
从上述对随机过程描述中可以看出，从信号分析的角度来看，一个随机过程就是针对一个实验测量无数个的样本，每个样本就是一个随机信号（要分清什么是随机过程和随机信号），而所有样本测量的时间都是0~T（这看起来很麻烦和不可思议），然后，我们把这些样本（每个随机信号）放在一起摆放，见图

其实，这就是一个随机过程，然后，我们怎样描述呢？我们从两个角度来描述它（和上述定义对应）。
（1）从图中纵向的角度（0~T），一个随机过程是你这次针对于某个任务进行测量的所有样本信号的集合。
（2）从横向的角度（在0-T区间的任意时刻），你可以拿把“菜刀”在任意时刻横着切一下，这时，你会看到所有你测量的样本信号在这个时刻的幅度值都不一样，它们是随机的，为了表征这些幅度值的统计特性，我们用此时刻的随机变量X(t)来表征，那么在0~T时间区间内有多少随机变量呢？当然是无数个（相当于你切了无数刀），这些无数个随机变量的集合就是随机过程。
如果所有工作都像这样去做的话，那么，所有随机信号分析就没有办法完成了，所以，人们将问题进行简化，尽量将随机过程说成平稳的并且具有各态历经的。

二、鞅（martingale）

在概率论中，鞅（martingale）是满足下述条件的随机过程：已知过去某一 时刻 s 以及之前所有时刻的观测值，若某一时刻 t 的观测值的条件期望等于过去某一时刻 s 的观测值，则称这一随机过程是鞅。而于博弈论中，鞅经常用来作为公平博弈的数学模型。

鞅的原名 martingale 原指一类于18世纪流行于法国的投注策略，称为加倍赌注法[1]。这类策略中最简单的一种策略是为博弈设计的。在博弈中，赌徒会掷硬币，若硬币正面向上，赌徒会赢得赌本，若硬币反面向上，赌徒会输掉赌本。这一策略使赌徒在输钱后加倍赌金投注，为的是在初次赢钱时赢回之前输掉的所有钱，同时又能另外赢得与最初赌本等值的收益。当赌徒的财产和可用时间同时接近无穷时，他掷硬币后硬币正面向上的概率会接近1，由此看来，加倍赌注法似乎是一种必然能赢钱的策略。然而，赌金的指数增长最终会导致使用这一策略的赌徒破产。

鞅的概念首先是由保罗·皮埃尔·莱维于 1934 年提出的，但他只提出了离散时间的版本，而且没有给予命名。直到 1939 年，约翰・维尔将此概念推广到连续时间的情况，并且首次提出 martingale 这个名称。约瑟夫·利奥·杜布（Joseph Leo Doob）等人在鞅的相关理论的初期发展做出重大贡献，而完成这些工作的部分动机是为了表明成功的投注策略不可能存在。此外，伊藤清在分析应用方面作出了重要的贡献。从 1970 年代开始，鞅论就在纯粹数学和应用数学的很多领域中有广泛的应用，特别是在数学物理和金融数学中。

以下我们讨论T个离散阶段下，鞅的相关定义和性质。

我们将一个随机试验的所有基本结果称为样本点，记为：
称其构成的集合文件样本空间：
其对应的概率空间记为：
记每个阶段的状态为S0，S1, …, ST，我们定义以下集合：
Ft可以理解为t时刻前所有状态的信息，即整个随机过程到达当前状态的所有可能路径。我们再定义:

这样一来F就包含了所有时刻，所有状态，所有路径的信息。以下的性质显而易见：

因为现t2时刻的信息肯定是要多于t1时刻的。
我们继续有以下定义：

接下来我们给出条件期望的基本性质：

2.1 Martingale的定义

我们开始讨论鞅的定义，鞅有以下三种等价的定义：

我们以下证明三个定义是等价的：
以上便是鞅的基本定义，为了更好的理解鞅，我们可以举个例子。
假设赌徒在t时刻有$M_t$ 的财富，如果赌局是公平的，那么他未来财富的期望应当是与当前财富相同，同时不受之前赌博策略的影响（因为我们求的是条件期望）。事实上martingale理论最早便是用于赌博策略的研究，我们现在讨论的鞅是公平赌局，也存在像上鞅和下鞅这样不公平赌局。至于为什么中文会翻译成鞅呢，因为在法文里，martingale有两层意思，一是“倍赌策略”（每输一局就赌资加倍，只要赢一局就扳回所有损失），二是“马缰绳”，将其翻译成鞅（马的缰绳）是用了它的第二层含义。

2.2 Martingale, supermartingale, submartingale

鞅：E[Xt | Fs] = Xs （t>s）；
上鞅：E[Xt | Fs] < Xs （t>s）；
下鞅：E[Xt | Fs] > Xs （t>s）。

通俗理解：Fs表示你在s时刻所掌握的所有信息，那么E[Xt | Fs]就表示你在s时刻所掌握的所有信息对X在未来t时刻预期的判断，如何和X在s时刻是一致的说明游戏是公平的；上鞅说明对过程X是看衰的，下鞅则是看好。
鞅是公平赌博，下鞅是有利赌博，上鞅是不利赌博。

鞅直观意义是一条（有固定趋势的）线。
考虑一个函数f(t)，横轴是时间轴，纵轴是函数值。这个函数可以是平的（就是$f(t)=c$），也可以单调增，也可以单调减。因此从图像上看是一条有固定趋势的线。
如果仅仅是线，那么未免太乏味无趣。概率就是把固定取值的变量$X=c$取代为可以同时取多个值、每个值都有一定权重的随机变量$X$，相当于引入了多个可能状态。随机过程进一步引入一个时间轴，考虑随机变量的序列$X_n$（或者连续版本的$X_t$），因此如果画在刚刚那个图里，$X_n(w)$可以看做某个大体上是线、局部有些抖动的点列，而$X_n$就是这些点列的集合。
那么离散鞅是一种点列集（随机过程），这个点列集大体趋势是$f(t)=c$，只是会有一些抖动，而且是考虑一个抖动围绕这条线的点列集整体。类似的，上鞅、下鞅就是单调减、单调增函数的点列集。

2.3 Stopping time

关于随机过程 X1，X2，X3，… 的停时是随机变量 τ，这一随机变量具有如下性质：对于每一个时间 ，事件 τ = t 的发生与否仅取决于 X1，X2，X3，…，Xt 的取值。从定义中可以感受到的直觉是在任一特定时刻 t，我们都可以知道在这一时刻随机过程是否到了停时。现实生活中停时的例子如赌徒离开赌桌的时刻，这一时刻可能是赌徒以前赢得钱财的函数（例如，仅当他没有钱时，他才可能离开赌桌），但是他不可能根据还未完成的博弈的结果来选择离开还是留下。

上述停时定义满足强条件，下面给出一个弱条件的停时定义：若事件 τ = t 的发生与否统计独立于 Xt+1，Xt+2，… 但并不是完全决定于时刻 t 以及之前的过程历史，则随机变量 τ 是停时。虽然这是一个弱条件，但在需要用到停时的证明中的一些情况也算是足够强的条件。

鞅的一个基本性质是若$(X_t{)}_{t>0}$是下\上鞅且 $\tau$ 是停时，由 X t τ : = X min ⁡ { τ , t } X_{t}^{\tau }:=X_{
{\min\{\tau ,t\}}} Xtτ​:=Xmin{
τ,t}​定义的对应停止过程$(X_t^{\tau }{)}_{t>0}$也是下\上鞅。

停时鞅的概念引出了一系列定理，例如可选停止定理（又称可选抽样定理）：在特定条件下，停时的鞅的期望值等于其初始值。利用这一定理，我们可以证明对于一个寿命有限且房产有限的赌徒，成功的投注策略不可能存在。

2.3.1 Optional stopping theorem

In probability theory, the optional stopping theorem (or Doob’s optional sampling theorem) says that, under certain conditions, the expected value of a martingale at a stopping time is equal to its initial expected value. Since martingales can be used to model the wealth of a gambler participating in a fair game, the optional stopping theorem says that, on average, nothing can be gained by stopping play based on the information obtainable so far (i.e., without looking into the future). Certain conditions are necessary for this result to hold true. In particular, the theorem applies to doubling strategies.

The optional stopping theorem is an important tool of mathematical finance in the context of the fundamental theorem of asset pricing.

2.3.1.1 discrete-time version of the theorem

2.3.1.2 Applications

2.3.1.3 Proof

三、markov chain

3.1 Markov Property

在学习马尔可夫链之前，我们首先了解一下什么是马尔可夫性质（或者叫马尔可夫特性，Markov Property）。
假设你有一个系统，这个系统具有$M$ 种可能的状态，并且在状态之间正在移动。真实世界中，我们常常见到这种例子，例如天气从炎热到温和，或者是股票市场从熊市到牛市再到停滞的状态（stagnant states）。

3.1.1 马尔可夫性的定义

关于马尔可夫性，我们给出了如下的Definition：

从上述的式子可以看出，t+1时刻的状态包含了1,…,t时刻状态的全部历史信息，并且当我们知道t时刻的状态后，我们只关注于环境的信息，而不用管之前所有状态的信息，这就是马尔可夫性，当论文中说某一状态或其他信息符合马尔可夫性时，我们也应当联想到这个性质。

3.1.2 State Transition Matrix状态传输矩阵

如果一个过程具有马尔可夫性质，它也被成为马尔可夫过程（Markov Process）。

3.2 Markov Processes马尔可夫过程

在概率论及统计学中，马尔可夫过程（英语：Markov process）是一个具备了马尔可夫性质的随机过程，因为俄国数学家安德雷·马尔可夫得名。马尔可夫过程是不具备记忆特质的（memorylessness）。换言之，马尔可夫过程的条件概率仅仅与系统的当前状态相关，而与它的过去历史或未来状态，都是独立、不相关的。

3.3 markov chain马尔可夫链

具备离散状态的马尔可夫过程，通常被称为马尔可夫链。

3.3.1 定义与解释

3.3.1 定义举例1

3.3.2 定义举例2

Markov Chain 体现的是状态空间的转换关系，下一个状态只决定于当前的状态.通俗点说就是：你下一时刻在哪只跟你现在在哪有关，于你以前去过哪没有关系. 如下图：

这个状态图的转换关系可以用一个转换矩阵 T 来表示：

举一个例子，如果当前状态为 u(x) = (0.5, 0.2, 0.3), 那么下一个矩阵的状态就是 u(x)T = (0.18, 0.64, 0.18), 依照这个转换矩阵一直转换下去，最后的系统就趋近于一个稳定状态 (0.22, 0.41, 0.37) 。而事实证明无论你从那个点出发，经过很长的 Markov Chain 之后都会汇集到这一点。

满足什么条件下经过很长的 Markov Chain 后系统会趋近一个稳定状态呢，大概的条件如下：

• Irreducibility. 即图是联通的，各个状态之间都有过去的办法，举个不联通的例子，比如爬虫爬不到内部局域网的网页
• Aperiodicity. 非周期性, 即图中遍历不会陷入到一个死圈里，进去了再也出不来，有些网站为了防机器人，会专门设置这种陷阱
• Detailed Balance，这是保证系统有稳态的一个重要条件，详细说明见下面。
  假设 p(x) 是最后的稳态，那么 detailed balance 可以用公式表示为：
  
  什么意思呢？假设上面状态图 x1 有 0.22 元， x2 有 0.41 元，x3 有 0.37 元，那么 0.22×1 表示 x1 需要给 x2 钱，以此类推，手动计算，可以发现下一个状态每个人手中的钱都没有变。值得说明的是，这里体现了一个很重要的特性，那就是从一个高概率状态 xi 向一个低概率状态 x(i-1) 转移的概率等于从这个低概率状态向高概率状态转移的概率（reversible，至于要不要转移又是另外一回事)
• ergodicity, 遍历性，是指不管事物现在出于什么状态，在较长时间内，马尔科夫过程逐渐趋于稳定状况，而且稳定状态与初始状况无关。

3.3.3 通俗解释

先说说我们村智商为0的王二狗，人傻不拉几的，见人就傻笑，每天中午12点的标配，仨状态：吃，玩，睡。这就是传说中的状态分布。

你想知道他n天后中午12点的状态么？是在吃，还是在玩，还是在睡？这些状态发生的概率分别都是多少？
先看个假设，他每个状态的转移都是有概率的，比如今天玩，明天睡的概率是几，今天玩，明天也玩的概率是几几，看图更清楚一点。

这个矩阵就是转移概率矩阵P，并且它是保持不变的，就是说第一天到第二天的转移概率矩阵跟第二天到第三天的转移概率矩阵是一样的。（这个叫时齐，有兴趣的同学自行百度）。
有了这个矩阵，再加上已知的第一天的状态分布，就可以计算出第N天的状态分布了。

S1 是4月1号中午12点的的状态分布矩阵 [0.6, 0.2, 0.2]，里面的数字分别代表吃的概率，玩的概率，睡的概率。
那么
4月2号的状态分布矩阵 S2 = S1 * P (俩矩阵相乘)。
4月3号的状态分布矩阵 S3 = S2 * P (看见没，跟S1无关，只跟S2有关)。
4月4号的状态分布矩阵 S4 = S3 * P (看见没，跟S1，S2无关，只跟S3有关)。
…
4月n号的状态分布矩阵 Sn = Sn-1 * P (看见没，只跟它前面一个状态Sn-1有关)。

总结：马尔可夫链就是这样一个任性的过程，它将来的状态分布只取决于现在，跟过去无关！

就把下面这幅图想象成是一个马尔可夫链吧。实际上就是一个随机变量随时间按照Markov性进行变化的过程。

附：S2 的计算过程 (没兴趣的同学自行略过)

3.3.2 n步转移概率

除了一步转移以外，马尔科夫链还有n步转移，即通过n次达到目标状态

其实是Markov Chain 的基本性质“无后效性”，即事物将来的状态及其出现的概率的大小，只取决于该事物现在所处的状态，而与以前时间的状态无关

3.3.2.1 C-K方程

查普曼-柯尔莫格洛夫方程(Chapman-Kolmogorov equation,C-K equation)给出了计算 [公式] 步转移概率的一个方法：

3.3.3 状态分类

①可达的：i到j的概率为正

②常返的：从i出发后还能再回到i

③非常返的：从i出发后不一定能回到j

常返和非常返的区别！：常返指出发后可以无限次回到自身；非常返并不是无法回到自身，而是只能有限次回到自身。

常返类：相互可达的状态组。

马尔科夫链的分解：

• 一个马尔科夫链可以分解成一个或者多个常返类，加上可能的非常返类。
• 一个常返态从它所属的类里任何一个状态出发都是可达的，但从其他类里的常返状态出法是不可达的。
• 从任何一个常返状态出发都不可到达非常返状态
• 从一个非常返状态出发，至少有一个常返态是可达的

常返类的周期：

若一个常返类能被分成几个组，这几个组之间的下一个状态转移必定不在自己组。则这个常返类被称为有周期的。

命题得证。

很显然，暂态也是一个类性质。而利用上述性质可以得到：有限马尔可夫链的所有状态不可能都是暂态，有限不可约马尔可夫链的所有状态都是常返态。

3.3.4 稳态

3.3.4.1 稳态概率举例1

3.3.4.2 稳态概率举例2

举个具体的例子。社会学家把人按其经济状况分为3类：下层，中层，上层，我们用1,2,3表示这三个阶层。社会学家发现决定一个人的收入阶层最重要的因素就是其父母的收入阶层。如果一个人的收入属于下层类别，则它的孩子属于下层收入的概率为0.65，属于中层收入的概率为0.28，属于上层收入的概率为0.07。从父代到子代，收入阶层转移概率如下

我们用P表示这个转移矩阵，则

假设第1代人的阶层比例为
则前10代人的阶层分布如下

我们可以看到，在相同的转移矩阵作用下，状态变化最终会趋于平稳。对于第n代人的阶层分布，我们有
从表达式上我们可以看到，π是一维向量，P是两维矩阵，P进行足够多次自乘后，值趋于稳定。

3.3.5 马氏链平稳分布

在转移矩阵P作用下达到的平稳状态，我们称之为马氏链平稳分布。对于这个特性，有如下精彩定理

我在这里直观的解释一下上面定理

条件

（1）非周期马氏链：马氏链转移要收敛，就一定不能是周期性的。不做特别处理，我们处理的问题基本上都是非周期性的，在此不做多余解释。
（2）存在概率转移矩阵P，任意两个状态是连通的：这里的连通可以不是直接相连，只要能够通过有限次转移到达即可。比如对于a, b, c状态，存在a->b, b->c，则我们认为a到c是可达的。

结论

（1）不论初始状态是什么，经过足够多次概率转移后，会存在一个稳定的状态π。
（2）概率转移矩阵自乘足够多次后，每行值相等。即

3.3.5.1 马尔科夫链平稳状态定理的物理解释

我们再用一个更加简单的例子来阐明这个定理的物理含义。假设城市化进程中，农村人转移为城市人的概率为0.5，城市人转移为农村人的概率为0.1。

假设一开始有100个农村人，0个城市人，每代转移人数如下

可以看到，城市化进程中马尔科夫平稳状态就是农村人转移为城市人的速度等于城市人转移为农村人的速度。对于上述转移矩阵P，平稳分布为农村人17%，城市人83%。如果我们可以得到当前中国城市化转移矩阵P，我们就可以算出中国最终城市化率大概为多少（这里不考虑P的变化）。同时如果我们知道了中国城市化人口比例，我们就能知道城市化进程还可以持续多少代人。

3.3.5.2 长期频率解释

3.3.6 Birth-Death process 生灭过程

3.3.6.1 平稳分布

3.3.6.2 性质

3.3.6.3 例子

3.3.6.4 Probability balance equations 概率平衡方程

3.3.7 吸收概率和吸收期望时间

3.3.7.1 平均时间

平均吸收时间：i的平均吸收时间为从状态i开始直到到达吸收态所需要的期望步数

平均吸收时间计算：解平均吸收时间方程组

3.3.8 为什么马尔可夫链这么重要

因为它的静态分布（或叫平稳性分布，stationary distribution）性质。
我们用一下Wiki上给出的股票市场的例子（见下图），这是一个具有静态空间的马尔可夫链（a continuous-time Markov chain with state-space），包含了 牛市、熊市、停滞市场三种状态（Bull market, Bear market, Stagnant market） 。

它各状态之间由转移概率组成的转移矩阵（ transition rate matrix）如下，记为 Q：


你可以尝试从其他状态出发，但最终的静态分布一定是相同的。那么得出这个静态分布有什么用处呢？

静止状态分布的重要性在于它可以让你在随机的时间里，对于一个系统的任意一个状态确定其概率。（The stationary state distribution is important because it lets you define the probability for every state of a system at a random time.）

对于上述例子，你可以说整个系统，62.5%的概率系统将处于牛市，31.25%的概率将处于熊市，而6.25%的概率将处于停滞态。

直观上说，你可以想想成在一个状态链上进行随机游走。你在一个状态出，然后你决定你前往下一个状态的方法是在给定当前状态下观察下一个状态的概率分布。

3.3.9 小结

马尔科夫链是一个数学对象，包含一系列状态以及状态之间的转移概率，如果每个状态转移到其他状态的概率只与当前状态相关，那么这个状态链就称为马尔科夫链。有这样一个马尔科夫链之后，我们可以任取一个初始点，然后根据状态转移概率进行随机游走。假设能够找到这样一个马尔科夫链，其状态转移概率正比于我们想要采样的分布（如贝叶斯分析中的后验分布），采样过程就变成了简单地在该状态链上移动的过程。

3.4 连续时间马尔科夫链CTMC

3.4.1 定义

3.4.1.1 定义1:一般定义

和起始时间t无关的话，我们称这是时间齐次的马尔科夫链。这个转移矩阵和离散时间不同的是，离散时间给出的是一步转移概率，但是连续马尔科夫链的转移概率给出的是和时间相关的。

3.4.1.2 定义2

我们考虑连续时间马尔科夫链从一个状态 i开始，到状态发生变化，比如变成j所经过的时间，由于马尔科夫链的马尔科夫性，这个时间是具有无记忆性的，所以这个时间是服从指数分布的。这和离散时间马尔科夫链是密切相关的，离散时间马尔科夫链中的时间是离散时间，因为由无记忆性，所以是服从几何分布的。

这样我们就可以这样定义连续时间马尔科夫链。马尔科夫链是这样的一个过程。

• (i)在转移到不同的状态 $i$前，它处于这个状态的时间是速率为$v_i$的指数分布。
• (ii)当离开状态 $i$时，以某种概率$P_{ij}$进入下一个状态$j$,当然$P_{ij}$满足
  

对比于半马尔科夫链，我们可以发现，连续时间马尔科夫链是一种特殊的半马尔科夫链，在一个状态所待的时间是只不过是一个具体的分布–指数分布，而半马尔科夫链只是说所待的时间是任意的一个随机时间。

3.4.1.3 定义3

3.4.1.4 举例

3.4.2 两个定义(1和2)之间的关系

3.4.3 CTMC的极限概率

3.4.4 CTMC极限分布的应用

3.4.5 时间可逆性

3.4.6 截止的马尔科夫链

四、martingale与markov chain的关系

二者都是随机过程中的一种过程。

Martingale的词本意是指马车上马夫控制马前进的缰绳（如果我记得没错的话），所以从词源来看刻画了一种过程前进（未来）与现在出发点关系的含义。具体来说缰绳的作用是使得马车的前进方向与现在所冲的方向一致，所以在概率上来解释就是未来发生的所有路径可能性的期望值与现在的出发点一致。

从这个意义上来说Matingale的核心是说明了一个过程在向前演化的过程中的稳定性性质。但它并没有说明这个过程是如何到达这一时间点的（是否由上一个时间点所在的位置决定，matingale并没有说明）。再用马车的例子来说，Martingale告诉了你马车在未来是怎么向前走的，中间会有左右的波动（比如马、车夫走神了，路上有坑马要绕开，etc.），但整体来说马是沿着一条直线向前走的。

而马尔科夫过程的核心在于点明了过程的演化是无记忆性的。还拿马车来举栗子，假设车夫喝醉了，他没有意识并在一个很大的广场，马车下一刻前进的方向并不需要是一条直线（经过车夫与马的直线，这种情况下缰绳是绷直的，是martingale），或者说缰绳由于车夫没绷紧是松垮的，这种情况下马车在下一刻可以去任何一个方向，整体上来说前进方向也不必须有什么稳定性规律可循，但整个过程唯一的共性是马迈出前腿的时候，能够到达的所有可能范围，是由它的后腿（你现在所在的位置）决定的（但马可能扭一扭屁股，身子弯曲一下，所以不必须走直线，不必需走直线，不必需走直线！），而并没有由上一步马所在的位置决定，这也就是所谓的无记忆性。

所以从这两个角度来理解，两个名词

有共性：

都从一个过程的全生命角度描画了一个过程的演进性质，

有重叠：
当还是马车例子的时候,martingale也是一个markov（因为虽然走直线，但下一刻的位置还是只由现在决定，只是马身子不能扭曲，不能改变方向），这个例子在概率上最熟悉的模型就是brownian motion了；而反过来，马车未来位置由现在决定，又走直线，所以此时markov process 也是一个martingale （例子还是brownian motion）；

但更重要的是两个过程本质上不是在讲一回事：比如还是马车车夫，喝醉了但走在一个三维空间，这时候它是一个markov process,但是由于方向不确定，此时已经不是martingale而变成了一个local martingale； 而反过来，假设有一个错帧宇宙，空间共享但时间差一天，这时候同一个马车走在不同的宇宙里（但行走轨迹独立），缰绳拉直，此时两架马车都走之前，两架马车组成的系统是一个martingale，但是由于下一时刻前进的方向与宇宙1中的此时有关，也与宇宙2中的昨天有关，所以两架马车组成的系统就不再是一个markov了。

总结一下，brownian motion (wiener process)既是markov process 又是 martingale; 而markov process 与martingale是相交而非包含与反包含关系。只能说你中有我我中有你，但你不属于我我也不属于你

五、Monte Carlo

5.1 蒙特卡罗方法引入

蒙特卡洛（Monte Carlo）方法来自于摩纳哥的蒙特卡洛赌场，许多纸牌类游戏需要计算其胜利的概率。用它作为名字大概是因为蒙特卡罗方法是一种随机模拟的方法，这很像赌博场里面的扔骰子的过程。最早的蒙特卡罗方法都是为了求解一些不太好求解的求和或者积分问题。比如积分：
如果我们很难求解出f(x)的原函数，那么这个积分比较难求解。当然我们可以通过蒙特卡罗方法来模拟求解近似值。(我们可以将蒙特卡洛理解为简单的模拟，通过模拟的情景来计算其发生的概率。)
如何模拟呢？假设我们函数图像如下图:

5.2 概率分布采样

上一节我们讲到蒙特卡罗方法的关键是得到x的概率分布。如果求出了x的概率分布，我们可以基于概率分布去采样基于这个概率分布的n个x的样本集，带入蒙特卡罗求和的式子即可求解。但是还有一个关键的问题需要解决，即如何基于概率分布去采样基于这个概率分布的n个x的样本集。

5.3 接受-拒绝采样

5.4 蒙特卡罗方法小结

使用接受-拒绝采样，我们可以解决一些概率分布不是常见的分布的时候，得到其采样集并用蒙特卡罗方法求和的目的。但是接受-拒绝采样也只能部分满足我们的需求，在很多时候我们还是很难得到我们的概率分布的样本集。比如：

从上面可以看出，要想将蒙特卡罗方法作为一个通用的采样模拟求和的方法，必须解决如何方便得到各种复杂概率分布的对应的采样样本集的问题。而马尔科夫链就是帮助找到这些复杂概率分布的对应的采样样本集的白衣骑士。

5.5 应用场景

通常情况下，许多概率的计算非常复杂甚至是无法计算的，但是我们可以通过计算机对期待发生的场景进行大量的模拟，从而计算出其发生的概率。简单的公式表达为：

学习过概率的同学应当知道这就是典型的频率学派的观点。

对于蒙特卡洛方法，经典的例子就是计算$\pi$值。在（-1，1）之间随机取两个数，如果在单位圆内，则记一次，在圆外则不计入次数。

import random
import numpy as np
from math import sqrt, pi

num1 = 1000 # simulation times 
freq1 = 0
x1 = []
y1 = []
for i in range(1,num1+1):
        x_i, y_i = random.uniform(-1, 1), random.uniform(-1, 1)
        square_distance = x_i**2+y_i**2
        x1.append(x_i)
        y1.append(y_i)
        if square_distance <= 1: # 圆的方程x^2+y^2=1
                freq1 += 1

simulation_pi1 = 4*freq1/num1;

num2 = 10000 # simulation times 
freq2 = 0
x2 = []
y2 = []
for i in range(num2):
        x_i, y_i = random.uniform(-1, 1), random.uniform(-1, 1)
        square_distance = x_i**2+y_i**2
        x2.append(x_i)
        y2.append(y_i)
        if square_distance <= 1: # 圆的方程x^2+y^2=1
                freq2 += 1

simulation_pi2 = 4*freq2/num2;

% matplotlib inline
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(8,8))

theta = np.linspace(0,2*pi,500)
x,y = np.cos(theta), np.sin(theta)
ax.plot(x, y, color='red', linewidth=1.0)

ax.scatter(x1,y1,color='blue',alpha=0.75)
ax.set_xlim(-1,1);
ax.set_ylim(-1,1);

结果为：

真实的pi值： 3.141592653589793
蒙特卡洛方法下得出的pi值（模拟1000次）： 3.188
相对误差（模拟1000次）： 1.4771917153924754 %
蒙特卡洛方法下得出的pi值（模拟10000次）： 3.1292
相对误差（模拟10000次）： 0.39447041536821975 %

2. 小结

蒙特卡洛就是一种模拟事情发生的方法。

六、MCMC

蒙特卡洛马尔可夫链 Markov Chain Monte Carlo简称MCMC，是一个抽样方法，用于解决难以直接抽样的分布的随机抽样模拟问题。用来在概率空间，通过随机采样估算兴趣参数的后验分布。

简而言之，就是用马尔科夫链方法来抽样。

抽样算法的主要任务是找到符合给定分布p(x)的一系列样本。对于简单的分布，可以通过基本的抽样算法进行抽样。大多数分布都是不容易直接抽样的，马尔可夫链蒙特卡罗算法解决了不能通过简单抽样算法进行抽样的问题，是一种重要的实用性很强的抽样算法。

马尔可夫链蒙特卡罗算法（简写为MCMC）的核心思想是找到某个状态空间的马尔可夫链，使得该马尔可夫链的稳定分布就是我们的目标分布p(x)。这样我们在该状态空间进行随机游走的时候，每个状态x的停留时间正比于目标概率p(x)。在用MCMC进行抽样的时候，我们首先引进一个容易抽样的参考分布q(x)，在每步抽样的过程中从q(x)里面得到一个候选样本y, 然后按照一定的原则决定是否接受该样本，该原则的确定就是要保证我们得到的原本恰好服从p(x)分布。

在了解什么是MCMC的时候，我们需要考虑一个情况。举例，我们已经知道一个分布（例如beta分布）的概率密度函数PDF，那么我们怎么样从这个分布中提取样本呢？

MCMC给我们提供了一种方式来从概率分布中进行采样，而这种方法在贝叶斯统计中的后验概率进行采样十分方便。

1. 动因

MCMC通常用于解决高维度的积分和最优化问题，这两种问题也是机器学习、物理、统计、计量等领域的基础。比如：

2. MCMC的定义

Wikipedia的解释如下:

Markov Chain Monte Carlo (MCMC) methods are a class of algorithms for sampling from a probability distribution based on constructing a Markov chain that has the desired distribution as its stationary distribution. The state of the chain after a number of steps is then used as a sample of the desired distribution. The quality of the sample improves as a function of the number of steps.
马尔可夫链蒙特卡罗方法是一类以期望分布为平稳分布的马尔可夫链为基础，对概率分布进行抽样的算法。经过一系列步骤之后，链的状态将用作期望分布的样本。样品的质量随着步骤数量的增加而提高。

MCMC方法提供了可以创建一个以Beta分布为其平稳分布（stationary distribution）的马尔科夫链的算法，从而使取样变得简单，因为我们可以从一个均匀分布中取样(这相对容易)。

如果我们基于一些算法，重复地从一个随机的状态开始，遍历到下一个状态，我们将会创建一个马尔可夫链，一个以beta分布作为其平稳分布的马尔可夫链。

这类算法，一个典型的代表就是Metropolis-Hastings Algorithm算法。

MCMC由梅特罗波利斯(Metropolis)于1949年基于马尔可夫链的基本性质提出，下面介绍一下与马尔可夫链相关的性质。

3. 相关性质

一、稳定分布

稳定分布是指当我们给定状态空间的一个初始分布${\pi }_0$以后，按照转移矩阵进行跳转最终达到的稳定状态。

即每个状态的流出概率等于该状态注入的概率，该方程称为全局平衡方程。满足该方程的分布称为转移矩阵A的稳定分布。

二、细致平衡

对于一个稳定分布，我们可以给其一个更强的条件限制，使得任意一个状态$i$满足如下条件，

下面我们介绍基本的梅特罗波利斯算法(Metropolis algorithm)。假设p(x)是我们的目标分布函数，我们想得到一系列服从该分布的样本。我们考虑样本的产生过程构成一个马尔可夫链，并且让p(x)是该马尔可夫链的稳定分布，那么该样本序列就服从p(x)。现在p(x)是已知的，问题的关键在于如何构造这个产生过程，也即如何构造马尔可夫链。


即目标概率p(x)满足细致平衡方程，因此p(x)是转移概率的稳定分布。

最后，回顾一下梅特罗波利斯抽样的主要思想。我们首先构造了一个马尔可夫链，接着证明了p(x)满足其细致平衡方程，进而说明p(x)是该链的稳定分布。然后将抽样的过程看成是在马尔可夫链状态空间进行跳转的过程，跳转的候选状态由参考分布q(x)产生。最后 得到一个的跳转序列，该序列在每个状态x的停留时间与p(x)成比，即服从p(x)分布。

4. Metropolis-Hastings (MH) Algorithm算法

梅特罗波利斯－黑斯廷斯算法（英语：Metropolis–Hastings algorithm）是统计学与统计物理中的一种马尔科夫蒙特卡洛（MCMC）方法，用于在难以直接采样时从某一概率分布中抽取随机样本序列。得到的序列可用于估计该概率分布或计算积分（如期望值）等。梅特罗波利斯－黑斯廷斯或其他MCMC算法一般用于从多变量（尤其是高维）分布中采样。对于单变量分布而言，常会使用自适应判别采样（adaptive rejection sampling）等其他能抽取独立样本的方法，而不会出现MCMC中样本自相关的问题。

该算法的名称源于美国物理学家尼古拉斯·梅特罗波利斯与加拿大统计学家W·K·黑斯廷斯

为了更形象地理解这个算法，我们用下面这个例子来类比。假设我们想知道某个湖的水容量以及这个湖中最深的点，湖水很浑浊以至于没法通过肉眼来估计深度，而且这个湖相当大，网格近似法显然不是个好办法。为了找到一个采样策略，我们请来了两个好朋友小马和小萌。经过讨论之后想出了如下办法，我们需要一个船（当然，也可以是竹筏）和一个很长的棍子，这比声呐可便宜多了，而且我们已经把有限的钱都花在了船上。

（1）随机选一个点，然后将船开过去。

（2）用棍子测量湖的深度。

（3）将船移到另一个地点并重新测量。

（4）按如下方式比较两点的测量结果。

• 如果新的地点比旧的地点水位深，那么在笔记本上记录下新的测量值并重复过程（2）。
• 如果新的地点比旧的地点水位浅，那么我们有两个选择：接受或者拒绝。接受意味着记录下新的测量值并重复过程（2）；拒绝意味着重新回到上一个点，再次记录下上一个点的测量值。

如何决定接受还是拒绝新的测量值呢？这里的一个技巧便是使用Metropolis-Hastings准则，即接受新的测量值的概率正比于新旧两点的测量值之比。

事实上，理论保证了在这种情形下，如果我们能采样无数次，最终能得到完整的后验。幸运地是，实际上对于很多问题而言，我们只需要相对较少地采样就可以得到一个相当准确的近似。

5. 数学表达

6. 算法介绍

7. 代码实现MCMC采样beta分布

import random
# Lets define our Beta Function to generate s for any particular state. 
# We don't care for the normalizing constant here.

def beta_s(w,a,b): # beta distribution pdf
    return w**(a-1)*(1-w)**(b-1)

# This Function returns True if the coin with probability P of heads comes heads when flipped.
def random_coin(p):
    unif = random.uniform(0,1)
    if unif>=p:
        return False
    else:
        return True

# This Function runs the MCMC chain for Beta Distribution.
def beta_mcmc(N_hops,a,b):
    states = []
    cur = random.uniform(0,1)
    for i in range(0,N_hops):
        states.append(cur)
        next = random.uniform(0,1)
        ap = min(beta_s(next,a,b)/beta_s(cur,a,b),1) # Calculate the acceptance probability
        if random_coin(ap):
            cur = next
    return states[-1000:] # Returns the last 1000 states of the chain

在定义了相关函数后，我们来检查一下我们采样的结果和实际beta概率分布之间的差距

import numpy as np
from matplotlib import pylab as pl
import scipy.special as ss
% matplotlib inline
plt.rcParams['figure.figsize'] = (17.0, 4.0)

# Actual Beta PDF.
def beta(a, b, i):
    e1 = ss.gamma(a + b)
    e2 = ss.gamma(a)
    e3 = ss.gamma(b)
    e4 = i ** (a - 1)
    e5 = (1 - i) ** (b - 1)
    return (e1/(e2*e3)) * e4 * e5

# Create a function to plot Actual Beta PDF with the Beta Sampled from MCMC Chain.
def plot_beta(a, b):
    Ly = []
    Lx = []
    i_list = np.mgrid[0:1:100j]
    for i in i_list:
        Lx.append(i)
        Ly.append(beta(a, b, i))
    pl.plot(Lx, Ly, label="Real Distribution: a="+str(a)+", b="+str(b))
    plt.hist(beta_mcmc(100000,a,b),density=True,bins =25, 
            histtype='step',label="Simulated_MCMC: a="+str(a)+", b="+str(b))
    pl.legend()
    pl.show()
    
plot_beta(0.1, 0.1)
plot_beta(1, 1)
plot_beta(2, 3)

相关图形如下所示：

从上面的采样值可以看出，我们的采样值还是很接近beta分布本身。虽然我们用了beta分布最为例子，但实际上其他的分布也是可以类似来进行采样。

针对那些我们无法直接利用共轭分布得到的后验分布，特别是高维的分布，MCMC方法更显的尤为重要。

8. MH Algorithm的进一步说明

从贝叶斯的角度看，Metropolis-Hastings算法使得我们能够从任意分布中以概率p(x)得到采样值，只要我们能算出某个与p(x)成比例的值。这一点很有用，因为在类似贝叶斯统计的许多问题中，最难的部分是计算归一化因子，也就是贝叶斯定理中的分母。Metropolis-Hastings算法的步骤如下：

最后，我们会得到一连串记录值，有时候也称采样链或者迹。如果一切都正常进行，那么这些采样值就是后验的近似。在采样链中出现次数最多的值就是对应后验中最可能的值。该过程的一个优点是：后验分析很简单，我们把对后验求积分的过程转化成了对采样链所构成的向量求和的过程。

强烈建议阅读以下博文，作者用一个简单的例子实现了metropolis方法，并将其用于求解后验分布，文中用非常好看的图展示了采样的过程，同时简单讨论了最初选取的步长是如何影响结果的。
https://github.com/findmyway/Bayesian-Analysis-with-Python/blob/master/MCMC-sampling-for-dummies.ipynb

9. 一句话总结MCMC

MCMC generates samples from the posterior distribution by constructing a reversible Markov-chain that has as its equilibrium distribution the target posterior distribution.
https://twiecki.io/blog/2015/11/10/mcmc-sampling/

10. MCMC如何解决具体问题

• 寻找p(x)最大值（simulated annealing）
• 转换核的混搭和循环
• Gibbs 抽样
• Monte Carlo EM（期望-最大化）
• 辅助变量抽样（Auxiliary variable samplers)
• 可逆跳跃MCMC（Reversible jump MCMC）

详见

———————————————–
链接：
https://zhuanlan.zhihu.com/p/83314877
https://www.zhihu.com/question/26887292/answer/310904395
https://zh.wikipedia.org/wiki/%E9%9E%85_(%E6%A6%82%E7%8E%87%E8%AE%BA)
https://www.zhihu.com/question/52778636/answer/133133860
https://www.zhihu.com/question/31173033/answer/113202955
https://zhuanlan.zhihu.com/p/116725922
https://zhuanlan.zhihu.com/p/86995916
https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm

今天的文章martingale、markov chain、Monte Carlo、MCMC分享到此就结束了，感谢您的阅读。