[MCSM] Slice Sampler

    之前介绍的MCMC算法都具有一般性和通用性（这里指Metropolis-Hasting 算法），但也存在一些特殊的依赖于仿真分布特征的MCMC方法。在介绍这一类算法（指Gibbs sampling）之前，本节将介绍一种特殊的MCMC算法。 我们重新考虑了仿真的理论基础，建立了Slice Sampler。

    考虑到[MCSM]伪随机数和伪随机数生成器中提到的产生服从f(x)密度分布随机数等价于在子图f上产生均匀分布，即

    类似笔记“[MCSM] Metropolis-Hastings 算法”（文章还没写好），考虑采用马尔可夫链的稳态分布来等价上的均匀分布，以此作为f分布的近似。很自然的想法是采用随机行走（random walk）。这样得到的稳态分布是在集合上的均匀分布。

2. 2D slice sample

    有很多方法实现在集合上的"random walk"，最简单的就是一次改变一个方向上的取值，每个方向的改变交替进行，由此得到的算法是 2D slice sampler

---

    在第t次迭代中，执行

    1.

    2. ， 其中

---

    举例

    其中，是归一化因子，代码如下，第一幅图是前10个点的变化轨迹，第二幅图表明初始点的选取影响不大

% p324
T = 0:10000;
T = T/10000;
% N(3,1)
y = exp(-(T+3).^2/2);
plot(T,y);
hold on;
x = 0.25;
u = rand *(exp(-(x+3).^2/2));
x_s = [x];
u_s = [u];
for k = 1:10;
    limit = -3 + sqrt(-2*log(u));
    limit = min([limit 1]);
    x = rand * limit;
    x_s = [x_s x];
    u_s = [u_s u];
    u = rand *(exp(-(x+3).^2/2));
    x_s = [x_s x];
    u_s = [u_s u];
end
plot(x_s,u_s,'-*');
hold off;

%%
x = 0.01;
u = 0.01;
x_s = [x];
u_s = [u];
for k = 1:50;
    limit = -3 + sqrt(-2*log(u));
    limit = min([limit 1]);
    x = rand * limit;
    x_s = [x_s x];
    u_s = [u_s u];
    u = rand *(exp(-(x+3).^2/2));
    x_s = [x_s x];
    u_s = [u_s u];
end
figure;
subplot(1,3,1);
plot(x_s,u_s,'*');hold on;plot(T,y);
x = 0.99;
u = 0.0001;
x_s = [x];
u_s = [u];
for k = 1:50;
    limit = -3 + sqrt(-2*log(u));
    limit = min([limit 1]);
    x = rand * limit;
    x_s = [x_s x];
    u_s = [u_s u];
    u = rand *(exp(-(x+3).^2/2));
    x_s = [x_s x];
    u_s = [u_s u];
end
subplot(1,3,2);
plot(x_s,u_s,'*');hold on;plot(T,y);
x = 0.25;
u = 0.0025;
x_s = [x];
u_s = [u];
for k = 1:50;
    limit = -3 + sqrt(-2*log(u));
    limit = min([limit 1]);
    x = rand * limit;
    x_s = [x_s x];
    u_s = [u_s u];
    u = rand *(exp(-(x+3).^2/2));
    x_s = [x_s x];
    u_s = [u_s u];
end
subplot(1,3,3);
plot(x_s,u_s,'*');hold on;plot(T,y);

View Code