luggage baggage

Machine learning, data analysis, web technologies and things around me.

『モンテカルロ統計計算』章末問題を Julia で解く：第2章「乱数」

julialang データ分析 ベイズ 確率・統計 機械学習 MCMC

Julia 言語と MCMC の勉強のため、R 向けに書かれている教科書『モンテカルロ統計計算』（講談社、鎌谷研吾 著、駒木文保 編）の章末問題を解いていこうかと思います。また、単に問題解くだけでは面白くなかったりもするので、適宜 Python との比較などを追加しているものもあります。

今回は、第2章「乱数」です。

2.1 一様乱数のプロット
using Plots

function scatter_2d(x)
    reshaped = reshape(x, (2, length(x) ÷ 2))'
    scatter(reshaped[:, 1], reshaped[:, 2])
end

x = rand(1000)
scatter_2d(x)
2.2 線形合同法による一様乱数の生成
function linear_congruential_generator(a, b, y, n, n_samples)
    x = zeros(n_samples)
    for i in 1:n_samples
        y = (a * y + b) % n
        x[i] = y / n
    end
    x
end

a = 13
b = 0
n = 67
y_init = 1234
n_samples = 1000
x = linear_congruential_generator(a, b, y_init, n, n_samples)
scatter_2d(x)
savefig("2.2.png")

f:id:yoshidabenjiro:20200617232815p:plain

2.3 逆変換法によるロジスティック乱数の生成
using BenchmarkTools, Distributions

function gen_logistic_variables(n_samples)
    x = -log.(1 ./ rand(n_samples) .- 1)
end

n_samples = 1000
dlogistic = Logistic(0, 1)
a = @benchmark gen_logistic_variables(n_samples)
b = @benchmark rand(dlogistic, n_samples)
println(a)
println(b)

# Trial(10.565 μs)
# Trial(17.482 μs)

ちなみに全く同じ実装を Python (Numpy) でやると：

In [1]: import numpy as np

In [2]: def gen_logistic_variables(n_samples):
   ...:     return -np.log(1 / np.random.rand(n_samples) - 1)
   ...:

In [3]: %timeit gen_logistic_variables(1000)
14.5 µs ± 55.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

となり、逆変換法を Julia で実装したもののほうが速い。また、np.random.logistic についても、

In [6]: %timeit np.random.logistic(0, 1, 1000)
17.6 µs ± 74.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

となり、手元の環境では Julia の逆変換法による実装が最速ということになった。なお numpy は conda 経由で入れており MKL にリンクされている。

2.4 逆変換法によるパレート乱数の生成

パレート分布の確率密度関数が p(x;a,b)=\frac{ab^a}{x^{a+1}} なので、b から x までの累積分布関数は F(x;a,b)=1-(\frac{b}{x})^a となる。したがって逆変換法は、 u を一様乱数として x=b(1-u)^{-\frac{1}{a}} とすればよい。コードは以下の通り：

function gen_pareto_variables(a, b, n_samples)
    u = rand(n_samples)
    x = similar(u)
    @. x = b * u ^ (-1 / a)
end
2.5 逆変換法によるレイリー乱数の生成

レイリー分布の確率密度関数が p(x;\sigma)=\frac{x}{\sigma^2}\exp(-\frac{x^2}{2\sigma^2}) なので、 x\ge 0 での累積分布関数は F(x;\sigma) = 1 - \exp (-\frac{x^2}{2\sigma^2}) となる。したがって逆変換法は、 u を一様乱数として x = \sqrt{(-2\sigma^2\ln(1-u))} とすればよい。コードは以下の通り：

function gen_rayleigh_variables(σ, n_samples)
    u = rand(n_samples)
    x = similar(u)
    @. x = sqrt(-2σ^2 * log(u))
end
2.6 逆変換法によるワイブル乱数の生成

ワイブル分布の確率密度関数が p(x;a,b)=\frac{a}{b^a}x^{a-1}\exp( (-\frac{x}{b})^{a}) なので、 x\ge 0 での累積分布関数は F(x;a,b)=1-\exp(-(\frac{x}{b})^a) となる。したがって逆変換は、 u を一様乱数として x=b(-\ln(1-u))^{1/a} とすればよい。コードは以下の通り：

function gen_weibull_variables(a, b, n_samples)
    u = rand(n_samples)
    x = similar(u)
    @. x = b * (-log(u)) ^ (1 / a)
end
2.7 逆変換法による二項乱数の生成
using Distributions, Plots

function gen_binomial_variables(N, θ, n_samples)
    x = zeros(Int, n_samples)
    for i in 1:n_samples
        u = rand()
        n = 0
        F = (1 - θ) ^ N
        while u > F
            n = n + 1
            F = F + binomial(N, n) * θ ^ n * (1 - θ) ^ (N - n)
        end
        x[i] = n
    end
    x
end

N = 10
θ = 0.3
n_samples = 1000
dbinom = Binomial(N, θ)
a = gen_binomial_variables(N, θ, n_samples)
b = rand(dbinom, n_samples)
histogram(a, label = "inversion")
histogram!(b, label = "Distributions", alpha = 0.5)
savefig("2.7.png")

f:id:yoshidabenjiro:20200618234656p:plain

2.8 負の二項分布の生成

次の手順により負の二項分布に従う乱数を生成する：

using Distributions, Plots

function gen_negative_binomial(r, θ, n_samples)
    dgamma = Gamma(r, (1 - θ) / θ)
    x = zeros(n_samples)
    @. x = rand(Poisson(rand(dgamma)))
end

r = 10
θ = 0.5
n_samples = 1000
x = gen_negative_binomial(r, θ, n_samples)
y = rand(NegativeBinomial(r, θ), n_samples)
histogram(x, normalized = true, label = "implemented")
histogram!(y, normalized = true, label = "Distributions", alpha = 0.5)
savefig("2.8.png")

f:id:yoshidabenjiro:20200619005225p:plain
この方法でよい理由は次の通り：上記の式は P(x|y)=\frac{y^x}{x!}e^{-y} および P(y)=\frac{\beta^r}{\Gamma(r)}y^{r-1}e^{-\beta y} を意味するので（ \beta=\frac{\theta}{1-\theta} とした）、 P(x)=\int_0^\infty dyP(x|y)P(y)=\frac{(x+r-1)!}{x!r!}\theta^r(1-\theta)^x となり負の二項分布が得られる。

2.9 TBA
2.10 ラプラス分布を提案分布とする標準正規乱数の棄却法による生成

r(x) = \frac{1}{\lambda}\sqrt{\frac{2}{\pi}}\exp(\frac{-x^2}{2}+\lambda|x|) なので、対数をとり微分した式=0を x の符号に気をつけて解くと |x|=\lambda が得られる。したがって R=r(|x|)=\frac{1}{\lambda}\sqrt{\frac{2}{\pi}}e^{\frac{\lambda^2}{2}} となる。さらに R の対数をとり \lambda について微分したものをゼロとおくと \lambda^2=1 となるが、 \lambda>0 のため \lambda=1 となる。結局、棄却法で用いる定数の最小値は R=\sqrt{\frac{2e}{\pi}} である。

2.11 TBA