随机数

Random number generation in Julia uses the Mersenne Twister library via MersenneTwister objects. Julia has a global RNG, which is used by default. Other RNG types can be plugged in by inheriting the AbstractRNG type; they can then be used to have multiple streams of random numbers. Besides MersenneTwister, Julia also provides the RandomDevice RNG type, which is a wrapper over the OS provided entropy.

大部分与随机数生成相关的函数都接受一个可选的 AbstractRNG 对象作为第一个参数，如果不指定则使用全局默认的。此外，某些函数还接受一个可选的维度参数 dims... (可以是元组）来生成随机数组。

一个 MersenneTwister 或 RandomDevice RNG 能够生成如下类型的随机数： 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 generation functions

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

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!(0); randstring()
"0IPrGg0J"

julia> randstring(MersenneTwister(0), 'a':'z', 6)
"aszvqk"

julia> randstring("ACGT")
"TATCGGTC"

Subsequences, permutations and shuffling

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

julia> rng = MersenneTwister(1234);

julia> randsubseq(rng, collect(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, collect(1:8), 0.3);

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

Generators (creation and seeding)

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

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()

Hooking into the Random API

有两种主要的正交的方式可以扩展 Random 的功能：

1. 为自定义类型生成随机值
2. 创造新的生成器

The API for 1) is quite functional, but is relatively recent so it may still have to evolve in subsequent releases of the Random module. For example, it's typically sufficient to implement one rand method in order to have all other usual methods work automatically.

The API for 2) is still rudimentary, and may require more work than strictly necessary from the implementor, in order to support usual types of generated values.

Generating random values of custom types

There are two categories: generating values from a type (e.g. rand(Int)), or from a collection (e.g. rand(1:3)). The simple cases are explained first, and more advanced usage is presented later. We assume here that the choice of algorithm is independent of the RNG, so we use AbstractRNG in our signatures.

Generating values from a type

Given a type T, it's currently assumed that if rand(T) is defined, an object of type T will be produced. In order to define random generation of values of type T, the following method can be defined: rand(rng::AbstractRNG, ::Random.SamplerType{T}) (this should return what rand(rng, T) is expected to return).

Let's take the following example: we implement a Die type, with a variable number n of sides, numbered from 1 to n. We want rand(Die) to produce a die with a random number of up to 20 sides (and at least 4):

struct Die
    nsides::Int # number of sides
end

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

# output

Die 的标量和数组方法现在可以正常工作了：

julia> rand(Die)
Die(18)

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

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

julia> a = Vector{Die}(undef, 3); rand!(a)
3-element Array{Die,1}:
 Die(18)
 Die(6)
 Die(8)

Generating values from a collection

Given a collection type S, it's currently assumed that if rand(::S) is defined, an object of type eltype(S) will be produced. In order to define random generation out of objects of type S, the following method can be defined: rand(rng::AbstractRNG, sp::Random.SamplerTrivial{S}). Here, sp simply wraps an object of type S, which can be accessed via sp[]. Continuing the Die example, we want now to define rand(d::Die) to produce an Int corresponding to one of d's sides:

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

julia> rand(Die(4))
3

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

In the last example, a Vector{Any} is produced; the reason is that eltype(Die) == Any. The remedy is to define Base.eltype(::Type{Die}) = Int.

Generating values for an AbstractFloat type

AbstractFloat types are special-cased, because by default random values are not produced in the whole type domain, but rather in [0,1). The following method should be implemented for T <: AbstractFloat: Random.rand(::AbstractRNG, ::Random.SamplerTrivial{Random.CloseOpen01{T}})

Optimizing generation with cached computation between calls

When repeatedly generating random values (with the same rand parameters), it happens for some types that the result of a computation is used for each call. In this case, the computation can be decoupled from actually generating the values. This is the case for example with the default implementation for AbstractArray. Assume that rand(rng, 1:20) has to be called repeatedly in a loop: the way to take advantage of this decoupling is as follows:

rng = MersenneTwister()
sp = Random.Sampler(rng, 1:20) # or Random.Sampler(MersenneTwister,1:20)
for x in X
    n = rand(rng, sp) # similar to n = rand(rng, 1:20)
    # use n
end

This mechanism is of course used by the default implementation of random array generation (like in rand(1:20, 10)). In order to implement this decoupling for a custom type, a helper type can be used. Going back to our Die example: rand(::Die) uses random generation from a range, so there is an opportunity for this optimization:

import Random: Sampler, rand

struct SamplerDie <: Sampler{Int} # generates values of type Int
    die::Die
    sp::Sampler{Int} # this is an abstract type, so this could be improved
end

Sampler(RNG::Type{<:AbstractRNG}, die::Die, r::Random.Repetition) =
    SamplerDie(die, Sampler(RNG, 1:die.nsides, r))
# the `r` parameter will be explained later on

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

It's now possible to get a sampler with sp = Sampler(rng, die), and use sp instead of die in any rand call involving rng. In the simplistic example above, die doesn't need to be stored in SamplerDie but this is often the case in practice.

This pattern is so frequent that a helper type named Random.SamplerSimple is available, saving us the definition of SamplerDie: we could have implemented our decoupling with:

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)

Here, sp.data refers to the second parameter in the call to the SamplerSimple constructor (in this case equal to Sampler(rng, 1:die.nsides, r)), while the Die object can be accessed via sp[].

Another helper type is currently available for other cases, Random.SamplerTag, but is considered as internal API, and can break at any time without proper deprecations.

Using distinct algorithms for scalar or array generation

In some cases, whether one wants to generate only a handful of values or a large number of values will have an impact on the choice of algorithm. This is handled with the third parameter of the Sampler constructor. Let's assume we defined two helper types for Die, say SamplerDie1 which should be used to generate only few random values, and SamplerDieMany for many values. We can use those types as follows:

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

Of course, rand must also be defined on those types (i.e. rand(::AbstractRNG, ::SamplerDie1) and rand(::AbstractRNG, ::SamplerDieMany)).

Note: Sampler(rng, x) is simply a shorthand for Sampler(rng, x, Val(Inf)), and Random.Repetition is an alias for Union{Val{1}, Val{Inf}}.

Creating new generators

The API is not clearly defined yet, but as a rule of thumb:

1. any rand method producing "basic" types (isbitstype integer and floating types in Base) should be defined for this specific RNG, if they are needed;
2. other documented rand methods accepting an AbstractRNG should work out of the box, (provided the methods from 1) what are relied on are implemented), but can of course be specialized for this RNG if there is room for optimization.

Concerning 1), a rand method may happen to work automatically, but it's not officially supported and may break without warnings in a subsequent release.

To define a new rand method for an hypothetical MyRNG generator, and a value specification s (e.g. s == Int, or s == 1:10) of type S==typeof(s) or S==Type{s} if s is a type, the same two methods as we saw before must be defined:

1. Sampler(::Type{MyRNG}, ::S, ::Repetition), which returns an object of type say SamplerS
2. rand(rng::MyRNG, sp::SamplerS)

It can happen that Sampler(rng::AbstractRNG, ::S, ::Repetition) is already defined in the Random module. It would then be possible to skip step 1) in practice (if one wants to specialize generation for this particular RNG type), but the corresponding SamplerS type is considered as internal detail, and may be changed without warning.

Specializing array generation

In some cases, for a given RNG type, generating an array of random values can be more efficient with a specialized method than by merely using the decoupling technique explained before. This is for example the case for MersenneTwister, which natively writes random values in an array.

To implement this specialization for MyRNG and for a specification s, producing elements of type S, the following method can be defined: rand!(rng::MyRNG, a::AbstractArray{S}, ::SamplerS), where SamplerS is the type of the sampler returned by Sampler(MyRNG, s, Val(Inf)). Instead of AbstractArray, it's possible to implement the functionality only for a subtype, e.g. Array{S}. The non-mutating array method of rand will automatically call this specialization internally.