1

私は Matlab と Julia の間で簡単な速度比較を行っていました (tic、toc コマンドを使用して) 拒否サンプリングで Gibbs サンプラーを実行していました。Matlab では、2 つの異なるバージョンのコードを実行しました。1 つは組み込みの unifrnd(a,b) を使用し、もう 1 つは次の関数を呼び出して間隔 (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_バー = 1.1426

経過時間は 255.201619 秒です。

ケース 2: rand_uniform(a,b) を呼び出す Matlab (上記の関数):

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 コードの速度が大幅に低下しましたが、問題はその理由です。密度は、社会科学者のための応用ベイジアン統計と推定の紹介の例の 1 つから取得されます。誰かが理論上の平均に興味がある場合は、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