1

我在运行带有拒绝采样的 Gibbs 采样器的 Matlab 和 Julia（使用 tic、toc 命令）之间进行了简单的速度比较。对于 Matlab，我运行了两个不同版本的代码。一个是使用内置的 unifrnd(a,b) 而另一个是通过调用以下函数在区间 (a,b) 中绘制均匀随机数：

% Draw a uniform random sample in the interval (a,b)
function f = rand_uniform(a,b)
f = a + rand*(b - a);
end

我在 Julia 代码中使用了与上面相同的函数。

以 1000000 次迭代运行代码的结果如下：

案例 1：使用 unifrnd(a,b) 的 Matlab：

x_bar = 1.0944

y_bar = 1.1426

经过的时间是 255.201619 秒。

案例 2：Matlab 调用 rand_uniform(a,b) （上面的函数）：

x_bar =1.0947

y_bar =1.1429

经过的时间是 38.704601 秒。

案例 3：Julia 调用 rand_uniform(a,b)（上面的函数）：

x_bar = 1.0951446303536603

y_bar = 1.142634615899686

经过时间：3.563854193 秒

显然使用 unifrnd(a,b) 会大大降低 Matlab 代码的速度，但问题是为什么？密度取自社会科学家应用贝叶斯统计和估计简介中的一个示例，如果有人对理论平均值感兴趣，mu_x = 1.095 而 mu_y = 1.143。

案例1：使用内置unifrnd(a,b)的Matlab代码为：

clear;
clc;
tic

% == Setting up matrices and preliminaries == %
d = 1000000;                        % No. of iterations                     
b = d*0.25;                         % Burn-in length. We discard 25% of the sample (not used here)
x = zeros(1,d);
x(1) = -1;
y = zeros(1,d);
y(1) = -1;
m = 25;                             % The constant multiplied on g(x)
count = 0;                          % Counting no. of iterations

% == The Gibbs sampler using rejection sampling == %
for i = 2:d

    % Sample from x|y
    z = 0;
    while z == 0
        u = unifrnd(0,2);           
        if ((2*u+3*y(i-1)+2) > unifrnd(0,2)*m) % Height is (1/(b-a))*m = 0.5*25 = 12.5
            x(i) = u;
            z = 1;
        end
    end

    % Sample from y|x
    z = 0;
    while z == 0
        u = unifrnd(0,2);
        if ((2*x(i)+3*u+2) > unifrnd(0,2)*m)
            y(i) = u;
            z = 1;
        end
    end     
    %count = count+1                 % For counting no. of total draws from m*g(x) 
end

x_bar = mean(x)
y_bar = mean(y)

toc

案例 2：调用 rand_uniform(a,b) 的 Matlab 代码（上面的函数）：

clear;
clc;
tic

% == Setting up matrices and preliminaries == %
d = 1000000;                        % No. of iterations                     
b = d*0.25;                         % Burn-in length. We discard 25% of the sample (not used here)
x = zeros(1,d);
x(1) = -1;
y = zeros(1,d);
y(1) = -1;
m = 25;                             % The constant multiplied on g(x)
count = 0;                          % Counting no. of iterations

% == The Gibbs sampler using rejection sampling == %
for i = 2:d

    % Sample from x|y
    z = 0;
    while z == 0
        u = rand_uniform(0,2);           
        if ((2*u+3*y(i-1)+2) > rand_uniform(0,2)*m) % Height is (1/(b-a))*m = 0.5*25 = 12.5
            x(i) = u;
            z = 1;
        end
    end

    % Sample from y|x
    z = 0;
    while z == 0
        u = rand_uniform(0,2);
        if ((2*x(i)+3*u+2) > rand_uniform(0,2)*m)
            y(i) = u;
            z = 1;
        end
    end     
    %count = count+1                 % For counting no. of total draws from m*g(x) 
end

x_bar = mean(x)
y_bar = mean(y)

toc

案例 3：调用 rand_uniform(a,b) 的 Julia 代码（上面的函数）：

# Gibbs sampling with rejection sampling
tic()

# == Return a uniform random sample from the interval (a, b) == #
function rand_uniform(a, b)
    a + rand()*(b - a)
end

# == Setup and preliminaries == #
d = 1000000                         # No. of iterations
b = d*0.25                          # Burn-in length. We discard 25% of the sample (not used here)
x = zeros(d)
x[1] = -1
y = zeros(d)
y[1] = -1
m = 25

# == The Gibbs sampler using rejection sampling == #
for i in 2:d

  #Sample from y|x
  z = 0
  while z==0
    u = rand_uniform(0,2)
    if ((2*u+3*y[i-1]+2) > rand_uniform(0,2)*m)   #Height is (1/(b-a))*m = 0.5*25 = 12.5
      x[i] = u
      z = 1
    end
  end

  #Sample from x|y
  z = 0
    while z == 0
        u = rand_uniform(0,2)
        if ((2*x[i]+3*u+2) > rand_uniform(0,2)*m)
            y[i] = u
            z = 1
        end
    end
end

println("x_bar = ", mean(x))
println("y_bar = ", mean(y))

toc()
4

0 回答 0