吉布斯采样

【吉布斯采样理论等我理解透彻在更新】
吉布斯采样是一种简单并且被广泛使用的马尔科夫蒙特卡洛算法，它可以被看成特殊情况的Metropolis-Hastings算法。
假设我们期望从p(Z)=p(z1,z2,...,zM)中进行采样，并且已经知道初始状态的马尔科夫链。吉布斯采样过程是替换其中一个变量，通过使用剩余变量作为条件对这个变量进行采样。
吉布斯 采样的流程
1、初始化zi:i=1,...,M
2、对于t=1,...T
采样zt+11服从p(z1|zt2,zt3,...,ztM)
采样zt+12服从p(z2|zt+11,zt3,...,ztM)
.
.
采样zt+1M服从p(zM|zt+11,zt+13,...,zt+1M−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;    endendtemp=pix(:,4);temp=temp(2:end);length(find(temp>0))*1.0/(num-1)

python程序

#Gibbsimport numpy as npimport numpy.random as nrnum=100000pix=np.zeros((num,3))pix[0,:]=[0,0,10]temp=0for 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+=1print temp*1.0/num

运行结果是在0.25左右