Merge branch 'master' into dh/paramsym

src/univariate/continuous/arcsine.jl

@@ -3,7 +3,7 @@ immutable Arcsine <: ContinuousUnivariateDistribution
     b::Float64
 
     function Arcsine(a::Float64, b::Float64)
-        a < b || error("a must be less than b.")
+        a < b || throw(ArgumentError("a must be less than b."))
         new(a, b)
     end
 
@@ -51,4 +51,3 @@ quantile(d::Arcsine, p::Float64) = location(d) + abs2(sin(halfπ * p)) * scale(d
 ### Sampling
 
 rand(d::Arcsine) = quantile(d, rand())
-

src/univariate/continuous/beta.jl

@@ -2,9 +2,10 @@ immutable Beta <: ContinuousUnivariateDistribution
     α::Float64
     β::Float64
 
-    function Beta(a::Real, b::Real)
-        (a > zero(a) && b > zero(b)) || error("α and β must be positive")
-        @compat new(Float64(a), Float64(b))
+    function Beta(α::Real, β::Real)
+        (α > zero(α) && β > zero(β)) ||
+            throw(ArgumentError("α and β must be positive"))
+        @compat new(Float64(α), Float64(β))
     end
 
     Beta(α::Real) = Beta(α, α)
@@ -13,7 +14,6 @@ end
 
 @distr_support Beta 0.0 1.0
 
-
 #### Parameters
 
 params(d::Beta) = (d.α, d.β)

src/univariate/continuous/betaprime.jl

@@ -1,18 +1,21 @@
 immutable BetaPrime <: ContinuousUnivariateDistribution
-    betad::Beta
+    α::Float64
+    β::Float64
 
-    BetaPrime(α::Float64, β::Float64) = new(Beta(α, β))
-    BetaPrime(α::Float64) = new(Beta(α))
-    BetaPrime() = new(Beta())
+    function BetaPrime(α::Float64, β::Float64)
+        (α > zero(α) && β > zero(β)) || error("α and β must be positive")
+        @compat new(Float64(α), Float64(β))
+    end
+
+    BetaPrime(α::Float64) = BetaPrime(α)
+    BetaPrime() = new(1.0, 1.0)
 end
 
 @distr_support BetaPrime 0.0 Inf
 
-show(io::IO, d::BetaPrime) = ((α, β) = params(d); show_oneline(io, d, [(:α, α), (:β, β)]))
-
 #### Parameters
 
-params(d::BetaPrime) = params(d.betad)
+params(d::BetaPrime) = (d.α, d.β)
 
 
 #### Statistics
@@ -46,21 +49,20 @@ end
 
 pdf(d::BetaPrime, x::Float64) = exp(logpdf(d, x))
 
-cdf(d::BetaPrime, x::Float64) = cdf(d.betad, x / (1.0 + x))
-ccdf(d::BetaPrime, x::Float64) = ccdf(d.betad, x / (1.0 + x))
-logcdf(d::BetaPrime, x::Float64) = logcdf(d.betad, x / (1.0 + x))
-logccdf(d::BetaPrime, x::Float64) = logccdf(d.betad, x / (1.0 + x))
+cdf(d::BetaPrime, x::Float64) = betacdf(d.α, d.β, x / (1.0 + x))
+ccdf(d::BetaPrime, x::Float64) = betaccdf(d.α, d.β, x / (1.0 + x))
+logcdf(d::BetaPrime, x::Float64) = betalogcdf(d.α, d.β, x / (1.0 + x))
+logccdf(d::BetaPrime, x::Float64) = betalogccdf(d.α, d.β, x / (1.0 + x))
+
+quantile(d::BetaPrime, p::Float64) = (x = betainvcdf(d.α, d.β, p); x / (1.0 - x))
+cquantile(d::BetaPrime, p::Float64) = (x = betainvccdf(d.α, d.β, p); x / (1.0 - x))
+invlogcdf(d::BetaPrime, p::Float64) = (x = betainvlogcdf(d.α, d.β, p); x / (1.0 - x))
+invlogccdf(d::BetaPrime, p::Float64) = (x = betainvlogccdf(d.α, d.β, p); x / (1.0 - x))
 
-quantile(d::BetaPrime, p::Float64) = (x = quantile(d.betad, p); x / (1.0 - x))
-cquantile(d::BetaPrime, p::Float64) = (x = cquantile(d.betad, p); x / (1.0 - x))
-invlogcdf(d::BetaPrime, p::Float64) = (x = invlogcdf(d.betad, p); x / (1.0 - x))
-invlogccdf(d::BetaPrime, p::Float64) = (x = invlogccdf(d.betad, p); x / (1.0 - x))
-
 
 #### Sampling
 
-function rand(d::BetaPrime) 
+function rand(d::BetaPrime)
     (α, β) = params(d)
     rand(Gamma(α)) / rand(Gamma(β))
 end
-

src/univariate/continuous/chi.jl

@@ -1,69 +1,71 @@
 immutable Chi <: ContinuousUnivariateDistribution
-    chisqd::Chisq
+    ν::Float64
 
-    Chi(df::Real) = new(Chisq(df))
+    function Chi(ν::Real)
+        ν > zero(ν) || throw(ArgumentError("Chi: ν must be positive."))
+        @compat new(Float64(ν))
+    end
 end
 
 @distr_support Chi 0.0 Inf
 
 
 #### Parameters
 
-dof(d::Chi) = dof(d.chisqd)
-params(d::Chi) = (dof(d),)
+dof(d::Chi) = d.ν
+params(d::Chi) = (d.ν,)
 
 
 #### Statistics
 
-mean(d::Chi) = (k = dof(d); sqrt2 * gamma((k + 1.0) / 2.0) / gamma(k / 2.0))
+mean(d::Chi) = (h = d.ν * 0.5; sqrt2 * gamma(h + 0.5) / gamma(h))
 
-var(d::Chi) = dof(d) - mean(d)^2
+var(d::Chi) = d.ν - mean(d)^2
+_chi_skewness(μ::Float64, σ::Float64) = (σ2 = σ^2; σ3 = σ2 * σ; (μ / σ3) * (1.0 - 2.0 * σ2))
 
 function skewness(d::Chi)
-    μ, σ = mean(d), std(d)
-    (μ / σ^3) * (1.0 - 2.0 * σ^2)
+    μ = mean(d)
+    σ = sqrt(d.ν - μ^2)
+    _chi_skewness(μ, σ)
 end
 
 function kurtosis(d::Chi)
-    μ, σ, γ = mean(d), std(d), skewness(d)
+    μ = mean(d)
+    σ = sqrt(d.ν - μ^2)
+    γ = _chi_skewness(μ, σ)
     (2.0 / σ^2) * (1 - μ * σ * γ - σ^2)
 end
 
-function entropy(d::Chi)
-    k = dof(d)
-    lgamma(k / 2.0) - log(sqrt(2.0)) -
-        ((k - 1.0) / 2.0) * digamma(k / 2.0) + k / 2.0
-end
+entropy(d::Chi) = (ν = d.ν;
+    lgamma(ν / 2.0) - 0.5 * logtwo - ((ν - 1.0) / 2.0) * digamma(ν / 2.0) + ν / 2.0)
 
 function mode(d::Chi)
-    k = dof(d)
-    k >= 1.0 || error("Chi distribution has no mode when df < 1")
-    sqrt(k - 1.0)
+    d.ν >= 1.0 || error("Chi distribution has no mode when ν < 1")
+    sqrt(d.ν - 1.0)
 end
 
 
 #### Evaluation
 
 pdf(d::Chi, x::Float64) = exp(logpdf(d, x))
 
-function logpdf(d::Chi, x::Float64) 
-    k = dof(d)
-    (1.0 - 0.5 * k) * logtwo + (k - 1.0) * log(x) - 0.5 * x^2 - lgamma(0.5 * k)
-end
+logpdf(d::Chi, x::Float64) = (ν = d.ν;
+    (1.0 - 0.5 * ν) * logtwo + (ν - 1.0) * log(x) - 0.5 * x^2 - lgamma(0.5 * ν)
+)
 
-gradlogpdf(d::Chi, x::Float64) = x >= 0.0 ? (dof(d) - 1.0) / x - x : 0.0
+gradlogpdf(d::Chi, x::Float64) = x >= 0.0 ? (d.ν - 1.0) / x - x : 0.0
 
-cdf(d::Chi, x::Float64) = cdf(d.chisqd, x^2)
-ccdf(d::Chi, x::Float64) = ccdf(d.chisqd, x^2)
-logcdf(d::Chi, x::Float64) = logcdf(d.chisqd, x^2)
-logccdf(d::Chi, x::Float64) = logccdf(d.chisqd, x^2)
+cdf(d::Chi, x::Float64) = chisqcdf(d.ν, x^2)
+ccdf(d::Chi, x::Float64) = chisqccdf(d.ν, x^2)
+logcdf(d::Chi, x::Float64) = chisqlogcdf(d.ν, x^2)
+logccdf(d::Chi, x::Float64) = chisqlogccdf(d.ν, x^2)
 
-quantile(d::Chi, p::Float64) = sqrt(quantile(d.chisqd, p))
-cquantile(d::Chi, p::Float64) = sqrt(cquantile(d.chisqd, p))
-invlogcdf(d::Chi, p::Float64) = sqrt(invlogcdf(d.chisqd, p))
-invlogccdf(d::Chi, p::Float64) = sqrt(invlogccdf(d.chisqd, p))
+quantile(d::Chi, p::Float64) = sqrt(chisqinvcdf(d.ν, p))
+cquantile(d::Chi, p::Float64) = sqrt(chisqinvccdf(d.ν, p))
+invlogcdf(d::Chi, p::Float64) = sqrt(chisqinvlogcdf(d.ν, p))
+invlogccdf(d::Chi, p::Float64) = sqrt(chisqinvlogccdf(d.ν, p))
 
 
 #### Sampling
 
-rand(d::Chi) = sqrt(rand(d.chisqd))
+rand(d::Chi) = sqrt(_chisq_rand(d.ν))

src/univariate/continuous/chisq.jl

@@ -55,4 +55,5 @@ gradlogpdf(d::Chisq, x::Float64) =  x >= 0.0 ? (dof(d) / 2.0 - 1) / x - 0.5 : 0.
 
 #### Sampling
 
-rand(d::Chisq) = StatsFuns.Rmath.chisqrand(d.df)
+_chisq_rand(ν::Float64) = StatsFuns.Rmath.chisqrand(ν)
+rand(d::Chisq) = _chisq_rand(d.df)

test/dirichlet.jl

@@ -3,6 +3,8 @@
 using Distributions
 using Base.Test
 
+srand(34567)
+
 d = Dirichlet(3, 2.0)
 
 @test length(d) == 3
@@ -71,6 +73,5 @@ x = x ./ sum(x, 1)
 r = fit_mle(Dirichlet, x)
 @test_approx_eq_eps r.alpha d.alpha 0.25
 
-r = fit_mle(Dirichlet, x, fill(2.0, n))
-@test_approx_eq_eps r.alpha d.alpha 0.25
-
+# r = fit_mle(Dirichlet, x, fill(2.0, n))
+# @test_approx_eq_eps r.alpha d.alpha 0.25