Markov Chain Monte Carlo（MCMC采样）

近似推断（确定性的变分VI，不确定的采样mcmc），而MCMC由两个MC组成，即蒙特卡罗方法（Monte Carlo Simulation，简称MC）和马尔科夫链（Markov Chain ，也简称MC）。

Monte Carlo方法：
为什么要使用它？
$$\theta ={\int }_a^{b}f(x)dx$$诸如上式求积分，但如果f（x）相当复杂时积分将会很难求解，但是我们知道定积分表示图形的面积，而且根据大数定理采样的最终分布往往会收敛为原分布，所以对于 $\theta$ 的估计往往可以采用 随机采样 的方法来模拟求解近似值。（Monte Carlo是某赌场的名字，实质上它表示的就是随机模拟方法）。比如最简单的扔硬币，在重复10000次甚至更多的时候，它最终会告诉我们正反的概率几乎50%。

Monte Carlo方法主要有以下：
1.概率分布采样
如果有了函数的pdf（probability density function，概率密度函数），而且可以求得该函数的cdf（cumulative distribution function，累积分布函数），那么可以直接在[a,b]上采样，并且可以得到相应的值。但是如果pdf非常复杂，cdf是求不出来的，所以这种方法有很大的局限性。

2.adjection sampling（接受–拒绝采样）
对于类似上图难以处理的p（x），接受–拒绝采样使用相对来说比较简单的分布来做概率分布采用，即构造一个q（x）的k倍使其 p(x)总在 kq(x)的下方（因为q和p的概率分布都为1，所以q需要乘以k来保证一定可以包含p）。然后对于随机采样 $z_0,z_1,...z_{n-1}$，如果样本落在了p（x）中就接受该采样，如果落入阴影部分则拒绝该采样。（某点的$\frac{p(x)}{kq(x)}$的值可以定义为接受率，低于该值就接受采样样本）。
但是这种方式非常依靠假设出的形态q，而假设出非常接近的函数q和合适的k都太难，实际中（特别是高维情况下）效率太低，在极端情况下可能大多数的样本点都无法被接受。

3.importance sampling（重要性采样）
在q和k都太难得到的情况下再怎么办？回顾采样的目的求得参数，诸如扔硬币，实际上我们想要求得它出现的期望E，那么就通过期望。$E[f(x)]=\int xf(x)p(x)dx$,其中x~p，从这个式子可以启发性的得到，如果q太难生成，那么就引入另一个函数，即$$g(x)=\int xf(x)p(x)/q(x)\ast q(x)dx=\int xg(x)q(x)dx=f(x)p(x)q(x)=f(x)w(x)$$其中w(x)=p(x)q(x) 就叫做Importance Weight。关键就在这里，通过转换将不好求的期望，变成了一个在另一个分布下相对好求的期望。通过他们之间的商即接收率来进行加权采样，即样本采样时按重要性加权了，就可以一定程度上避免接受–拒绝采样的那种极端情况。

4.sampling-importance-sampling
先采样，按权重分布（看成概率值）再采样，可以使最后的结果本身的权重一样没有差别，而它们之间的采样个数不一样了。

Markov Chain：
但是通用的Monte Carlo方法想要得到合适的样本集进行随机模拟求解还是很困难。

承认客观世界中有这样一种现象，其未来由现在决定的程度，使得我们关于过去的知识丝毫不影响这种决定性。这种在已知“现在”的条件下，“未来”与“过去”彼此独立的特性就被称为马尔科夫性，具有这种性质的随机过程就叫做马尔科夫过程。–辛钦

即马尔科夫链的定义是某一时刻状态转移的概率只依赖于它的前一个状态。即假定今天的心情只与昨天的心情有关，而和前天，大前天都没有关系。
即对于序列链$...X_{t-2},X_{t-1},X_t,X_{t+1},...$…t+1时刻的概率：
$$P(X_{t+1}|...X_{t-2},X_{t-1},X_t)=P(X_{t+1}|X_t)$$只与 t 时刻的状态有关。于是可以得到某一时刻的状态转移矩阵$$P=\left(\begin{array}{ccc}0& 1& 0.6\\ 0& 0.1& 0.9\\ 0.6& 0.4& 0\end{array}\right)$$故它的采样方法就是在服从f（x）分布的条件下在样本之间“跳来跳去”（由前一时刻所决定）。再次回顾使用MCMC的目的，重要在于求期望难求，才要通过采样方法，而对于复杂的函数又无法直接采样，所以需要间接的用MCMC采样。那么如果得到了某个平稳分布所对应的马尔科夫链状态转移矩阵，就很容易求出样本集了。（状态转移矩阵收敛到的稳定概率分布将会与初始状态概率分布无关，就可以得到符合对应稳定概率分布的样本了。）

故问题变为如何求符合条件的转移矩阵，条件就在于需要满足马尔科夫链的细致平稳条件（detailed balance）。
设任意时刻的概率分布为$\pi （x）$，经过n次采样后，马尔可夫链收敛平稳状态，那么有：
$${\pi }_n(x)={\pi }_{n+1}(x)={\pi }_{n+2}(x)=...=\pi (x)$$
而且每个分布：
$${\pi }_i(x)={\pi }_{i-1}(x)P={\pi }_{i-2}(x)P^2={\pi }_0(x)P^i$$
P为状态转移矩阵，而细致平稳条件的定义就是：$$\pi (i)P(i,j)=\pi (j)P(j,i)$$即 i 跳到 j 和 j 跳回到 i 的值相等。
@@@跳来与跳去的概率是一样的就是细致平稳分布，所以满足detailed balance的一定是平稳分布，但是反过来不一定。
那么如何寻找转移矩阵（设为Q）满足细致平稳条件呢？$P(i,j)=Q(i,j)\alpha (i,j)$，即先引入一个$\alpha (i,j)$
$$\pi (i)Q(i,j)\alpha (i,j)=\pi (j)Q(j,i)\alpha (j,i)$$
且：$$\alpha (i,j)=\pi (j)Q(j,i)$$

$$\alpha (j,i)=\pi (i)Q(i,j)$$$\alpha (i,j)$为为接受率。取值在[0,1]之间，可以理解为一个概率值。即目标矩阵P可以通过任意一个马尔科夫链状态转移矩阵Q以一定的接受率获得。（类似于接受–拒绝）所以类似的如果α过小，将会导致频繁拒绝转移，马氏链难以收敛，所以考虑放大α，使得α(i,j)和α(j,i)中大的值为1，小的值放大为$\frac{\pi (j)P(j,i)}{\pi (i)P(i,j)}$，那么 α ( i , j ) = m i n { 1 , π ( j ) P ( j , i ) π ( i ) P ( i , j ) } \alpha(i,j) = min\{1,\frac{\pi(j)P(j,i)}{\pi(i)P(i,j)}\} α(i,j)=min{1,π(i)P(i,j)π(j)P(j,i)​}。

M-H采样
使用 α ( i , j ) = m i n { 1 , π ( j ) P ( j , i ) π ( i ) P ( i , j ) } \alpha(i,j) = min\{1,\frac{\pi(j)P(j,i)}{\pi(i)P(i,j)}\} α(i,j)=min{1,π(i)P(i,j)π(j)P(j,i)​}就是M-H采样了，它基本解决了手动取参所带来的低效率性。但在数据特征非常的多，M-H采样的接受率计算太复杂，算法效率将会变低。

Gibbs采样
先看它的采样方式：$(x_1^{(1)},x_2^{(1)})\to (x_1^{(1)},x_2^{(2)})\to (x_1^{(2)},x_2^{(2)})\to ...\to (x_1^{(n_1+n_2-1)},x_2^{(n_1+n_2-1)})$每次仅仅变动一维，而且使用当前状态值。这样做的好处是：对于$A(x_1^{(1)},x_2^{(1)})$和$B(x_1^{(1)},x_2^{(2)})$有：

$$\pi (x_1^{(1)},x_2^{(1)})\pi (x_2^{(2)}|x_1^{(1)})=\pi (x_1^{(1)})\pi (x_2^{(1)}|x_1^{(1)})\pi (x_2^{(2)}|x_1^{(1)})$$

$$\pi (x_1^{(1)},x_2^{(2)})\pi (x_2^{(1)}|x_1^{(1)})=\pi (x_1^{(1)})\pi (x_2^{(2)}|x_1^{(1)})\pi (x_2^{(1)}|x_1^{(1)})$$
所以：
$$\pi (x_1^{(1)},x_2^{(1)})\pi (x_2^{(2)}|x_1^{(1)})=\pi (x_1^{(1)},x_2^{(2)})\pi (x_2^{(1)}|x_1^{(1)})$$
那么：
$$\pi (A)\pi (x_2^{(2)}|x_1^{(1)})=\pi (B)\pi (x_2^{(1)}|x_1^{(1)})$$
天然满足平稳细致条件。换句话说即Gibbs采样是特殊的M-H，且它的接受率恒为100%。所以当前时代基本都是Gibbs的世界。

import pymc3 as pm

X, y = linear_training_data()
with pm.Model() as linear_model:
    weights = pm.Normal('weights', mu=0, sigma=1)
    noise = pm.Gamma('noise', alpha=2, beta=1)
    y_observed = pm.Normal('y_observed',
                mu=X.dot(weights),
                sigma=noise,
                observed=y)

    prior = pm.sample_prior_predictive()
    posterior = pm.sample()
    posterior_pred = pm.sample_posterior_predictive(posterior)

贝叶斯网络与马尔科夫随机场
两者都是在联合概率、条件概率、边际概率上来回驰骋。贝叶斯网络是有向图（利用条件概率分布进行转移），而马尔可夫是无向图（利用转移矩阵进行转移）。

能量函数
对于大变化的取值都建矩阵对建模来说太困难，就出现了能量函数。在物理学中，能量越大的物质存在概率就越小，能量越小存在的概率就越大。在玻尔兹曼机中就用了这样的思想。

条件随机场CRF
之所以叫条件，就是参与计算的两部分，一个是可观察，一个是隐含变量。