matlab 产生任意概率密度的联合分布

用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,生成单个拉普拉斯分布的采样

%%

laplaceDis1=@(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没有关联,所以可以直接生成联合分布

laplaceDis1=@(y)exp(-abs(y))/2; %mu=0,b=1的拉普拉斯分布
laplaceDis2=@(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');

如图: