正規乱数の生成法

正規乱数を生成する方法として，rand( ) の第1引数に Normal(μ, σ) を指定する方法と，randn( ) で標準正規乱数を生成して線形変換する方法がある。

rand( ) による方法

using Distributions
d = Normal(50.0, 10.0)

Normal{Float64}(μ=50.0, σ=10.0)

x = rand(d, 1000000);

using Statistics

mean(x)

49.992122684571775

std(x, corrected=false)

10.011372632539898

randn( ) による方法

y = randn(1000000) * 10.0 .+ 50.0;

mean(y)

49.998784228855804

std(y, corrected=false)

10.001876233453125

速度比較

前述の 2 方法を a, b の関数とし，それぞれを 100 回繰り返して呼びだし，累積実行時間を求める関数 f を書き，実行する。

a(n) = rand(Normal(50.0, 10.0), n)

b(n) = randn(n) * 10.0 .+ 50.0;

function f(FUNC; n=1000000, loop=100)
    for i = 1:loop
        FUNC(n)
    end
end

f (generic function with 1 method)

以下に示すように，randn( ) の結果を線形変換する方法のほうが 1.7 倍ほど速い。

@time f(a)

  1.336329 seconds (200 allocations: 762.947 MiB, 2.74% gc time)

@time f(b)

  0.754293 seconds (600 allocations: 2.235 GiB, 13.41% gc time)

1.336329 / 0.754293

1.7716311831078906

完全に指定した母平均，母標準偏差になる正規乱数を生成する

生成された正規乱数の平均，標準偏差が，指定した母平均，母標準偏差になるようにする。

母平均値=0，母標準偏差=1 を指定しても，生成された標準正規乱数の平均値=0，標準偏差=1 ではない。

そこで，まずこの正規乱数を標準化する。つまり，平均値=0，標準偏差=1 になる。

次に，この標準正規乱数を指定された母平均，母標準偏差になるように線形変換する（標準化の逆）。

function c(n, μ, σ)
    x = randn(n)
    ((x .- mean(x)) ./ std(x, corrected=false)) .* σ .+ μ 
end

c (generic function with 1 method)

x = c(10000, 50, 10);
println("mean = $(mean(x)), sd = $(std(x, corrected=false))")

mean = 50.0, sd = 9.999999999999998

蛇足であるが，関数 c の中で randn( ) を使ったが，ここは rand( ) でも，極端な場合には 1:n でもなんでもよい。それぞれ得られる乱数の分布の形は違うが，指定された母平均，母標準偏差と全く同じ乱数を作ることができる。

n, μ, σ = 100, 50, 10
y = rand(n) # 一様乱数
y = ((y .- mean(y)) ./ std(y, corrected=false)) .* σ .+ μ
println("mean = $(mean(y)), sd = $(std(y, corrected=false))")

mean = 50.0, sd = 9.999999999999998

n, μ, σ = 100, 50, 10
z = 1:n # 等差数列
z = ((z .- mean(z)) ./ std(z, corrected=false)) .* σ .+ μ
println("mean = $(mean(z)), sd = $(std(z, corrected=false))")

mean = 50.0, sd = 9.999999999999996