Random

Julia の乱数生成は MersenneTwister オブジェクトを通してdSFMTというメルセンヌツイスターライブラリを利用します。Julia はグローバルな乱数生成器 (random number generator, RNG) を持ち、デフォルトではこれを使って乱数が生成されます。他の RNG 型は AbstractRNG を継承することで追加でき、そういった型は異なる種類の乱数列を生成できます。MersenneTwister 以外に Julia が提供する RNG 型の一つが RandomDevice であり、これは OS が提供するエントロピープールのラッパーです。

乱数生成に関する関数の多くは AbstractRNG オブジェクトを省略可能な第一引数に受け取ります。この引数を渡さなければグローバルな RNG が使われます。さらに一部の関数では可変長引数 dims... (もしくはタプルの引数 dims) に次元を指定でき、指定すると乱数の配列が生成されます。

MersenneTwister と RandomDevice が生成できる乱数の型は Float16, Float32, Float64, BigFloat, Bool, Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128, BigInt (およびこれらの複素数型) です。浮動小数点数型の乱数は [0,1) の範囲で一様に生成されます。BigInt 型はいくらでも大きい整数を表せるので、乱数を生成するときは rand(big.(1:6)) などとして区間を指定する必要があります。

AbstractFloat と Complex の一部に対しては正規分布や指数分布も実装されます。詳細は randn と randexp を参照してください。

Random モジュール

Random

乱数生成関数

rand([rng=GLOBAL_RNG], [S], [dims...])

julia> rand(Int, 2)
2-element Array{Int64,1}:
 1339893410598768192
 1575814717733606317

julia> using Random

julia> rand(MersenneTwister(0), Dict(1=>2, 3=>4))
1=>2

julia> rand((2, 3))
3

julia> rand(Float64, (2, 3))
2×3 Array{Float64,2}:
 0.999717  0.0143835  0.540787
 0.696556  0.783855   0.938235

rand!([rng=GLOBAL_RNG], A, [S=eltype(A)])

julia> rng = MersenneTwister(1234);

julia> rand!(rng, zeros(5))
5-element Array{Float64,1}:
 0.5908446386657102
 0.7667970365022592
 0.5662374165061859
 0.4600853424625171
 0.7940257103317943

bitrand([rng=GLOBAL_RNG], [dims...])

julia> rng = MersenneTwister(1234);

julia> bitrand(rng, 10)
10-element BitArray{1}:
 0
 1
 1
 1
 1
 0
 1
 0
 0
 1

randn([rng=GLOBAL_RNG], [T=Float64], [dims...])

julia> using Random

julia> rng = MersenneTwister(1234);

julia> randn(rng, ComplexF64)
0.6133070881429037 - 0.6376291670853887im

julia> randn(rng, ComplexF32, (2, 3))
2×3 Array{Complex{Float32},2}:
 -0.349649-0.638457im  0.376756-0.192146im  -0.396334-0.0136413im
  0.611224+1.56403im   0.355204-0.365563im  0.0905552+1.31012im

randn!([rng=GLOBAL_RNG], A::AbstractArray) -> A

julia> rng = MersenneTwister(1234);

julia> randn!(rng, zeros(5))
5-element Array{Float64,1}:
  0.8673472019512456
 -0.9017438158568171
 -0.4944787535042339
 -0.9029142938652416
  0.8644013132535154

randexp([rng=GLOBAL_RNG], [T=Float64], [dims...])

julia> rng = MersenneTwister(1234);

julia> randexp(rng, Float32)
2.4835055f0

julia> randexp(rng, 3, 3)
3×3 Array{Float64,2}:
 1.5167    1.30652   0.344435
 0.604436  2.78029   0.418516
 0.695867  0.693292  0.643644

randexp!([rng=GLOBAL_RNG], A::AbstractArray) -> A

julia> rng = MersenneTwister(1234);

julia> randexp!(rng, zeros(5))
5-element Array{Float64,1}:
 2.4835053723904896
 1.516703605376473
 0.6044364871025417
 0.6958665886385867
 1.3065196315496677

randstring([rng=GLOBAL_RNG], [chars], [len=8])

julia> Random.seed!(3); randstring()
"4zSHdXlw"

julia> randstring(MersenneTwister(3), 'a':'z', 6)
"bzlhqn"

julia> randstring("ACGT")
"AGGACATT"

部分列/置換/シャッフル

randsubseq([rng=GLOBAL_RNG,] A, p) -> Vector

julia> rng = MersenneTwister(1234);

julia> randsubseq(rng, 1:8, 0.3)
2-element Array{Int64,1}:
 7
 8

randsubseq!([rng=GLOBAL_RNG,] S, A, p)

julia> rng = MersenneTwister(1234);

julia> S = Int64[];

julia> randsubseq!(rng, S, 1:8, 0.3)
2-element Array{Int64,1}:
 7
 8

julia> S
2-element Array{Int64,1}:
 7
 8

randperm([rng=GLOBAL_RNG,] n::Integer)

julia> randperm(MersenneTwister(1234), 4)
4-element Array{Int64,1}:
 2
 1
 4
 3

randperm!([rng=GLOBAL_RNG,] A::Array{<:Integer})

julia> randperm!(MersenneTwister(1234), Vector{Int}(undef, 4))
4-element Array{Int64,1}:
 2
 1
 4
 3

randcycle([rng=GLOBAL_RNG,] n::Integer)

julia> randcycle(MersenneTwister(1234), 6)
6-element Array{Int64,1}:
 3
 5
 4
 6
 1
 2

randcycle!([rng=GLOBAL_RNG,] A::Array{<:Integer})

julia> randcycle!(MersenneTwister(1234), Vector{Int}(undef, 6))
6-element Array{Int64,1}:
 3
 5
 4
 6
 1
 2

shuffle([rng=GLOBAL_RNG,] v::AbstractArray)

julia> rng = MersenneTwister(1234);

julia> shuffle(rng, Vector(1:10))
10-element Array{Int64,1}:
  6
  1
 10
  2
  3
  9
  5
  7
  4
  8

shuffle!([rng=GLOBAL_RNG,] v::AbstractArray)

julia> rng = MersenneTwister(1234);

julia> shuffle!(rng, Vector(1:16))
16-element Array{Int64,1}:
  2
 15
  5
 14
  1
  9
 10
  6
 11
  3
 16
  7
  4
 12
  8
 13

RNG の作成とシード

seed!([rng=GLOBAL_RNG], seed) -> rng
seed!([rng=GLOBAL_RNG]) -> rng

julia> Random.seed!(1234);

julia> x1 = rand(2)
2-element Array{Float64,1}:
 0.590845
 0.766797

julia> Random.seed!(1234);

julia> x2 = rand(2)
2-element Array{Float64,1}:
 0.590845
 0.766797

julia> x1 == x2
true

julia> rng = MersenneTwister(1234); rand(rng, 2) == x1
true

julia> MersenneTwister(1) == Random.seed!(rng, 1)
true

julia> rand(Random.seed!(rng), Bool) # not reproducible
true

julia> rand(Random.seed!(rng), Bool)
false

julia> rand(MersenneTwister(), Bool) # not reproducible either
true

AbstractRNG

MersenneTwister(seed)
MersenneTwister()

julia> rng = MersenneTwister(1234);

julia> x1 = rand(rng, 2)
2-element Array{Float64,1}:
 0.5908446386657102
 0.7667970365022592

julia> rng = MersenneTwister(1234);

julia> x2 = rand(rng, 2)
2-element Array{Float64,1}:
 0.5908446386657102
 0.7667970365022592

julia> x1 == x2
true

RandomDevice()

乱数 API に対するフック

Random の機能はほぼ直交する次の二つの方法で拡張できます:

1. 独自型のランダムな値を生成する。
2. 新しい RNG を作成する。

(1) に対する API は非常に関数的ですが、比較的最近に追加された API なので、将来の Random モジュールのリリースで変更が必要になるかもしれません。通常は rand メソッドを一つ実装すれば他のメソッドも自動的に全て使えるようになります。

(2) に対する API は現在でも初歩的であり、通常の型の値を生成するだけでも厳密には必要としてはいけない量の作業が実装者に課される場合があります。

乱数生成のカスタマイズ

何らかの分布に沿う乱数を生成する処理には様々なトレードオフがあります。事前に計算された値を使う方法、例えば離散分布に対する別名法や単変量分布に対する"圧縮" 関数を使った棄却法を使えば、サンプリングを大きく高速化できます。こういった手法では事前に計算すべき情報の量が分布から生成するサンプル数に依存する場合もあれば、RNG が持つ特定の性質をアルゴリズムが利用できるという場合もあります。

Random モジュールはこういった問題に対処しつつ乱数を生成するためのカスタマイズ可能なフレームワークを定義しています。rand が呼び出しのたびに Sampler 型のサンプラーを生成し、Sampler 型に対するメソッドを追加することで上述のトレードオフを念頭に置いたカスタマイズを行うというものです。その後 Sampler 型に対するメソッドは RNG に対してディスパッチして、RNG の持つ分布や反復回数のヒントに関する情報を利用します。反復回数のヒントとして現在使えるのは Val{1} (一度のサンプル) と Val{Inf} (任意回数のサンプル) だけです。この二つの型は Random.Repetition の部分型です。

乱数の生成は Sampler オブジェクトを使って行われます。値 X を使う乱数生成インターフェースを実装するとき、実装すべきメソッドは Sampler(rng, X, repetition) が返す sampler に対する rand(rng, sampler) です。

rand(rng, sampler) を実装してさえいればどんな値でもサンプラーにできますが、ほとんどのアプリケーションでは次の事前に定義されるサンプラーで十分です:

1. SamplerType{T}() は型 T の値を生成するサンプラーの実装で利用されます。これは rand(Int) のような型を使った乱数生成におけるデフォルトの Sampler です。

2. SamplerTrivial(self) は self の簡単なラッパーです。self には [] でアクセスできます。事前の計算が必要ないサンプラーではこのサンプラーが推奨されます。これは rand(1:3) のような値を使った乱数生成におけるデフォルトの Sampler です。

3. SamplerSimple(self, data) は data という追加のフィールドを持つサンプラーです。data は任意の事前に計算された値を格納するために使われます。data の値は Sampler の独自メソッドで計算されるべきです。

後でそれぞれの例を示します。アルゴリズムの選択は RNG と独立していると仮定するので、シグネチャでは AbstractRNG を使います。

Sampler(rng, x, repetition = Val(Inf))

SamplerType{T}()

SamplerTrivial(x)

SamplerSimple(x, data)

事前計算と値の生成の分離は API の一部であり、ユーザーからも利用できます。例えばループの中で rand(rng, 1:20) を何度も呼び出す必要があるときは、次のようなコードを書けば事前計算と値の生成の分離を活用できます:

rng = MersenneTwister()
sp = Random.Sampler(rng, 1:20) # あるいは Random.Sampler(MersenneTwister, 1:20)
for x in X
    n = rand(rng, sp) # n = rand(rng, 1:20) と同じ処理
    # n を使った処理...
end

これは標準ライブラリで使われているメカニズムでもあります。例えば rand(1:20, 10) でランダムな配列生成をするとこのテクニックが使われます。

型を使った値の生成

現在 Julia は型 T に対して rand(T) が定義されているなら、それは T 型のオブジェクトを生成すると仮定します。型を使った乱数生成でフォールバックとして使われるデフォルトのサンプラーは SamplerType です。T 型の値のランダムな生成を定義するには、メソッド rand(rng::AbstractRNG, ::Random.SamplerType{T}) を定義して、そこから rand(rng, T) が返すべき値を返してください。

例として、1 から n までの n 個の目を持つサイコロを表す Die 型に対する rand(Die) を実装します。rand(Die) は目の個数を 4 から 20 の間でランダムに選んだ Die を返すとします:

struct Die
    nsides::Int # 目の個数
end

Random.rand(rng::AbstractRNG, ::Random.SamplerType{Die}) = Die(rand(rng, 4:20))

これで Die に対するスカラーと配列の rand メソッドが動作するようになります:

julia> rand(Die)
Die(10)

julia> rand(MersenneTwister(0), Die)
Die(16)

julia> rand(Die, 3)
3-element Array{Die,1}:
 Die(5)
 Die(20)
 Die(9)

julia> a = Vector{Die}(undef, 3); rand!(a)
3-element Array{Die,1}:
 Die(11)
 Die(20)
 Die(10)

事前に計算されたデータを持たない簡単なサンプラー

続いてコレクションを使った乱数生成で使われるサンプラーを定義します。事前に計算されるデータが無いなら、このサンプラーは値を使った乱数生成におけるデフォルトのフォールバック SamplerTrivial を使って実装できます。

型 S のオブジェクトからのランダムな値の生成を定義するには、メソッド rand(rng::AbstractRNG, sp::Random.SamplerTrivial{S}) の定義が必要です。ここで sp は S 型のオブジェクトをラップしているだけであり、そのオブジェクトには sp[] でアクセスできます。Die の例に戻って、次は rand(d::Die) が d の目からランダムに選んだ Int を返すようにしたいとします。この処理は次のコードで実装できます:

julia> Random.rand(rng::AbstractRNG, d::Random.SamplerTrivial{Die}) =
         rand(rng, 1:d[].nsides);

julia> rand(Die(4))
2

julia> rand(Die(4), 3)
3-element Array{Any,1}:
 1
 4
 2

現在 Julia は S 型のコレクションに対して rand(::S) が定義されているなら、rand(::S) は eltype(S) 型のオブジェクトを生成すると仮定します。最後の例では Vector{Any} が生成されていますが、これは eltype(Die) = Any であるためです。実際のコードでは Base.eltype(::Type{Die}) = Int と定義することが推奨されます。

AbstractFloat 型の値の生成

AbstractFloat 型の乱数生成は特別扱いされており、型が表せる範囲全体ではなく [0,1) の範囲から乱数を取るのがデフォルトになっています。そのため T <: AbstractFloat に対してはメソッド Random.rand(::AbstractRNG, ::Random.SamplerTrivial{Random.CloseOpen01{T}}) を定義するべきです。

事前に計算されたデータを持つ最適化されたサンプラー

与えられた確率に従って 1:n の整数が生成される離散分布を考えます。この分布から多くの値を生成する必要があるときは、別名テーブルを使った方法が高速です。このテーブルを構築するアルゴリズムはここに示しませんが、このアルゴリズムを実装した make_alias_table(probabilities) 関数が定義されていて、さらに別名テーブルを使って値を生成する関数 draw_number(rng, alias_table) が利用できるとします。

離散分布は次の型で表します:

struct DiscreteDistribution{V <: AbstractVector}
    probabilities::V
end

ここでは必要な値の個数に関わらず常に別名テーブルを構築するとします (このカスタマイズ方法は後述します)。この場合は、事前に計算されたデータを持ったサンプラーを返す次のメソッドを定義することになります:

Random.eltype(::Type{<:DiscreteDistribution}) = Int

function Random.Sampler(::Type{<:AbstractRNG},
                        distribution::DiscreteDistribution,
                        ::Repetition)
    SamplerSimple(disribution, make_alias_table(distribution.probabilities))
end

そしてランダムな値を生成するメソッドの定義は次のようになります:

function rand(rng::AbstractRNG, sp::SamplerSimple{<:DiscreteDistribution})
    draw_number(rng, sp.data)
end

独自のサンプラー型

事前に計算されたデータを使う場合はたいてい SamplerSimple 型で十分です。ただ、ここでは独自のサンプラー型を使う方法を示すために、SamplerSimple に似た型を実装します。

Die の例をここでも使います: rand(::Die) は固定された区間から乱数生成を行うので、この部分は最適化が可能です。独自のサンプラー型は SamplerDie と呼ぶことにします:

import Random: Sampler, rand

struct SamplerDie <: Sampler{Int} # 生成されるのは Int 型の値
    die::Die
    sp::Sampler{Int} # 抽象型なので、改善が可能
end

Sampler(RNG::Type{<:AbstractRNG}, die::Die, r::Random.Repetition) =
    SamplerDie(die, Sampler(RNG, 1:die.nsides, r))
# パラメータ r については後述

rand(rng::AbstractRNG, sp::SamplerDie) = rand(rng, sp.sp)

この定義があれば、rng を使う rand の呼び出しに die ではなく sp = Sampler(rng, die) と構築したサンプラー sp を渡せるようになります。この非常に単純な例では die を SamplerDie に格納する必要はありませんが、実際のコードで独自のサンプラー型を書くときは通常必要になります。

もちろん、このパターンは非常に頻繁に表れるので、上述のヘルパー型 Random.SamplerSimple が用意されています。この型を使えば SamplerDie を新しく定義する必要はありません。つまりサンプラーの生成とランダムな値の生成の分離は次のようにも実装できます:

Sampler(RNG::Type{<:AbstractRNG}, die::Die, r::Random.Repetition) =
    SamplerSimple(die, Sampler(RNG, 1:die.nsides, r))

rand(rng::AbstractRNG, sp::SamplerSimple{Die}) = rand(rng, sp.data)

このコードで sp.data が指すのは SamplerSimple コンストラクタの第二パラメータに渡した Sampler(rng, 1:die.nsides, r) です。第一パラメータに渡した Die オブジェクトには sp[] でアクセスできます。

SamplerDie と同様に、任意の独自サンプラー型は Sampler{T} の部分型でなければならず、ここで T は生成される値を表します。なお SamplerSimple(x, data) isa Sampler{eltype(x)} であり、SamplerSimple の第一引数はこの事実で制限されます (上述の Die の例のように、Sampler の定義では x を SamplerSimple の引数としてそのまま使うことが推奨されます)。同様に SamplerTrivial(x) isa Sampler{eltype(x)} が成り立ちます。

現在 Random.SamplerTag というヘルパー型も利用可能ですが、これは内部 API とみなされており、任意の時点で非推奨の警告無しに変更される可能性があります。

生成する値の個数でアルゴリズムを切り替える

ランダムな値を数個だけ生成するのか、それとも大量に生成するのかがアルゴリズムの選択に影響する場合があります。これは Sampler コンストラクタの第三パラメータで処理されます。Die に対するヘルパー型として、少数のランダムな値を生成するのに使う SamplerDie1 と大量のランダムな値を生成する SamplerDieMany の二つを次のように定義したとすれば、この型は次のよう利用できます:

Sampler(RNG::Type{<:AbstractRNG}, die::Die, ::Val{1}) = SamplerDie1(...)
Sampler(RNG::Type{<:AbstractRNG}, die::Die, ::Val{Inf}) = SamplerDieMany(...)

rand は型ごとに (rand(::AbstractRNG, ::SamplerDie1) と rand(::AbstractRNG, ::SamplerDieMany) の) 定義がもちろん必要です。また通常通り、必要なら SamplerTrivial と SamplerSimple を使うこともできます。

注意: Sampler(rng, x) は Sampler(rng, x, Val(Inf)) の省略形、Random.Repetition は Union{Val{1}, Val{Inf}} の別名です。

新しい RNG の作成

新しい RNG を作る API はまだ固まっていませんが、大まかな要件は次の通りです:

1. "基礎的な" 型 (整数あるいは浮動小数点数を表す Base に含まれるビット型) を返す rand メソッドは、必要なら新しい RNG に対して個別に定義されるべきである。
2. 他のドキュメントされた rand メソッドで AbstractRNG を受け取るものは (1. のメソッドが存在するなら) そのまま動作する。新しい RNG 型で最適化が可能なら、もちろん特殊化したメソッドを定義することもできる。
3. RNG の copy はコピー元の RNG がその時点から生成するのと全く同じ乱数列を生成する独立したコピーを返す。これが不可能な場合 (例えばハードウェアベースの RNG を使う場合) には copy を実装してはいけない。

1. に関して、個別に定義していない rand メソッドであっても自動的に動作する可能性もありますが、この振る舞いは公式にサポートされておらず、将来のリリースで変更される場合でも警告は出されません。

例えば MyRNG 型に対して新しい rand メソッドを定義するときは、上述の通り次の二つのメソッドの定義が必要です。ここで S は生成する値を指定する値 s (例えば s == Int や s == 1:10) の型 typeof(s) (s が型なら S == Type{s}) です:

1. Sampler(::Type{MyRNG}, ::S, ::Repetition)
2. rand(rng::MyRNG, sp::SamplerS) (SamplerS は 1. の関数が返すサンプラーの型)

Sampler(rng::AbstractRNG, ::S, ::Repetition) が Random モジュールで既に定義されている可能性もありますが、その場合は 1. をスキップしても動作します (新しい RNG 型を使った乱数生成部分だけの定義で済みます) が、そのとき対応する SamplerS 型は内部の実装詳細とみなされるので、警告無しに変更される場合があります。

配列生成の特殊化

ランダムな値からなる配列を生成するとき、RNG 型によっては前述した事前計算と値の生成の分離する方法より特殊化されたメソッドの方が高速になる場合があります。例えば MersenneTwister はランダムな値をネイティブに配列として生成するので、この特徴を持ちます。

要素型 S を持つ値の指定 s と MyRNG 型の RNG に対する特殊化を実装するには、rand!(rng::MyRNG, a::AbstractArray{S}, ::SamplerS) を実装します。ここで SamplerS は Sampler(MyRNG, s, Val(Inf)) が返すサンプラーの型を表します。AbstractArray の代わりに Array{S} のような部分型に対してだけ機能を実装することもできます。上書きを行わないバージョンの rand はこの特殊化を内部で自動的に呼び出します。

再現性

指定されたシードで初期化された RNG パラメータを使えば、プログラムを複数回実行したときにも同じ疑似乱数列を生成できます。ただし、Julia のマイナーリリース (例えば 1.3 から 1.4 への変更) は同じシードから生成される疑似乱数列の変更を許されています (rand のような低水準関数が生成する列は変わらなくても、randsubseq のような高水準関数でアルゴリズムが更新されて疑似乱数列が変わることもあり得ます)。疑似乱数列が変わらないことを保証すると、アルゴリズム的な改善が不可能になるためです。

ランダムなデータの正確な再現性が必要なら、ランダムなデータを保存してしまうというのも一つの手です。例えば科学的な出版物の補助資料として付けてください。もちろん、Julia とパッケージのバージョンを指定するだけでもある程度の再現性は保証できます。

特定の "ランダムな" データにだけ依存するソフトウェアテストでは、一般にそのデータは保存するかテストコードに埋め込むべきです。一方で、ほとんどのランダムなデータで通過するはずのテスト (例えばランダム行列 A = rand(n,n) に対する A \ (A*x) = x) では、固定されたシードを持つ RNG を使って何度もテストを実行したときに非常に確率の低いデータ (例えば非常に悪条件な行列) によってテストが失敗しないことを確認してしまって構いません。

生成されたランダムなサンプルの統計的な分布が Julia のマイナーリリースで変わらないことは保証されます。

1. 訳注: 複素整数型の情報を追加した。^[return]