我们常常用期望值来估计一些未知的参数，然而期望值的积分有时候很难算，此时就出现了蒙特卡洛方法求取近似数值解。但是蒙特卡洛遇到了采样难的问题，因此出现了各种采样的方法，MCMC就是其中之一。
对于一个随机变量
X
X
X而言，它可以取许多许多数值，但是取到每个数值是有一定概率的。取到某些值的概率比较大，取到某些值的概率比较小，取到每个值的概率画成一张图就是概率密度图。通常，我们也会用随机变量的期望
E
(
x
)
E(x)
E(x)来描述这个随机变量的性质，期望
E
(
x
)
E(x)
E(x)的计算公式如下：
E
(
X
)
=
∫
x
p
(
x
)
d
x
E(X) = \int xp(x)dx
E(X)=∫xp(x)dx
如果这时候存在一个关于随机变量的函数
g
(
x
)
g(x)
g(x)，实际上
Y
=
g
(
x
)
Y=g(x)
Y=g(x)也是个随机变量。此时期望
E
(
Y
)
=
E
[
g
(
x
)
]
E(Y)=E[g(x)]
E(Y)=E[g(x)]可以用如下公式计算：
E
(
Y
)
=
E
[
g
(
x
)
]
=
∫
g
(
x
)
p
(
x
)
d
x
E(Y) = E[g(x)] = \int g(x)p(x)dx
E(Y)=E[g(x)]=∫g(x)p(x)dx
其中，
p
(
x
)
p(x)
p(x)是随机变量
X
X
X的概率密度函数。
在这里，期望是一个积分的形式，当这个积分难以写出显式解（表达式）的时候，就可以用蒙特卡洛模拟的方法去近似计算这个积分值，从而计算出期望。
蒙特卡洛方法给出的计算积分的方法是：生成一系列服从
p
(
x
)
p(x)
p(x)分布的
x
i
x_i
xi，然后算出一系列
g
(
x
i
)
g(x_i)
g(xi)，然后用这些值的平均值近似代替我想要的积分值。这篇博客里面写了证明，用公式描述就是：
E
(
Y
)
≈
1
N
Σ
i
=
1
N
g
(
x
i
)
w
h
e
r
e
x
i
∼
p
(
x
)
E(Y) \approx \frac{1}{N}\Sigma_{i=1}^{N}g(x_i) \\ where \ x_i \sim p(x)
E(Y)≈N1Σi=1Ng(xi)where xi∼p(x)
那么问题来了，应该怎么取服从
p
(
x
)
p(x)
p(x)分布的
x
i
x_i
xi呢？这里的
p
(
x
)
p(x)
p(x)是概率密度函数（Probability Density Function, PDF），实际上如果能把概率密度函数写成概率分布函数（Cumulative Distribution Function, CDF）就方便了，CDF实际上就是PDF的积分。举个例子：
假设我们要取服从正态分布 N ( 0 , 1 ) \mathcal{N}(0,1) N(0,1)的100个样本 [ x 1 , . . . , x 100 ] [x_1,...,x_{100}] [x1,...,x100]。已知其PDF，由PDF计算出CDF如下图：
从CDF纵坐标上我们均匀取100个点（一般是用随机数均匀取值，见下图红色箭头），这100个点对应的横坐标值就是我们想要的 x i x_i xi（见下图蓝色箭头）。
上文中提到的"生成一系列服从 p ( x ) p(x) p(x)分布的 x i x_i xi"就叫做采样，而写不出PDF，写不出CDF等情况，无法生成这样的 x i x_i xi就叫做无法采样。CDF是PDF的积分，如果这个PDF本身就很复杂，或者无法写出表达式，那CDF就写不出来了，就没法用上面的方法取服从服从 p ( x ) p(x) p(x)分布的 x i x_i xi了。举个无法写出CDF的例子：
下面这里例子里面，需要注意的是，上文中的 x i x_i xi变成了高维的随机变量。也就是说 x i = [ a 1 i , a 2 i , . . . , a 100 i ] x_i=[a_{1i},a_{2i},...,a_{100i}] xi=[a1i,a2i,...,a100i]。
已知一个联合概率分布 p ( a 1 , a 2 , . . . , a 100 ) p(a_1,a_2,...,a_{100}) p(a1,a2,...,a100)，其中每个 a a a的取值范围是 1 , 2 , 3 , . . . 50 1,2,3,...50 1,2,3,...50，其计算公式如下：
KaTeX parse error: No such environment: align at position 8: \begin{̲a̲l̲i̲g̲n̲}̲ & p(a_1,a_2,..…这个联合概率分布，分子上的 Z Z Z是一个归一化的常数，但是想算出这个常数就要把所有 a a a的取值都算一遍。分母上的 E E E写的比较简单，想表达的意思就是，这个函数我们知道，并且可以算出来，实际上可以用任何一个简单或者复杂的函数替代。
我们来分析一下这个联合概率分布：如果我们知道了一组 [ a 1 , a 2 , . . . , a 100 ] [a_1,a_2,...,a_{100}] [a1,a2,...,a100]，那么直接带入公式。分子可以计算出来，但是分母要遍历所有 a a a的取值。实际上，分母硬算是可以算出来的，但是太复杂了。况且，把维度加一加，把取值范围加一加，把 E E E整一整，还能更加复杂，一年两年不知道能不能搞出来。这个例子里面，PDF都算不出，更别说CDF了，就更别提采样了。
为了解决无法采样的情况，就有了接受-拒绝采样、重要性采样、MCMC等方法。下面就写写这些方法是怎么解决这个问题的。
接受-拒绝采样并不能解决上面那个例子中的问题，但是它能解决“可以写出PDF，但是写不出CDF”的问题。也就是说，我们已知一个 x i x_i xi，就可以算出对应的 p ( x i ) p(x_i) p(xi)的数值，但是 ∫ p ( x ) d x \int p(x)dx ∫p(x)dx算不出来。
回到我们最开始求期望的那个问题：
E
(
Y
)
=
E
[
g
(
x
)
]
=
∫
g
(
x
)
p
(
x
)
d
x
E(Y) = E[g(x)] = \int g(x)p(x)dx
E(Y)=E[g(x)]=∫g(x)p(x)dx
既然我们无法"生成一系列服从
p
(
x
)
p(x)
p(x)分布的
x
i
x_i
xi"，那我们"生成一系列服从
q
(
x
)
q(x)
q(x)分布的
x
i
x_i
xi"，然后从这些
x
i
x_i
xi中挑一部分，使得挑出来的
x
i
x_i
xi服从
p
(
x
)
p(x)
p(x)分布。其中
q
(
x
)
q(x)
q(x)是一个我们已知且可采样的分布（例如
N
(
0
,
1
)
\mathcal{N}(0,1)
N(0,1)）。具体的操作如下：
以某某概率接受其实很简单。假设以0.8的概率接受，那么就生成一个随机数（0-1之间均匀分布），如果这个数大于0.8，我就不接受，小于0.8就接受。
因为生成的 x i x_i xi有一部分被被扔掉了（被拒绝了），所以这个方法叫做接受-拒绝采样。
重要性采样并不能解决上面那个例子中的问题，但是它能解决“可以写出PDF，但是写不出CDF”的问题。也就是说，我们已知一个 x i x_i xi，就可以算出对应的 p ( x i ) p(x_i) p(xi)的数值，但是 ∫ p ( x ) d x \int p(x)dx ∫p(x)dx算不出来。
回到我们最开始求期望的那个问题，我们把求期望的公式换个写法：
KaTeX parse error: No such environment: align at position 8: \begin{̲a̲l̲i̲g̲n̲}̲ E(Y) = E[g(x)]…
这里的
s
(
x
)
s(x)
s(x)是一个我们已知的概率密度函数PDF（例如
N
(
0
,
1
)
\mathcal{N}(0,1)
N(0,1)），而给定一个
x
x
x，
h
(
x
)
h(x)
h(x)可以被算出来。因此只需要生成一系列“服从
h
(
x
)
h(x)
h(x)分布的
x
i
x_i
xi”，然后根据：
E
(
Y
)
≈
1
N
Σ
i
=
1
N
h
(
x
i
)
w
h
e
r
e
x
i
∼
s
(
x
)
E(Y) \approx \frac{1}{N}\Sigma_{i=1}^{N}h(x_i) \\ where \ x_i \sim s(x)
E(Y)≈N1Σi=1Nh(xi)where xi∼s(x)
即可算出期望值。
因为 h ( x ) = g ( x ) p ( x ) s ( x ) h(x) = g(x) \frac{p(x)}{s(x)} h(x)=g(x)s(x)p(x)中存在 p ( x ) s ( x ) \frac{p(x)}{s(x)} s(x)p(x)，就像是给原来的 g ( x ) g(x) g(x)加了个权重，所以叫做重要性采样。
MCMC方法可以解决上面那个例子中的问题，虽然分母算不出来，但是分子能算啊。MCMC中使用 p ( x ′ ) p ( x ) \frac{p(x')}{p(x)} p(x)p(x′)带入计算，恰好使得分母抵消了。
MCMC方法很鸡贼，它构造了一个马尔科夫链。这个马尔科夫链在进入稳态后，虽然在各个状态之间跳转，但是跳转有一定规律：落在各个状态的概率服从某一分布。这个分布就是我们希望采样的 p ( x ) p(x) p(x)。而我们只需要让这个马尔科夫链不断跳转，然后把每一时刻它所在的状态记录下来，记录下来的状态 x i x_i xi就是我们想要的服从 p ( x ) p(x) p(x)分布的 x i x_i xi。
MCMC主要包括Metropolis方法、Metropolis-Hasting方法以及Gibbs Sampling方法。其中M-H方法是对M方法的改进，而Gibbs是在做图像去噪的时候恰好采用了类似的思路，顺带改进了一下。下面就以H-C算法讲解一下其实现思路，算法原理可以在各种各样的博客里面找到。
首先准备一个提议分布 q ( x m ∣ x n ) q(x_{m}|x_{n}) q(xm∣xn)，这个分布可以随便取，它代表的是从状态 x n x_{n} xn转移到状态 x m x_m xm的概率。换句话说就是，我们假设从状态 x n x_{n} xn转移到状态 x m x_m xm的概率是 q ( x m ∣ x n ) q(x_{m}|x_{n}) q(xm∣xn)。
这个假设的分布
q
q
q显然不是我们想要的分布
p
p
p，因此这里借鉴一下接收-拒采样的想法。接收-拒绝抛弃一部分
x
i
x_i
xi，而H-C不接受一部分状态转移。换句话说就是，按照
α
\alpha
α的概率接受状态转移。其中：
α
=
m
i
n
(
1
,
p
(
x
′
)
q
(
x
∣
x
′
)
p
(
x
)
q
(
x
′
∣
x
)
)
\alpha = min(1, \frac{p(x')q(x|x')}{p(x)q(x'|x)})
α=min(1,p(x)q(x′∣x)p(x′)q(x∣x′))
[1] 如何理解 重要性采样(importance sampling)
[2] particle filtering—粒子滤波（讲的很通俗易懂）
[3] N. Metropolis et al (1953). Equations of state calculations by fast computing machines（Metropolis算法原始论文）
[4] Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications（Metropolis-Hasting算法原始论文）
[5] Geman, S. and Geman, D. (1984). Stochastic relaxation, Gibbs distributions and the Bayesian restoration of images（Gibbs Sampling算法原始论文）
