声明：此方法建立在著名的马尔可夫链蒙特卡洛采样算法（Markov chain Monte Carlo，MCMC）之上，并一改巴普诺夫把妹法和薛定谔把妹法的送餐设定，而是虚构了一个真实的故事场景，令学习者更加感同身受，可以说是一种更加科学的追女生方法。
正所谓得人心者得天下，送餐什么的太乃衣服了，懂得女生心理才是王道，下面我便将此法传授给大家。
话说….故事是这样展开的…网络上有一个很好的姑娘。正因为是很好的姑娘，所以她的追求者多呀，于是该选择谁作为伴侣便成了问题。终于有一天呢，她心生一计：“我应当选择一个最懂我的人，那什么样的人是最懂我的呢？当然是最了解我心情的吧”。于是呢，就在微博上发了一个英雄帖，内容如下：小女之心情，变换莫测，诡谲云涌，喜怒无常，可谓是风云变化几席上，蛟鼋出波澜杆前。然则每日有定值，为实数，知其概率密度，却不知其故也，今寻有识之士，可日得一心情值之概率密度，为期一载，知我心情之期望者，吾嫁之，白首不悔。
显然，这是一个十分棘手的问题，这是要通过我们自身的采样来估计姑娘心情值的概率分布呀… 这可真是太难了，许多的追求者因此望而却步…但是咱作为有(老)识(色)者(批)怎可轻言放弃，于是打算用科学的方法破解之。
你仔细分析了一番，发现这个问题相当于要知道小姐姐心情分布$\pi(x)$的均值。但是吧，正所谓女人心海底针，瞎猜肯定是不行的了，于是你马上想到了一种通过随机抽样来近似解决计算问题的方法——蒙特卡洛算法。
既然直接算不出来，那我可以根据$\pi(x)$来进行随机采样，得到一个样本集$x_1,x_2,…x_n$，然后计算这个样本集的期望来近似不就可以了吗？ \(\mathbb{E} = \frac{1}{n} \sum_{i=1}^n x_i\)
但是很显然，你想多了，基于$\pi$ 采样，你采样个屁呀，你懂少女心吗？还采样…
于是你又灵机一动，基于马尔科夫链能收敛到平稳分布这个性质，有了一个绝妙的想法：如果我们能够构造一个转移矩阵为P的马尔可夫链，使得该链的平稳分布恰好是$\pi(x)$ ,那么我们就可以从任何一个初始状态$x_0$ 出发，沿着转移矩阵转移，得到一个转移序列$x_0,x_1,x_2,…,x_n,x_{n+1}…$ ,如果在第$n$ 步已经收敛了，不就得到了$\pi(x) $ 的样本$x_n, x_{n+1}，x_{n+2} …$ 了吗？
正当你要着手设计这个状态转移矩阵的时候，不幸的又被现实给打脸了，你根本就不知道小姐姐每天的心情值的变化，而只知道某天某个心情值的概率密度而已… 这根本就不可能得出状态转移矩阵，况且，要是你连转移矩阵都知道了，那你直接可以就能算出稳定的概率分布了…
于是你又绞尽脑汁的想啊想…
就在你一筹莫展的时候，突然灵光乍现，想起一个以前用常见分布采样来估计不常见分布的方法——接受拒绝采样：比如我们难以对分布$p(x)$进行采样，就先基于常见的好采样的分布$q(x)$来进行采样，并选择一个常数$k$，使得对于任意的$x$都有$kq(x) > p(x)$，然后就好办了…
如上图所示，我们要直接基于$p(x)$采样很难，可以先构造这样的一个$q(x)$， 随机采样到一个样本$x_0$， 然后再从均匀分布$(0, kq(x_0))$ 上采样，得到一个值$u$, 如果$u > p(x_0)$ 则拒绝，否则就接受这个样本$x_0$… 重负这个过程，就可以得到一系列基于$p(x)$ 随机抽样的样本，x_0, x_1,x_2,…
你兴致勃勃的采样了一个$x_i$去询问小姐姐相应的概率密度，打算大干一场… 却无奈的发现这个满足条件的$q(x)$和$k$ 实在是不好选… 这少女的心思也太难猜了吧，于是你大叹了三声：孤为之奈何？孤为之奈何？孤为之奈何？
也许是上天垂帘你的执著，在你苦思冥想了几日不得其解之后，偶然的发现了一个定理：
定理：细致平稳条件 如果非周期马尔可夫链的转移矩阵 $ P $ 和分布$ \pi(x) $ 满足 \(\pi(i)P_{ij}=\pi(j)P_{ji} ~~~ for \ all \ i,j \ in \ X\)
则$\pi(x)$ 是马尔可夫链的平稳分布，上式被称为细致平稳条件 (detailed balance condition)。 凭借你的数学直觉，发现这个定理是很显然的，这里的要求要比稳态向量的定义$qP=q$ 要严格的多，$qP=q$ 只是要求转出到其他所有状态的概率密度等于其他所有状态转入的概率密度。而细致平稳条件则要求针对任意两个状态之间的转入和转出概率密度都相等。
当然，数学证明也是很简单的，由细致平稳条件易得： \(\sum_{i=1}^\infty\pi(i)P_{ij}=\sum_{i=1}^\infty\pi(j)P_{ji}=\pi(j)\sum_{i=1}^\infty P_{ji}=\pi({j})\) 能够推出$\pi$ 是方程$\pi P=\pi$ 的解，所以$\pi$ 是平稳分布。
我们用 $P_{ij}$ 表示转移矩阵$P$ 中的从状态$i$转到状态$j$ 的概率，$\pi$表示稳态向量，用$\pi_i$ 来表示处于状态$i$ 的概率，显然对于一个一般的状态转移矩阵$P$下的稳态向量$\pi$ ,细致平稳条件是不满足的，即 $\pi_iP_{ij}\neq \pi_jP_{ji}$
因此，我们需要对马尔可夫链进行改造，使得细致平稳条件成立。比如我们引入一个$\alpha_{ij}$ ,使得
\[\pi(i)P_{ij}\alpha_{ij}=\pi(j)P_{ji}\alpha_{ji}\]那么$\alpha_{ij}$ 应该取啥呢？根据对称性，我们可以简单的取 \(\alpha_{ij}=\pi(j)P_{ji}\) \(\alpha_{ji}=\pi(i)P_{ij}\) 这样一来上面的等式就成立了。于是我们改造之后的状态转移矩阵$P’=P \alpha$ , 而转移矩阵$P’$ 的稳态向量则是 $\pi$ ,我们将其概率分布表示为$p(x)$.
在构造$P’$ 的过程中，我们引入的$\alpha_{ij}$ 称之为接受率，物理意义可以理解为在原来的马尔可夫链上，从状态$i$, 以$P_{ij}$ 的概率跳转到状态$j$ 的时候，以$\alpha_{ij}$ 的概率接受这个跳转，以$1- \alpha_{ij}$ 的概率拒绝跳转。于是这个接受跳转的部分，以$P_{ij}\alpha_{ij}$ 的概率实现了转移，那么拒绝跳转的部分概率到哪里去了呢？我的理解是拒绝转移的这部分，都留给了$P_{ii}$ ,也就是转移到了当前的状态，只需要用1减去所有跳转到其他状态的概率即可。
整理一下上面的思路，便得到了原始的MCMC采样算法：
- 任意选择一个状态转移矩阵$P$， 平稳分布$\pi(x)$
- 基于任意概率分布采样初始状态$x_0$
- 对$t=0,1,2,…$， 循环以下过程采样： 第t个时刻的马尔可夫链状态为$X_t=x_t$ ,采样$ y \sim P(x | x_t) $ 从均匀分布采样$\mu \sim Uniform[0,1] $ 如果 $ \mu < \alpha_{x_ty}=\pi(y)P(x_t | y)$ 则接受转移$ x_t \rightarrow y $ , 即$ X_{t+1} = y $， 否则不接受转移，即$X_{t+1}=X_t$
仔细观察上面的算法，容易发现这个算法存在一个问题：就是马尔可夫链在转移过程中的接受率$\alpha_{ij}$ 可能偏小，采样过程中拒绝了大量的转移，导致收敛到稳定状态的速度很慢。于是就有人想，有没有什么办法能够提高接受率呢？
根据细致平稳条件， \(\pi(i)P_{ij}\alpha_{ij}=\pi(j)P_{ji}\alpha_{ji}\) 于是， \(\frac{\alpha_{ij} }{\alpha_{ji} }=\frac{\pi(j)P_{ji}}{\pi(i)P_{ij}}\) 可见分子和分母同比例放大或缩小等式恒成立，所以我们只需要让分子和分母中较大的一个值放大到1，就可以将接受率提高到最大，也就是说让 \(\alpha_{ij} = min \left\{\frac{\pi(j)P_{ji}}{\pi(i)P_{ij} },1 \right\}\) 于是，经过改造的MCMC采样算法就变成了教科书中最常见的 Metropolis-Hastings 算法：
- 任意选择一个状态转移矩阵$P$， 平稳分布$\pi(x)$
- 基于任意概率分布采样初始状态$x_0$
- 对$t=0,1,2,…$， 循环以下过程采样： 第t个时刻的马尔可夫链状态为$X_t=x_t$ ,采样$ y \sim P(x | x_t) $ 从均匀分布采样$\mu \sim Uniform[0,1] $ 如果$ \mu < \alpha_{x_ty}=min \left{\frac{\pi(j)P_{ji}}{\pi(i)P_{ij} },1 \right} $ 则接受转移$x_t\rightarrow y $ ，,即$ X_{t+1} = y $ 否则不接受转移，即$X_{t+1}=X_t$
到这里，我们的MCMC把妹法就算是大功告成了，我们可以用它来估计小姐姐心情值的概率分布：
- 任意选择一个少女状态转移函数$P$, 少女心情值概率分布记为$\pi(x)$
- 基于随机观察随便估计一个小姐姐的心情初始值$x_0$
- 然后如此循环一年的采样： 第t个天的马尔可夫链状态为$X_t=x_t$ ,采样$ y \sim P(x | x_t) $ 从均匀分布采样$\mu \sim Uniform[0,1] $ 如果$ \mu < \alpha_{x_ty}=min \left{\frac{\pi(j)P_{ji}}{\pi(i)P_{ij} },1 \right} $ 则接受转移$x_t\rightarrow y $ ,即$ X_{t+1} = y $ 否则不接受转移，即$X_{t+1}=X_t$
显然，上面的值我们都是容易取到计算的，$\pi(x)$也可以每日从小姐姐处获得。然后我们姑且取最后100天的均值方差作为姑娘心情分布的均值和方差大抵也所差无几了，如果假设前200多天已经收敛，那根据概率不等式，精度大概在0.1左右，已经是很好的估计了~
当你成功的运用MCMC把妹法估算出姑娘心情值分布之后，你将结果告诉了她，她计算并对比了一下自己记在自己小本子上的每日心情均值，和你给的差值尽然连1都不到，顿时感动的痛哭流涕，相见恨晚，扑向了你的怀抱，你从此过上了从此君王不早朝的幸福生活。
备注：实际上，上述的M-H采样有个很大的缺点，因为接受率的存在，在高维的情形下不容易收敛，而且一些高维数据的特征的联合分布很不好求，因此一种更好的MCMC采样方法是Gibbs采样,鉴于跟今天所讲的MCMC把妹法关系不大，改日单独再讲，到时把连接放这里… Gibbs采样