MCMC && gibbs sampling

MCMC 蒙特卡洛方法

1）输入马尔可夫链状态转移矩阵P，设定状态转移次数阈值n_1，需要的样本个数n_2
2）从任意简单概率分布采样得到初始状态值x_0
3）for t = 0 to n_1+n_2-1，从条件概率分布P(x|x_t)中采样得到样本x_t+1
样本集(x_n1,x_n1+1,x_n1+2,...,x_n1+n2-1)为我们需要的平稳分布对应的样本集

对于一个随意的平稳分布，如何得到其对应的马尔可夫链状态转移矩阵P呢：
http://www.cnblogs.com/pinard/p/6638955.html
1. MCMC采样
目标矩阵P可以通过任意一个马尔可夫状态转移矩阵Q乘上a(i,j)得到。
a(i,j) = phi(j)Q(j,i)
a(j,i) = phi(i)Q(i,j)
其中，a[i,j]表示为接受率，[0,1]之间，一个概率值，P通过Q以一定的接受率获得
1）输入任意选定的马尔可夫转移矩阵Q，平稳分布phi(x),设定状态转移次数阈值n_1，
需要的样本个数n_2
2）从任意简单概率分布采样得到初始状态值x_0
3）for t = 0 to n1+n2-1:
    a）从条件概率分布Q(x|x_t)中采样得到样本x_*
    b）从均匀分布采样u~uniform[0,1]
    c）如果u<a(x_t,x_*) = phi(x_*)Q(x_*,x_t),则接受转移x_t ——> x_*,
       即就是x_t+1=x_*
    d）否则不接受转移，t = max(t-1,0)//本次转移不成功，重新转移本次
样本集合(x_n1,x_n2,...x_n1+n2-1)为我们需要的平稳分布的对应的样本集合
缺点：
a(x_t,x_*)可能性非常小，导致大部分采样值被拒绝，采样效率很低，可能进行了上百万次
还没有收敛，n_1要非常的大，难以接受。

2. M-H采样
phi(i)Q(i,j)*0.1 = phi(j)Q(j,i)*0.2
扩大5倍，细致平稳条件满足：
phi(i)Q(i,j)*0.5 = phi(j)Q(j,i)*1
接受率进行改进：
a(i,j) = min{(π(j)Q(j,i))/ (π(i)Q(i,j)),1}
M-H采样算法过程：
1）输入任意选定的马尔可夫链转移矩阵Q，平稳分布phi(x)，设定状态转移次数阈值
n1，需要的样本个数n2
2）从任意简单的概率分布采样得到初始的状态值x_0
3）for t=0 to n1+n2-1:
    a）从条件概率分布Q(x|xt)中采样得到样本x_*
    b）从均匀分布采样u~uniform[0,1]
    c）如果u< a(x_t,x_*) = min{(π(j)Q(j,i))/(π(i)Q(i,j)),1}，则接受
    x_t ——> x_*,即x_t+1 = x_*
    d）否则不接受转移，t=max(t-1,0)
(x_n1,x_n2,...x_n1+n2-1)为需要的平稳分布对应的样本集
如果选择的Q是对称的，那么：
α(i,j)=min{π(j)/π(i),1}
M-H总结：
1）数据特征非常多，M-H采样由于计算π(j)Q(j,i)/(π(i)Q(i,j))，在高维计算
时间非常可观，算法效率很低，a(i,j)一般小于1，有时候辛苦计算处理被拒绝
2）特征维度大，很难求出各个特征维度的联合分布，但是可以求出各个特征之间的
条件概率分布。

3. Gibbs采样
http://www.cnblogs.com/pinard/p/6645766.html
二维Gibbs采样

P(A→B)=π(x2^(B)|x1^(1)) if x1^(A)=x1^(B)=x1^(1)
P(A→C)=π(x1^(C)|x2^(1)) ifx2^(A)=x2^(C)=x2^(1)
P(A→D)=0 else
有了这个状态转移矩阵，平面上任意两点E，F，满足细致平稳条件
π(E)P(E→F)=π(F)P(F→E)

1）输入平稳分布phi(x1,x2),设定状态转移次数阈值n_1，需要的样本个数n2
2）随机初始化状态值x1^0和x2^0
3）for t = 0 to n1+n2-1:
    a)从条件分布P(x2|x1^t)中采样得到x2^t+1
    b)从条件分布P(x1|x2^t+1)中采样得到样本x1^t+1
样本集{(x1^n1,x2^n1),(x1^n1+1,x2^n1+1),...,x1^n1+n2-1,x2^n1+n2-1}为得到
的平稳分布，通过坐标轴轮换，采样过程为：
(x1^(1),x2^(1))→(x1^(1),x2^(2))→(x1^(2),x2^(2))→...→(x1^(n1+n2−1),x2^(n1+n2−1))
坐标轴不停的轮换的，坐标轮换不是必须的，我们也可以每次随机选择一个坐标轴
进行采样，常用的Gibbs采样是基于坐标轴轮换的。

多维Gibbs采样
对于一个n维度的分布phi(x_1,x_2,...,x_n)，通过在n个坐标轴上轮换采样，得到新的样本
对于一个轮换到的任意一个坐标轴x_i上的转移，马尔可夫链的转移概率为：
P(xi|x1,x2,...,xi−1,xi+1,...,xn)，即固定n-1个坐标轴，在某一个坐标轴上进行移动
具体过程：
1）输入平稳分布π(x1,x2，...,xn)或者对应的所有特征的条件概率分布，设定转移次数n1，需要
的样本个数为n2
2）随机初始状态值(x1^(0),x2^(0),...,xn^(0))
3）for t=0 to n1+n2-1:
    a)从条件概率分布P(x1|x2^(t),x3^(t),...,xn^(t)) 中采样得到样本x1^t+1
    b)从条件概率分布P(x2|x1^(t+1),x3^(t),x4^(t),...,xn^(t))中采样得到样本x2^t+1
    c)从条件概率分布P(xn|x1^(t+1),x2^(t+1),...,xn-1^(t+1))中采样得到样本 xn^t+1
样本集{(x1^(n1),x2^(n1),...,xn^(n1)),...,(x1^(n1+n2−1),x2^(n1+n2−1),...,xn^(n1+n2−1))}
即为我们需要的平稳分布对应的样本集