你去 - Gibbs Sampler：燃烧期是为了在样本中达到一定的平稳性

burning_period=5000
iterations=1000
y=matrix(nrow=(burning_period+iterations),ncol=3)
a=matrix(nrow=iterations,ncol=3)

y[1,1]=0.5        #Initial Sample 
y[1,2]=0.6
y[1,3]=0.2   

 for(i in 2:(burning_period+iterations)){

     #I have put 3,3,4 as my theta's. You may make the code generic for any choice of theta's.
     t= 1+3*y[i-1,2]+4*y[i-1,3]
     # Use t or t*-1 based on your requirement.
     t=rexp(1,t)  

     y[i,1]=t
     t= 1+3*y[i,1]+5*y[i-1,3]   
     y[i,2]=rexp(1,t)

     t=1+4*y[i,1]+5*y[i,2]
     y[i,3]=rexp(1,t)  
    }
posteriorSample=y[(burning_period+1) : (burning_period+iterations), ]

完整的条件分布可以通过固定两个“其他”变量的值来找到，然后结合所有可以结合的项并查看你得到的结果：

$p(y_1|y_2,y_3) \propto \exp\{-(1+\theta_{12}y_2+\theta_{13}y_3)y_1\}$

这显然是参数为的指数分布，对于也是类似的，它是参数为和，它是参数为的指数分布。 $1+\theta_{12}y_2+\theta_{13}y_3$$p(y_2|y_1,y_3)$$1 + \theta_{12}y_1 + \theta_{23}y_3$$p(y_3|y_1,y_2)$$1 + \theta_{13}y_1 + \theta_{23}y_2$

的一些值开始，然后一次又一次地循环三个条件分布，从适当的分布生成指数变量，并用当前值替换适当的。自然地，您必须丢弃一些初始变量块以消除老化效果。$y_1, y_2, y_3$$y_i$