Added support for passing in a AbstractRNG to some of the univariat… (

#603)

* Added support for passing in a `AbstractRNG` to some of the univariate `rand` calls.

* Only for univariate distributions which don't rely on RFunctions for their `rand`.
* Doesn't alter the default behaviour (ie: for distributions which don't define their own `rand`)

TODO:
* Update default behaviour to accept `AbstractRNG` and add error case for those which don't support it.
* Update discreteuniform once the new version of StatsBase is released.
* Update truncated, samplers and multivariate.
* Look into replacing the RFunctions requirement.

* Explicitly import GLOBAL_RNG.

src/Distributions.jl

@@ -15,6 +15,7 @@ import Base: sum, mean, median, maximum, minimum, quantile, std, var, cov, cor
 import Base: +, -
 import Base.Math.@horner
 import Base.LinAlg: Cholesky
+import Base.Random: GLOBAL_RNG
 
 if isdefined(Base, :scale)
     import Base: scale

src/samplers/exponential.jl

@@ -2,10 +2,12 @@ immutable ExponentialSampler <: Sampleable{Univariate,Continuous}
     scale::Float64
 end
 
-rand(s::ExponentialSampler) = s.scale * randexp()
+rand(s::ExponentialSampler) = rand(GLOBAL_RNG, s)
+rand(rng::AbstractRNG, s::ExponentialSampler) = s.scale * randexp(rng)
 
 immutable ExponentialLogUSampler <: Sampleable{Univariate,Continuous}
     scale::Float64
 end
 
-rand(s::ExponentialLogUSampler) = -s.scale * log(rand())
+rand(s::ExponentialLogUSampler) = rand(GLOBAL_RNG, s)
+rand(rng::AbstractRNG, s::ExponentialLogUSampler) = -s.scale * log(rand(rng))

src/univariate/continuous/arcsine.jl

@@ -87,4 +87,5 @@ quantile(d::Arcsine, p::Real) = location(d) + abs2(sin(halfπ * p)) * scale(d)
 
 ### Sampling
 
-rand(d::Arcsine) = quantile(d, rand())
+rand(d::Arcsine) = rand(GLOBAL_RNG, d)
+rand(rng::AbstractRNG, d::Arcsine) = quantile(d, rand(rng))

src/univariate/continuous/exponential.jl

@@ -85,8 +85,8 @@ cf(d::Exponential, t::Real) = 1/(1 - t * im * scale(d))
 
 
 #### Sampling
-
-rand(d::Exponential) = xval(d, randexp())
+rand(d::Exponential) = rand(GLOBAL_RNG, d)
+rand(rng::AbstractRNG, d::Exponential) = xval(d, randexp(rng))
 
 
 #### Fit model

src/univariate/continuous/frechet.jl

@@ -137,4 +137,5 @@ end
 
 ## Sampling
 
-rand(d::Frechet) = d.θ * randexp() ^ (-1 / d.α)
+rand(d::Frechet) = rand(GLOBAL_RNG, d)
+rand(rng::AbstractRNG, d::Frechet) = d.θ * randexp(rng) ^ (-1 / d.α)

src/univariate/continuous/generalizedextremevalue.jl

@@ -256,12 +256,12 @@ ccdf(d::GeneralizedExtremeValue, x::Real) = - expm1(logcdf(d, x))
 
 
 #### Sampling
-
-function rand(d::GeneralizedExtremeValue)
+rand(d::GeneralizedExtremeValue) = rand(GLOBAL_RNG, d)
+function rand(rng::AbstractRNG, d::GeneralizedExtremeValue)
     (μ, σ, ξ) = params(d)
 
     # Generate a Float64 random number uniformly in (0,1].
-    u = 1 - rand()
+    u = 1 - rand(rng)
 
     if abs(ξ) < eps(one(ξ)) # ξ == 0
         rd = - log(- log(u))

src/univariate/continuous/generalizedpareto.jl

@@ -179,9 +179,10 @@ end
 
 #### Sampling
 
-function rand(d::GeneralizedPareto)
+rand(d::GeneralizedPareto) = rand(GLOBAL_RNG, d)
+function rand(rng::AbstractRNG, d::GeneralizedPareto)
     # Generate a Float64 random number uniformly in (0,1].
-    u = 1 - rand()
+    u = 1 - rand(rng)
 
     if abs(d.ξ) < eps()
         rd = -log(u)

src/univariate/continuous/gumbel.jl

@@ -92,4 +92,5 @@ gradlogpdf(d::Gumbel, x::Real) = - (1 + exp((d.μ - x) / d.θ)) / d.θ
 
 #### Sampling
 
-rand(d::Gumbel) = quantile(d, rand())
+rand(d::Gumbel) = rand(GLOBAL_RNG, d)
+rand(rng::AbstractRNG, d::Gumbel) = quantile(d, rand(rng))

src/univariate/continuous/inversegaussian.jl

@@ -148,13 +148,14 @@ end
 #   John R. Michael, William R. Schucany and Roy W. Haas (1976)
 #   Generating Random Variates Using Transformations with Multiple Roots
 #   The American Statistician , Vol. 30, No. 2, pp. 88-90
-function rand(d::InverseGaussian)
+rand(d::InverseGaussian) = rand(GLOBAL_RNG, d)
+function rand(rng::AbstractRNG, d::InverseGaussian)
     μ, λ = params(d)
-    z = randn()
+    z = randn(rng)
     v = z * z
     w = μ * v
     x1 = μ + μ / (2λ) * (w - sqrt(w * (4λ + w)))
     p1 = μ / (μ + x1)
-    u = rand()
+    u = rand(rng)
     u >= p1 ? μ^2 / x1 : x1
 end

src/univariate/continuous/kolmogorov.jl

@@ -98,20 +98,21 @@ end
 # Alternating series method, from:
 #   Devroye, Luc (1986) "Non-Uniform Random Variate Generation"
 #   Chapter IV.5, pp. 163-165.
-function rand(d::Kolmogorov)
+rand(d::Kolmogorov) = rand(GLOBAL_RNG, d)
+function rand(rng::AbstractRNG, d::Kolmogorov)
     t = 0.75
-    if rand() < 0.3728329582237386 # cdf(d,t)
+    if rand(rng) < 0.3728329582237386 # cdf(d,t)
         # left interval
         while true
-            g = rand_trunc_gamma()
+            g = rand_trunc_gamma(rng)
 
             x = pi/sqrt(8g)
             w = 0.0
             z = 1/(2g)
             p = exp(-g)
             n = 1
             q = 1.0
-            u = rand()
+            u = rand(rng)
             while u >= w
                 w += z*q
                 if u >= w
@@ -125,8 +126,8 @@ function rand(d::Kolmogorov)
         end
     else
         while true
-            e = randexp()
-            u = rand()
+            e = randexp(rng)
+            u = rand(rng)
             x = sqrt(t*t+e/2)
             w = 0.0
             n = 1
@@ -146,11 +147,11 @@ end
 
 # equivalent to
 # rand(Truncated(Gamma(1.5,1),tp,Inf))
-function rand_trunc_gamma()
+function rand_trunc_gamma(rng::AbstractRNG)
     tp = 2.193245422464302 #pi^2/(8*t^2)
     while true
-        e0 = rand(Exponential(1.2952909208355123))
-        e1 = rand(Exponential(2))
+        e0 = rand(rng, Exponential(1.2952909208355123))
+        e1 = rand(rng, Exponential(2))
         g = tp + e0
         if (e0*e0 <= tp*e1*(g+tp)) || (g/tp - 1 - log(g/tp) <= e1)
             return g

src/univariate/continuous/laplace.jl

@@ -106,7 +106,8 @@ end
 
 #### Sampling
 
-rand(d::Laplace) = d.μ + d.θ*randexp()*ifelse(rand(Bool), 1, -1)
+rand(d::Laplace) = rand(GLOBAL_RNG, d)
+rand(rng::AbstractRNG, d::Laplace) = d.μ + d.θ*randexp(rng)*ifelse(rand(rng, Bool), 1, -1)
 
 
 #### Fitting

src/univariate/continuous/levy.jl

@@ -90,4 +90,5 @@ end
 
 #### Sampling
 
-rand(d::Levy) = d.μ + d.σ / randn()^2
+rand(d::Levy) = rand(GLOBAL_RNG, d)
+rand(rng::AbstractRNG, d::Levy) = d.μ + d.σ / randn(rng)^2

src/univariate/continuous/logistic.jl

@@ -101,4 +101,5 @@ end
 
 #### Sampling
 
-rand(d::Logistic) = quantile(d, rand())
+rand(d::Logistic) = rand(GLOBAL_RNG, d)
+rand(rng::AbstractRNG, d::Logistic) = quantile(d, rand(rng))

src/univariate/continuous/lognormal.jl

@@ -117,7 +117,8 @@ end
 
 #### Sampling
 
-rand(d::LogNormal) = exp(randn() * d.σ + d.μ)
+rand(d::LogNormal) = rand(GLOBAL_RNG, d)
+rand(rng::AbstractRNG, d::LogNormal) = exp(randn(rng) * d.σ + d.μ)
 
 ## Fitting
 

src/univariate/continuous/normal.jl

@@ -78,7 +78,8 @@ cf(d::Normal, t::Real) = exp(im * t * d.μ - d.σ^2/2 * t^2)
 
 #### Sampling
 
-rand(d::Normal) = d.μ + d.σ * randn()
+rand(d::Normal) = rand(GLOBAL_RNG, d)
+rand(rng::AbstractRNG, d::Normal) = d.μ + d.σ * randn(rng)
 
 
 #### Fitting

src/univariate/continuous/normalcanon.jl

@@ -71,5 +71,7 @@ invlogccdf(d::NormalCanon, lp::Real) = xval(d, norminvlogccdf(lp))
 
 #### Sampling
 
-rand(cf::NormalCanon) = cf.μ + randn() / sqrt(cf.λ)
-rand!{T<:Real}(cf::NormalCanon, r::AbstractArray{T}) = rand!(convert(Normal, cf), r)
+rand(cf::NormalCanon) = rand(GLOBAL_RNG, cf)
+rand(rng::AbstractRNG, cf::NormalCanon) = cf.μ + randn(rng) / sqrt(cf.λ)
+rand!{T<:Real}(cf::NormalCanon, r::AbstractArray{T}) = rand!(GLOBAL_RNG, cf, r)
+rand!{T<:Real}(rng::AbstractRNG, cf::NormalCanon, r::AbstractArray{T}) = rand!(rng, convert(Normal, cf), r)

src/univariate/continuous/pareto.jl

@@ -109,7 +109,8 @@ quantile(d::Pareto, p::Real) = cquantile(d, 1 - p)
 
 #### Sampling
 
-rand(d::Pareto) = d.θ * exp(randexp() / d.α)
+rand(d::Pareto) = rand(GLOBAL_RNG, d)
+rand(rng::AbstractRNG, d::Pareto) = d.θ * exp(randexp(rng) / d.α)
 
 
 ## Fitting

src/univariate/continuous/rayleigh.jl

@@ -83,4 +83,5 @@ quantile(d::Rayleigh, p::Real) = sqrt(-2d.σ^2 * log1p(-p))
 
 #### Sampling
 
-rand(d::Rayleigh) = d.σ * sqrt(2 * randexp())
+rand(d::Rayleigh) = rand(GLOBAL_RNG, d)
+rand(rng::AbstractRNG, d::Rayleigh) = d.σ * sqrt(2 * randexp(rng))

src/univariate/continuous/symtriangular.jl

@@ -136,4 +136,5 @@ end
 
 #### Sampling
 
-rand(d::SymTriangularDist) = xval(d, rand() - rand())
+rand(d::SymTriangularDist) = rand(GLOBAL_RNG, d)
+rand(rng::AbstractRNG, d::SymTriangularDist) = xval(d, rand(rng) - rand(rng))

src/univariate/continuous/triangular.jl

@@ -143,10 +143,11 @@ end
 
 #### Sampling
 
-function rand(d::TriangularDist)
+rand(d::TriangularDist) = rand(GLOBAL_RNG, d)
+function rand(rng::AbstractRNG, d::TriangularDist)
     (a, b, c) = params(d)
     b_m_a = b - a
-    u = rand()
+    u = rand(rng)
     b_m_a * u < (c - a) ? d.a + sqrt(u * b_m_a * (c - a)) :
                           d.b - sqrt((1 - u) * b_m_a * (b - c))
 end

src/univariate/continuous/uniform.jl

@@ -103,7 +103,8 @@ end
 
 #### Evaluation
 
-rand(d::Uniform) = d.a + (d.b - d.a) * rand()
+rand(d::Uniform) = rand(GLOBAL_RNG, d)
+rand(rng::AbstractRNG, d::Uniform) = d.a + (d.b - d.a) * rand(rng)
 
 
 #### Fitting

src/univariate/continuous/weibull.jl

@@ -133,4 +133,5 @@ end
 
 #### Sampling
 
-rand(d::Weibull) = xv(d, randexp())
+rand(d::Weibull) = rand(GLOBAL_RNG, d)
+rand(rng::AbstractRNG, d::Weibull) = xv(d, randexp(rng))

src/univariate/discrete/bernoulli.jl

@@ -105,7 +105,8 @@ cf(d::Bernoulli, t::Real) = failprob(d) + succprob(d) * cis(t)
 
 #### Sampling
 
-rand(d::Bernoulli) = round(Int, rand() <= succprob(d))
+rand(d::Bernoulli) = rand(GLOBAL_RNG, d)
+rand(rng::AbstractRNG, d::Bernoulli) = round(Int, rand(rng) <= succprob(d))
 
 
 #### MLE fitting

src/univariate/discrete/geometric.jl

@@ -139,7 +139,8 @@ end
 
 ### Sampling
 
-rand(d::Geometric) = floor(Int,-randexp() / log1p(-d.p))
+rand(d::Geometric) = rand(GLOBAL_RNG, d)
+rand(rng::AbstractRNG, d::Geometric) = floor(Int,-randexp(rng) / log1p(-d.p))
 
 
 ### Model Fitting