Merge pull request #512 from JuliaStats/aa/rmath

Replace references to Rmath with StatsFuns

REQUIRE

@@ -1,6 +1,6 @@
 julia 0.4
 PDMats 0.4.2 0.5
-StatsFuns 0.1.1
+StatsFuns 0.3.0
 Calculus
-StatsBase 0.7.0
+StatsBase 0.8.3
 Compat 0.8.4

src/Distributions.jl

@@ -10,11 +10,6 @@ using Compat
 import Compat.view
 import Base.Random
 import Base: size, eltype, length, full, convert, show, getindex, scale, scale!, rand, rand!
-if VERSION < v"0.5.0-"
-    export shape
-else
-    import Base.shape
-end
 import Base: sum, mean, median, maximum, minimum, quantile, std, var, cov, cor
 import Base: +, -, .+, .-
 import Base.Math.@horner
@@ -241,6 +236,7 @@ export
     sampler,            # create a Sampler object for efficient samples
     scale,              # get the scale parameter
     scale!,             # provide storage for the scale parameter (used in multivariate distribution mvlognormal)
+    shape,              # shape parameter for applicable distributions
     skewness,           # skewness of the distribution
     span,               # the span of the support, e.g. maximum(d) - minimum(d)
     std,                # standard deviation of distribution

src/samplers/binomial.jl

@@ -4,8 +4,7 @@ immutable BinomialRmathSampler <: Sampleable{Univariate,Discrete}
     prob::Float64
 end
 
-rand(s::BinomialRmathSampler) =
-    round(Int, ccall((:rbinom, "libRmath-julia"), Float64, (Float64, Float64), s.n, s.prob))
+rand(s::BinomialRmathSampler) = round(Int, StatsFuns.RFunctions.binomrand(s.n, s.prob))
 
 
 # compute probability vector of a Binomial distribution

src/samplers/gamma.jl

@@ -3,8 +3,7 @@ immutable GammaRmathSampler <: Sampleable{Univariate,Continuous}
     d::Gamma
 end
 
-rand(s::GammaRmathSampler) =
-    ccall((:rgamma, "libRmath-julia"), Float64, (Float64, Float64), shape(s.d), scale(s.d))
+rand(s::GammaRmathSampler) = StatsFuns.RFunctions.gammarand(shape(s.d), scale(s.d))
 
 
 # "Generating gamma variates by a modified rejection technique"

src/samplers/poisson.jl

@@ -3,8 +3,7 @@ immutable PoissonRmathSampler <: Sampleable{Univariate,Discrete}
     mu::Float64
 end
 
-rand(s::PoissonRmathSampler) =
-    round(Int, ccall((:rpois, "libRmath-julia"), Float64, (Float64,), s.mu))
+rand(s::PoissonRmathSampler) = round(Int, StatsFuns.RFunctions.poisrand(s.mu))
 
 
 function poissonpvec(μ::Float64, n::Int)

src/univariate/continuous/beta.jl

@@ -101,7 +101,7 @@ gradlogpdf(d::Beta, x::Float64) =
 
 #### Sampling
 
-rand(d::Beta) = StatsFuns.Rmath.betarand(d.α, d.β)
+rand(d::Beta) = StatsFuns.RFunctions.betarand(d.α, d.β)
 
 
 #### Fit model

src/univariate/continuous/chisq.jl

@@ -73,5 +73,5 @@ gradlogpdf(d::Chisq, x::Float64) =  x > 0.0 ? (d.ν * 0.5 - 1) / x - 0.5 : 0.0
 
 #### Sampling
 
-_chisq_rand(ν::Float64) = StatsFuns.Rmath.chisqrand(ν)
+_chisq_rand(ν::Float64) = StatsFuns.RFunctions.chisqrand(ν)
 rand(d::Chisq) = _chisq_rand(d.ν)

src/univariate/continuous/erlang.jl

@@ -60,4 +60,4 @@ cf(d::Erlang, t::Real)  = (1.0 - im * t * d.θ)^(-d.α)
 
 @_delegate_statsfuns Erlang gamma α θ
 
-rand(d::Erlang) = StatsFuns.Rmath.gammarand(d.α, d.θ)
+rand(d::Erlang) = StatsFuns.RFunctions.gammarand(d.α, d.θ)

src/univariate/continuous/fdist.jl

@@ -85,4 +85,4 @@ end
 
 @_delegate_statsfuns FDist fdist ν1 ν2
 
-rand(d::FDist) = StatsFuns.Rmath.fdistrand(d.ν1, d.ν2)
+rand(d::FDist) = StatsFuns.RFunctions.fdistrand(d.ν1, d.ν2)

src/univariate/continuous/gamma.jl

@@ -84,7 +84,7 @@ cf(d::Gamma, t::Real) = (1.0 - im * t * d.θ)^(-d.α)
 gradlogpdf(d::Gamma, x::Float64) =
     insupport(Gamma, x) ? (d.α - 1.0) / x - 1.0 / d.θ : 0.0
 
-rand(d::Gamma) = StatsFuns.Rmath.gammarand(d.α, d.θ)
+rand(d::Gamma) = StatsFuns.RFunctions.gammarand(d.α, d.θ)
 
 
 #### Fit model

src/univariate/continuous/noncentralchisq.jl

@@ -35,4 +35,4 @@ end
 
 @_delegate_statsfuns NoncentralChisq nchisq ν λ
 
-rand(d::NoncentralChisq) = StatsFuns.Rmath.nchisqrand(d.ν, d.λ)
+rand(d::NoncentralChisq) = StatsFuns.RFunctions.nchisqrand(d.ν, d.λ)

src/univariate/continuous/tdist.jl

@@ -64,7 +64,7 @@ end
 
 @_delegate_statsfuns TDist tdist ν
 
-rand(d::TDist) = StatsFuns.Rmath.tdistrand(d.ν)
+rand(d::TDist) = StatsFuns.RFunctions.tdistrand(d.ν)
 
 function cf(d::TDist, t::Real)
     t == 0 && return complex(1.0)

src/univariate/discrete/binomial.jl

@@ -92,7 +92,7 @@ end
 
 @_delegate_statsfuns Binomial binom n p
 
-rand(d::Binomial) = convert(Int, StatsFuns.Rmath.binomrand(d.n, d.p))
+rand(d::Binomial) = convert(Int, StatsFuns.RFunctions.binomrand(d.n, d.p))
 
 immutable RecursiveBinomProbEvaluator <: RecursiveProbabilityEvaluator
     n::Int

src/univariate/discrete/hypergeometric.jl

@@ -72,7 +72,7 @@ end
 
 @_delegate_statsfuns Hypergeometric hyper ns nf n
 
-rand(d::Hypergeometric) = convert(Int, StatsFuns.Rmath.hyperrand(d.ns, d.nf, d.n))
+rand(d::Hypergeometric) = convert(Int, StatsFuns.RFunctions.hyperrand(d.ns, d.nf, d.n))
 
 immutable RecursiveHypergeomProbEvaluator <: RecursiveProbabilityEvaluator
     ns::Float64

src/univariate/discrete/negativebinomial.jl

@@ -66,7 +66,7 @@ mode(d::NegativeBinomial) = (p = succprob(d); floor(Int,(1.0 - p) * (d.r - 1.) /
 
 @_delegate_statsfuns NegativeBinomial nbinom r p
 
-rand(d::NegativeBinomial) = convert(Int, StatsFuns.Rmath.nbinomrand(d.r, d.p))
+rand(d::NegativeBinomial) = convert(Int, StatsFuns.RFunctions.nbinomrand(d.r, d.p))
 
 immutable RecursiveNegBinomProbEvaluator <: RecursiveProbabilityEvaluator
     r::Float64

src/univariate/discrete/poisson.jl

@@ -82,7 +82,7 @@ end
 
 @_delegate_statsfuns Poisson pois λ
 
-rand(d::Poisson) = convert(Int, StatsFuns.Rmath.poisrand(d.λ))
+rand(d::Poisson) = convert(Int, StatsFuns.RFunctions.poisrand(d.λ))
 
 immutable RecursivePoissonProbEvaluator <: RecursiveProbabilityEvaluator
     λ::Float64

test/mvlognormal.jl

@@ -116,5 +116,5 @@ for (g, μ, Σ) in [
 
     m,s = params(g)
     @test_approx_eq full(m) μ
-    test_mvlognormal(g)
+    test_mvlognormal(g, 10^4)
 end

test/mvnormal.jl

@@ -81,7 +81,7 @@ for (T, g, μ, Σ) in [
     @test isa(g, T)
     @test_approx_eq mean(g) μ
     @test_approx_eq cov(g) Σ
-    test_mvnormal(g)
+    test_mvnormal(g, 10^4)
 
     # conversion between mean form and canonical form
     if isa(g, MvNormal)