用Slice sampling 來得到某單個分佈的樣本演算法:

已知分佈f(x)

1.      選擇一個初始點，x₀.

2.      在[0,f(x0)] 之間，通過均分分佈產生一個取樣y

3.      在 f(x) 上畫一條通過 y的 水平線

4.      線上段內再次均勻取樣得到一個點x ,此x就是此次迭代得到的新取樣點。

5.      最後用新得的x作為初始值,重複step 2~step5

Y截得f(x)兩點(L,R) ,稱為f(x)的切片，經過不斷迭代，得到很多切片，其中，切片的長度越長，被選取到的概率越大。

難點在於：對於多峰函式，切片所在的邊界往往是不連續的，解決方法是拒絕更短的切片。

Matlab 程式1,生成單個拉普拉斯分佈的取樣

%%

[email protected](y)exp(-abs(y))/2; % mu=0,b=1的拉斯分佈

N=2000; %產生N個取樣點

x = slicesample(1,N,'pdf',laplaceDis1,'thin',5,'burnin',1000);%thin, burn-in 具體設定看介紹

histogram(x,50);%顯示取樣點的直方圖

如圖:

%%matlab 程式2

%% 生成拉普拉斯分佈的聯合分佈，每個隨機變數成分都是獨立的(方便演示)

%  類似Gibbs 多變數取樣的方法，對某個確定的x1,在條件概率密度下P(x2|x1)得到x2的取樣,然後對此x2,得到P(x1|x2)的一個取樣

% 由於此處x1,x2獨立，p(x2|x1)=p(x2),且slicesampling 得到的對概率密度積分和是否為1沒有關聯,所以可以直接生成聯合分佈

[email protected](y)exp(-abs(y))/2; %mu=0,b=1的拉普拉斯分佈
[email protected](y)exp(-sqrt(2)*abs(y-2))/sqrt(2);%mu=2,b=1/sqrt(2) 拉普拉斯分佈


limitL=-5;% 概率密度函式左邊界
limitR=5;

x=zeros(N,2);


for i=1:N
    start=unifrnd(limitL,limitR,1,2);%均布產生一個初值
    x(i,1)=slicesample(start(1),1,'pdf',laplaceDis1);
    x(i,2)=slicesample(start(2),1,'pdf',laplaceDis2);
end

figure
plot(x(:,1),x(:,2),'ro');

如圖: