吉布斯采样

【吉布斯采样理论等我理解透彻在更新】
吉布斯采样是一种简单并且被广泛使用的马尔科夫蒙特卡洛算法，它可以被看成特殊情况的Metropolis-Hastings算法。
假设我们期望从$p(Z)=p(z_1,z_2,...,z_M)$中进行采样，并且已经知道初始状态的马尔科夫链。吉布斯采样过程是替换其中一个变量，通过使用剩余变量作为条件对这个变量进行采样。
吉布斯 采样的流程
1、初始化${z_i : i=1,..., M}$
2、对于$t=1,...T$
采样$z_1^{t+1}$服从$p(z_1|z_2^t,z_3^t,...,z_M^t)$
采样$z_2^{t+1}$服从$p(z_2|z_1^{t+1},z_3^t,...,z_M^t)$
.
.
采样$z_M^{t+1}$服从$p(z_M|z_1^{t+1},z_3^{t+1},...,z_{M-1}^{t+1})$

解决问题，在一个绳子（假设为10米）上面剪两刀，求能构成三角形的概率。
吉布斯采样的大概流程是不知道联合概率分布，只知道每一个分量的条件概率分布。
在这个题目里面，条件概率很简单了。然后一次根据前面状态的分量来采样当前状态的分量。

matlab程序

num = 100000;
pix=zeros(num,4);
pix(1,:) = [0,0,10,0];
for i = 2:num
    for j = 1:2
        if j == 1
            pix(i,1) = rand(1)*(10-pix(i-1,2));
        else
            pix(i,2) = rand(1)*(10-pix(i,1));
        end 
    end
    pix(i,3) = 10-pix(i,1)-pix(i,2);
    if pix(i,1)+pix(i,2) > pix(i,3) && pix(i,3)+pix(i,2) > pix(i,1) && pix(i,1)+pix(i,3) > pix(i,2)
        pix(i,4) = 1;
    else
        pix(i,4) = 0;
    end
end
temp=pix(:,4);
temp=temp(2:end);
length(find(temp>0))*1.0/(num-1)

python程序

#Gibbs
import numpy as np
import numpy.random as nr

num=100000
pix=np.zeros((num,3))
pix[0,:]=[0,0,10]
temp=0
for idx in range(1,num):
    for j in range(0,2):
        if j==0:
            pix[idx,0]=nr.uniform(0,10-pix[idx-1,1],1)
        else:
            pix[idx,1]=nr.uniform(0,10-pix[idx,0],1)
    pix[idx,2]=10-pix[idx,0]-pix[idx,1]
    if pix[idx,0]+pix[idx,1]>pix[idx,2] and pix[idx,2]+pix[idx,1]>pix[idx,0] and pix[idx,0]+pix[idx,2]>pix[idx,1]:
        temp+=1
print temp*1.0/num

运行结果是在0.25左右