diff --git a/src/edgeworth.jl b/src/edgeworth.jl index faa921c1c..944d2ed81 100644 --- a/src/edgeworth.jl +++ b/src/edgeworth.jl @@ -12,10 +12,10 @@ kurtosis(d::EdgeworthAbstract) = kurtosis(d.dist) / d.n immutable EdgeworthZ{D<:UnivariateDistribution} <: EdgeworthAbstract dist::D n::Float64 + function EdgeworthZ{T<:UnivariateDistribution}(d::T, n::Real) - n > zero(n) || - error("n must be positive") - @compat new(d, Float64(n)) + @check_args(EdgeworthZ, n > zero(n)) + new(d, n) end end EdgeworthZ(d::UnivariateDistribution,n::Real) = EdgeworthZ{typeof(d)}(d,n) @@ -78,9 +78,8 @@ immutable EdgeworthSum{D<:UnivariateDistribution} <: EdgeworthAbstract dist::D n::Float64 function EdgeworthSum{T<:UnivariateDistribution}(d::T, n::Real) - n > zero(n) || - error("n must be positive") - @compat new(d, Float64(n)) + @check_args(EdgeworthSum, n > zero(n)) + new(d, n) end end EdgeworthSum(d::UnivariateDistribution, n::Real) = EdgeworthSum{typeof(d)}(d,n) diff --git a/src/matrix/inversewishart.jl b/src/matrix/inversewishart.jl index d3625d950..4e9cee2ec 100644 --- a/src/matrix/inversewishart.jl +++ b/src/matrix/inversewishart.jl @@ -69,7 +69,7 @@ function _logpdf(d::InverseWishart, X::DenseMatrix{Float64}) -0.5 * ((df + p + 1) * logdet(Xcf) + trace(Xcf \ Ψ)) - d.c0 end -@compat _logpdf{T<:Real}(d::InverseWishart, X::DenseMatrix{T}) = _logpdf(d, Float64(X)) +_logpdf{T<:Real}(d::InverseWishart, X::DenseMatrix{T}) = _logpdf(d, convert(Matrix{Float64}, X)) #### Sampling diff --git a/src/matrix/wishart.jl b/src/matrix/wishart.jl index 73b5cd4ea..46fa1f668 100644 --- a/src/matrix/wishart.jl +++ b/src/matrix/wishart.jl @@ -55,7 +55,7 @@ function mode(d::Wishart) end end -function meanlogdet(d::Wishart) +function meanlogdet(d::Wishart) p = dim(d) df = d.df v = logdet(d.S) + p * logtwo @@ -81,7 +81,7 @@ function _logpdf(d::Wishart, X::DenseMatrix{Float64}) 0.5 * ((df - (p + 1)) * logdet(Xcf) - trace(d.S \ X)) - d.c0 end -@compat _logpdf{T<:Real}(d::Wishart, X::DenseMatrix{T}) = _logpdf(d, Float64(X)) +_logpdf{T<:Real}(d::Wishart, X::DenseMatrix{T}) = _logpdf(d, convert(Matrix{Float64}, X)) #### Sampling diff --git a/src/multivariate/dirichlet.jl b/src/multivariate/dirichlet.jl index f810f68c4..8824ddac0 100644 --- a/src/multivariate/dirichlet.jl +++ b/src/multivariate/dirichlet.jl @@ -8,7 +8,7 @@ immutable Dirichlet <: ContinuousMultivariateDistribution lmnB::Float64 = 0.0 for i in 1:length(alpha) ai = alpha[i] - ai > 0 || throw(ArgumentError("alpha must be a positive vector.")) + ai > 0 || throw(ArgumentError("Dirichlet: alpha must be a positive vector.")) alpha0 += ai lmnB += lgamma(ai) end diff --git a/src/testutils.jl b/src/testutils.jl index 4b6888fb3..d394dae93 100644 --- a/src/testutils.jl +++ b/src/testutils.jl @@ -435,7 +435,7 @@ end function test_stats(d::DiscreteUnivariateDistribution, vs::AbstractVector) # using definition (or an approximation) - @compat vf = Float64[Float64(v) for v in vs] + vf = Float64[v for v in vs] p = pdf(d, vf) xmean = dot(p, vf) xvar = dot(p, abs2(vf .- xmean)) diff --git a/src/univariate/continuous/arcsine.jl b/src/univariate/continuous/arcsine.jl index a738b627f..212a67f87 100644 --- a/src/univariate/continuous/arcsine.jl +++ b/src/univariate/continuous/arcsine.jl @@ -2,13 +2,8 @@ immutable Arcsine <: ContinuousUnivariateDistribution a::Float64 b::Float64 - function Arcsine(a::Float64, b::Float64) - a < b || throw(ArgumentError("a must be less than b.")) - new(a, b) - end - - @compat Arcsine(a::Real, b::Real) = Arcsine(Float64(a), Float64(b)) - @compat Arcsine(b::Real) = Arcsine(0.0, Float64(b)) + Arcsine(a::Real, b::Real) = (@check_args(Arcsine, a < b); new(a, b)) + Arcsine(b::Real) = (@check_args(Arcsine, b > zero(b)); new(0.0, b)) Arcsine() = new(0.0, 1.0) end diff --git a/src/univariate/continuous/beta.jl b/src/univariate/continuous/beta.jl index 6cbce08b4..dfce68a81 100644 --- a/src/univariate/continuous/beta.jl +++ b/src/univariate/continuous/beta.jl @@ -3,11 +3,9 @@ immutable Beta <: ContinuousUnivariateDistribution β::Float64 function Beta(α::Real, β::Real) - (α > zero(α) && β > zero(β)) || - throw(ArgumentError("α and β must be positive")) - @compat new(Float64(α), Float64(β)) + @check_args(Beta, α > zero(α) && β > zero(β)) + new(α, β) end - Beta(α::Real) = Beta(α, α) Beta() = new(1.0, 1.0) end diff --git a/src/univariate/continuous/betaprime.jl b/src/univariate/continuous/betaprime.jl index 69a493446..47e0f19ae 100644 --- a/src/univariate/continuous/betaprime.jl +++ b/src/univariate/continuous/betaprime.jl @@ -2,12 +2,11 @@ immutable BetaPrime <: ContinuousUnivariateDistribution α::Float64 β::Float64 - function BetaPrime(α::Float64, β::Float64) - (α > zero(α) && β > zero(β)) || error("α and β must be positive") - @compat new(Float64(α), Float64(β)) + function BetaPrime(α::Real, β::Real) + @check_args(BetaPrime, α > zero(α) && β > zero(β)) + new(α, β) end - - BetaPrime(α::Float64) = BetaPrime(α) + BetaPrime(α::Real) = BetaPrime(α) BetaPrime() = new(1.0, 1.0) end diff --git a/src/univariate/continuous/biweight.jl b/src/univariate/continuous/biweight.jl index 10ea035cd..a770e9d9d 100644 --- a/src/univariate/continuous/biweight.jl +++ b/src/univariate/continuous/biweight.jl @@ -1,53 +1,58 @@ immutable Biweight <: ContinuousUnivariateDistribution - location::Float64 - scale::Float64 - function Biweight(l::Real, s::Real) - s > zero(s) || error("scale must be positive") - @compat new(Float64(l), Float64(s)) - end -end + μ::Float64 + σ::Float64 -Biweight(location::Real) = Biweight(location, 1.0) -Biweight() = Biweight(0.0, 1.0) + Biweight(μ::Real, σ::Real) = (@check_args(Biweight, σ > zero(σ)); new(μ, σ)) + Biweight(μ::Real) = new(μ, 1.0) + Biweight() = new(0.0, 1.0) +end -@distr_support Biweight d.location-d.scale d.location+d.scale +@distr_support Biweight d.μ - d.σ d.μ + d.σ ## Parameters -params(d::Biweight) = (d.location, d.scale) +params(d::Biweight) = (d.μ, d.σ) ## Properties -mean(d::Biweight) = d.location -median(d::Biweight) = d.location -mode(d::Biweight) = d.location +mean(d::Biweight) = d.μ +median(d::Biweight) = d.μ +mode(d::Biweight) = d.μ -var(d::Biweight) = d.scale*d.scale/7.0 +var(d::Biweight) = d.σ^2 / 7.0 skewness(d::Biweight) = 0.0 -kurtosis(d::Biweight) = 1/21-3 +kurtosis(d::Biweight) = -2.9523809523809526 # = 1/21-3 ## Functions -function pdf(d::Biweight, x::Real) - u = abs(x - d.location)/d.scale - u >= 1 ? 0.0 : 0.9375*(1-u*u)^2/d.scale +function pdf(d::Biweight, x::Float64) + u = abs(x - d.μ) / d.σ + u >= 1.0 ? 0.0 : 0.9375 * (1 - u^2)^2 / d.σ end -function cdf(d::Biweight, x::Real) - u = (x - d.location)/d.scale - u <= -1 ? 0.0 : u >= 1 ? 1.0 : 0.0625*(1+u)^3*@horner(u,8.0,-9.0,3.0) +function cdf(d::Biweight, x::Float64) + u = (x - d.μ) / d.σ + u <= -1.0 ? 0.0 : + u >= 1.0 ? 1.0 : + 0.0625 * (u + 1.0)^3 * @horner(u,8.0,-9.0,3.0) end -function ccdf(d::Biweight, x::Real) - u = (d.location - x)/d.scale - u <= -1 ? 1.0 : u >= 1 ? 0.0 : 0.0625*(1+u)^3*@horner(u,8.0,-9.0,3.0) + +function ccdf(d::Biweight, x::Float64) + u = (d.μ - x) / d.σ + u <= -1.0 ? 1.0 : + u >= 1.0 ? 0.0 : + 0.0625 * (u + 1.0)^3 * @horner(u,8.0,-9.0,3.0) end @quantile_newton Biweight -function mgf(d::Biweight, t::Real) - a = d.scale*t - a2 = a*a - a == 0 ? one(a) : 15.0*exp(d.location*t)*(-3.0*cosh(a)+(a+3.0/a)*sinh(a))/(a2*a2) +function mgf(d::Biweight, t::Float64) + a = d.σ*t + a2 = a^2 + a == 0 ? 1.0 : + 15.0 * exp(d.μ * t) * (-3.0 * cosh(a) + (a + 3.0/a) * sinh(a)) / (a2^2) end -function cf(d::Biweight, t::Real) - a = d.scale*t - a2 = a*a - a == 0 ? complex(one(a)) : -15.0*cis(d.location*t)*(3.0*cos(a)+(a-3.0/a)*sin(a))/(a2*a2) + +function cf(d::Biweight, t::Float64) + a = d.σ * t + a2 = a^2 + a == 0 ? 1.0+0.0im : + -15.0 * cis(d.μ * t) * (3.0 * cos(a) + (a - 3.0/a) * sin(a)) / (a2^2) end diff --git a/src/univariate/continuous/cauchy.jl b/src/univariate/continuous/cauchy.jl index 2f89a0f37..129244476 100644 --- a/src/univariate/continuous/cauchy.jl +++ b/src/univariate/continuous/cauchy.jl @@ -1,13 +1,9 @@ immutable Cauchy <: ContinuousUnivariateDistribution μ::Float64 - β::Float64 + σ::Float64 - function Cauchy(μ::Real, β::Real) - β > zero(β) || error("Cauchy: scale must be positive") - @compat new(Float64(μ), Float64(β)) - end - - @compat Cauchy(μ::Real) = new(Float64(μ), 1.0) + Cauchy(μ::Real, σ::Real) = (@check_args(Cauchy, σ > zero(σ)); new(μ, σ)) + Cauchy(μ::Real) = new(μ, 1.0) Cauchy() = new(0.0, 1.0) end @@ -16,54 +12,54 @@ end #### Parameters location(d::Cauchy) = d.μ -scale(d::Cauchy) = d.β +scale(d::Cauchy) = d.σ -params(d::Cauchy) = (d.μ, d.β) +params(d::Cauchy) = (d.μ, d.σ) #### Statistics mean(d::Cauchy) = NaN -median(d::Cauchy) = location(d) -mode(d::Cauchy) = location(d) +median(d::Cauchy) = d.μ +mode(d::Cauchy) = d.μ var(d::Cauchy) = NaN skewness(d::Cauchy) = NaN kurtosis(d::Cauchy) = NaN -entropy(d::Cauchy) = log(scale(d)) + log4π +entropy(d::Cauchy) = log4π + log(d.σ) #### Functions -zval(d::Cauchy, x::Float64) = (x - d.μ) / d.β -xval(d::Cauchy, z::Float64) = d.μ + z * d.β +zval(d::Cauchy, x::Float64) = (x - d.μ) / d.σ +xval(d::Cauchy, z::Float64) = d.μ + z * d.σ -pdf(d::Cauchy, x::Float64) = 1.0 / (π * scale(d) * (1 + zval(d, x)^2)) -logpdf(d::Cauchy, x::Float64) = - (logπ + log(scale(d)) + log1psq(zval(d, x))) +pdf(d::Cauchy, x::Float64) = 1.0 / (π * scale(d) * (1.0 + zval(d, x)^2)) +logpdf(d::Cauchy, x::Float64) = - (log1psq(zval(d, x)) + logπ + log(d.σ)) function cdf(d::Cauchy, x::Float64) - μ, β = params(d) - invπ * atan2(x - μ, β) + 0.5 + μ, σ = params(d) + invπ * atan2(x - μ, σ) + 0.5 end function ccdf(d::Cauchy, x::Float64) - μ, β = params(d) - invπ * atan2(μ - x, β) + 0.5 + μ, σ = params(d) + invπ * atan2(μ - x, σ) + 0.5 end function quantile(d::Cauchy, p::Float64) - μ, β = params(d) - μ + β * tan(π * (p - 0.5)) + μ, σ = params(d) + μ + σ * tan(π * (p - 0.5)) end function cquantile(d::Cauchy, p::Float64) - μ, β = params(d) - μ + β * tan(π * (0.5 - p)) + μ, σ = params(d) + μ + σ * tan(π * (0.5 - p)) end mgf(d::Cauchy, t::Real) = t == zero(t) ? 1.0 : NaN -cf(d::Cauchy, t::Real) = exp(im * (t * d.μ) - d.β * abs(t)) +cf(d::Cauchy, t::Real) = exp(im * (t * d.μ) - d.σ * abs(t)) #### Fitting diff --git a/src/univariate/continuous/chi.jl b/src/univariate/continuous/chi.jl index 9b98e8aa6..cdd3bd4bd 100644 --- a/src/univariate/continuous/chi.jl +++ b/src/univariate/continuous/chi.jl @@ -1,15 +1,11 @@ immutable Chi <: ContinuousUnivariateDistribution ν::Float64 - function Chi(ν::Real) - ν > zero(ν) || throw(ArgumentError("Chi: ν must be positive.")) - @compat new(Float64(ν)) - end + Chi(ν::Real) = (@check_args(Chi, ν > zero(ν)); new(ν)) end @distr_support Chi 0.0 Inf - #### Parameters dof(d::Chi) = d.ν diff --git a/src/univariate/continuous/chisq.jl b/src/univariate/continuous/chisq.jl index 60bb2ffd7..d56eff73b 100644 --- a/src/univariate/continuous/chisq.jl +++ b/src/univariate/continuous/chisq.jl @@ -1,59 +1,55 @@ immutable Chisq <: ContinuousUnivariateDistribution - df::Float64 + ν::Float64 - function Chisq(k::Real) - k > zero(k) || error("The degree of freedom k must be positive") - @compat new(Float64(k)) - end + Chisq(ν::Real) = (@check_args(Chisq, ν > zero(ν)); new(ν)) end @distr_support Chisq 0.0 Inf #### Parameters -dof(d::Chisq) = d.df -params(d::Chisq) = (d.df,) +dof(d::Chisq) = d.ν +params(d::Chisq) = (d.ν,) #### Statistics -mean(d::Chisq) = dof(d) +mean(d::Chisq) = d.ν -var(d::Chisq) = 2.0 * dof(d) +var(d::Chisq) = 2.0 * d.ν -skewness(d::Chisq) = sqrt(8.0 / dof(d)) +skewness(d::Chisq) = sqrt(8.0 / d.ν) -kurtosis(d::Chisq) = 12.0 / dof(d) +kurtosis(d::Chisq) = 12.0 / d.ν -mode(d::Chisq) = d.df > 2.0 ? d.df - 2.0 : 0.0 +mode(d::Chisq) = d.ν > 2.0 ? d.ν - 2.0 : 0.0 function median(d::Chisq; approx::Bool=false) if approx - k = dof(d) - return k * (1.0 - 2.0 / (9.0 * k))^3 + return d.ν * (1.0 - 2.0 / (9.0 * d.ν))^3 else return quantile(d, 0.5) end end function entropy(d::Chisq) - hk = 0.5 * dof(d) - hk + logtwo + lgamma(hk) + (1.0 - hk) * digamma(hk) + hν = 0.5 * d.ν + hν + logtwo + lgamma(hν) + (1.0 - hν) * digamma(hν) end #### Evaluation -@_delegate_statsfuns Chisq chisq df +@_delegate_statsfuns Chisq chisq ν -mgf(d::Chisq, t::Real) = (1.0 - 2.0 * t)^(-dof(d) / 2.0) +mgf(d::Chisq, t::Real) = (1.0 - 2.0 * t)^(-d.ν * 0.5) -cf(d::Chisq, t::Real) = (1.0 - 2.0 * im * t)^(-dof(d) / 2.0) +cf(d::Chisq, t::Real) = (1.0 - 2.0 * im * t)^(-d.ν * 0.5) -gradlogpdf(d::Chisq, x::Float64) = x >= 0.0 ? (dof(d) / 2.0 - 1) / x - 0.5 : 0.0 +gradlogpdf(d::Chisq, x::Float64) = x > 0.0 ? (d.ν * 0.5 - 1) / x - 0.5 : 0.0 #### Sampling _chisq_rand(ν::Float64) = StatsFuns.Rmath.chisqrand(ν) -rand(d::Chisq) = _chisq_rand(d.df) +rand(d::Chisq) = _chisq_rand(d.ν) diff --git a/src/univariate/continuous/cosine.jl b/src/univariate/continuous/cosine.jl index 2142b1e73..9f3ee9364 100644 --- a/src/univariate/continuous/cosine.jl +++ b/src/univariate/continuous/cosine.jl @@ -5,26 +5,22 @@ immutable Cosine <: ContinuousUnivariateDistribution μ::Float64 - s::Float64 + σ::Float64 - function Cosine(μ::Real, s::Real) - s > 0.0 || error("s must be positive.") - @compat new(Float64(μ), Float64(s)) - end - - @compat Cosine(μ::Real) = new(Float64(μ), 1.0) + Cosine(μ::Real, σ::Real) = (@check_args(Cosine, σ > zero(σ)); new(μ, σ)) + Cosine(μ::Real) = new(μ, 1.0) Cosine() = new(0.0, 1.0) end -@distr_support Cosine d.μ - d.s d.μ + d.s +@distr_support Cosine d.μ - d.σ d.μ + d.σ #### Parameters location(d::Cosine) = d.μ -scale(d::Cosine) = d.s +scale(d::Cosine) = d.σ -params(d::Cosine) = (d.μ, d.s) +params(d::Cosine) = (d.μ, d.σ) #### Statistics @@ -35,8 +31,7 @@ median(d::Cosine) = d.μ mode(d::Cosine) = d.μ -const _cosined_varcoef = 0.13069096604865779 # 1 / 3 - 2 / π^2 -var(d::Cosine) = d.s^2 * _cosined_varcoef +var(d::Cosine) = d.σ^2 * 0.13069096604865779 # 0.130... = 1/3 - 2 / π^2 skewness(d::Cosine) = 0.0 @@ -47,9 +42,8 @@ kurtosis(d::Cosine) = -0.59376287559828102362 function pdf(d::Cosine, x::Float64) if insupport(d, x) - μ, s = params(d) - z = (x - μ) / s - return (1.0 + cospi(z)) / (2 * s) + z = (x - d.μ) / d.σ + return (1.0 + cospi(z)) / (2 * d.σ) else return 0.0 end @@ -58,14 +52,12 @@ end logpdf(d::Cosine, x::Float64) = insupport(d, x) ? log(pdf(d, x)) : -Inf function cdf(d::Cosine, x::Float64) - μ, s = params(d) - z = (x - μ) / s + z = (x - d.μ) / d.σ 0.5 * (1.0 + z + sinpi(z) * invπ) end function ccdf(d::Cosine, x::Float64) - μ, s = params(d) - nz = (μ - x) / s + nz = (d.μ - x) / d.σ 0.5 * (1.0 + nz + sinpi(nz) * invπ) end diff --git a/src/univariate/continuous/epanechnikov.jl b/src/univariate/continuous/epanechnikov.jl index c56302563..4e9b066bd 100644 --- a/src/univariate/continuous/epanechnikov.jl +++ b/src/univariate/continuous/epanechnikov.jl @@ -1,51 +1,59 @@ immutable Epanechnikov <: ContinuousUnivariateDistribution - location::Float64 - scale::Float64 - function Epanechnikov(l::Real, s::Real) - s > zero(s) || error("scale must be positive") - @compat new(Float64(l), Float64(s)) - end -end + μ::Float64 + σ::Float64 -Epanechnikov(location::Real) = Epanechnikov(location, 1.0) -Epanechnikov() = Epanechnikov(0.0, 1.0) + Epanechnikov(μ::Real, σ::Real) = (@check_args(Epanechnikov, σ > zero(σ)); new(μ, σ)) + Epanechnikov(μ::Real) = new(μ, 1.0) + Epanechnikov() = new(0.0, 1.0) +end -@distr_support Epanechnikov d.location-d.scale d.location+d.scale +@distr_support Epanechnikov d.μ - d.σ d.μ + d.σ ## Parameters -params(d::Epanechnikov) = (d.location, d.scale) + +location(d::Epanechnikov) = d.μ +scale(d::Epanechnikov) = d.σ +params(d::Epanechnikov) = (d.μ, d.σ) ## Properties -mean(d::Epanechnikov) = d.location -median(d::Epanechnikov) = d.location -mode(d::Epanechnikov) = d.location +mean(d::Epanechnikov) = d.μ +median(d::Epanechnikov) = d.μ +mode(d::Epanechnikov) = d.μ -var(d::Epanechnikov) = d.scale*d.scale/5 +var(d::Epanechnikov) = d.σ^2 / 5 skewness(d::Epanechnikov) = 0.0 -kurtosis(d::Epanechnikov) = 3/35-3 +kurtosis(d::Epanechnikov) = -2.914285714285714 # 3/35-3 ## Functions -function pdf(d::Epanechnikov, x::Real) - u = abs(x - d.location)/d.scale - u >= 1 ? 0.0 : 0.75*(1-u*u)/d.scale +function pdf(d::Epanechnikov, x::Float64) + u = abs(x - d.μ) / d.σ + u >= 1 ? 0.0 : 0.75 * (1 - u^2) / d.σ end -function cdf(d::Epanechnikov, x::Real) - u = (x - d.location)/d.scale - u <= -1 ? 0.0 : u >= 1 ? 1.0 : 0.5+u*(0.75-0.25*u*u) + +function cdf(d::Epanechnikov, x::Float64) + u = (x - d.μ) / d.σ + u <= -1 ? 0.0 : + u >= 1 ? 1.0 : + 0.5 + u * (0.75 - 0.25 * u^2) end -function ccdf(d::Epanechnikov, x::Real) - u = (d.location - x)/d.scale - u <= -1 ? 1.0 : u >= 1 ? 0.0 : 0.5+u*(0.75-0.25*u*u) + +function ccdf(d::Epanechnikov, x::Float64) + u = (d.μ - x) / d.σ + u <= -1 ? 1.0 : + u >= 1 ? 0.0 : + 0.5 + u * (0.75 - 0.25 * u^2) end @quantile_newton Epanechnikov -function mgf(d::Epanechnikov, t::Real) - a = d.scale*t - a == 0 ? one(a) : 3.0*exp(d.location*t)*(cosh(a)-sinh(a)/a)/(a*a) +function mgf(d::Epanechnikov, t::Float64) + a = d.σ * t + a == 0 ? 1.0 : + 3.0 * exp(d.μ * t) * (cosh(a) - sinh(a) / a) / a^2 end -function cf(d::Epanechnikov, t::Real) - a = d.scale*t - a == 0 ? complex(one(a)) : -3.0*exp(im*d.location*t)*(cos(a)-sin(a)/a)/(a*a) +function cf(d::Epanechnikov, t::Float64) + a = d.σ * t + a == 0 ? 1.0+0.0im : + -3.0 * exp(im * d.μ * t) * (cos(a) - sin(a) / a) / a^2 end diff --git a/src/univariate/continuous/erlang.jl b/src/univariate/continuous/erlang.jl index df7fe38f9..90ef7575d 100644 --- a/src/univariate/continuous/erlang.jl +++ b/src/univariate/continuous/erlang.jl @@ -2,10 +2,9 @@ immutable Erlang <: ContinuousUnivariateDistribution gammad::Gamma function Erlang(α::Real, θ::Real) - isinteger(α) || error("Erlang shape parameter must be an integer") + @check_args(Erlang, isinteger(α) && α >= zero(α)) new(Gamma(α, θ)) end - Erlang(α::Real) = Erlang(α, 1.0) Erlang() = new(Gamma()) end @@ -56,4 +55,3 @@ cf(d::Erlang, t::Real) = cf(d.gammad, t) #### Sampling rand(d::Erlang) = rand(d.gammad) - diff --git a/src/univariate/continuous/exponential.jl b/src/univariate/continuous/exponential.jl index 31ac92f27..d112a9b43 100644 --- a/src/univariate/continuous/exponential.jl +++ b/src/univariate/continuous/exponential.jl @@ -1,11 +1,7 @@ immutable Exponential <: ContinuousUnivariateDistribution - β::Float64 # note: scale not rate - - function Exponential(β::Real) - β > zero(β) || error("scale must be positive") - @compat new(Float64(β)) - end + θ::Float64 # note: scale not rate + Exponential(θ::Real) = (@check_args(Exponential, θ > zero(θ)); new(θ)) Exponential() = new(1.0) end @@ -14,33 +10,33 @@ end #### Parameters -scale(d::Exponential) = d.β -rate(d::Exponential) = 1.0 / d.β +scale(d::Exponential) = d.θ +rate(d::Exponential) = 1.0 / d.θ -params(d::Exponential) = (d.β,) +params(d::Exponential) = (d.θ,) #### Statistics -mean(d::Exponential) = scale(d) +mean(d::Exponential) = d.θ -median(d::Exponential) = logtwo * scale(d) +median(d::Exponential) = logtwo * d.θ mode(d::Exponential) = 0.0 -var(d::Exponential) = scale(d)^2 +var(d::Exponential) = d.θ^2 skewness(d::Exponential) = 2.0 kurtosis(d::Exponential) = 6.0 -entropy(d::Exponential) = 1.0 + log(scale(d)) +entropy(d::Exponential) = 1.0 + log(d.θ) #### Evaluation -zval(d::Exponential, x::Float64) = x / d.β -xval(d::Exponential, z::Float64) = z * d.β +zval(d::Exponential, x::Float64) = x / d.θ +xval(d::Exponential, z::Float64) = z * d.θ pdf(d::Exponential, x::Float64) = (λ = rate(d); x < 0.0 ? 0.0 : λ * exp(-λ * x)) logpdf(d::Exponential, x::Float64) = (λ = rate(d); x < 0.0 ? -Inf : log(λ) - λ * x) @@ -72,7 +68,7 @@ immutable ExponentialStats <: SufficientStats sx::Float64 # (weighted) sum of x sw::Float64 # sum of sample weights - @compat ExponentialStats(sx::Real, sw::Real) = new(Float64(sx), Float64(sw)) + ExponentialStats(sx::Real, sw::Real) = new(sx, sw) end suffstats{T<:Real}(::Type{Exponential}, x::AbstractArray{T}) = ExponentialStats(sum(x), length(x)) diff --git a/src/univariate/continuous/fdist.jl b/src/univariate/continuous/fdist.jl index d56df8de4..79630108e 100644 --- a/src/univariate/continuous/fdist.jl +++ b/src/univariate/continuous/fdist.jl @@ -1,45 +1,46 @@ immutable FDist <: ContinuousUnivariateDistribution - d1::Float64 - d2::Float64 + ν1::Float64 + ν2::Float64 - function FDist(d1::Real, d2::Real) - d1 > zero(d1) && d2 > zero(d2) || error("Degrees of freedom must be positive") - @compat new(Float64(d1), Float64(d2)) + function FDist(ν1::Real, ν2::Real) + @check_args(FDist, ν1 > zero(ν1) && ν2 > zero(ν2)) + new(ν1, ν2) end end @distr_support FDist 0.0 Inf + #### Parameters -params(d::FDist) = (d.d1, d.d2) +params(d::FDist) = (d.ν1, d.ν2) #### Statistics -mean(d::FDist) = (d2 = d.d2; d2 > 2.0 ? d2 / (d2 - 2.0) : NaN) +mean(d::FDist) = (ν2 = d.ν2; ν2 > 2.0 ? ν2 / (ν2 - 2.0) : NaN) -mode(d::FDist) = ((d1, d2) = params(d); d1 > 2.0 ? ((d1 - 2.0)/d1) * (d2 / (d2 + 2.0)) : 0.0) +mode(d::FDist) = ((ν1, ν2) = params(d); ν1 > 2.0 ? ((ν1 - 2.0)/ν1) * (ν2 / (ν2 + 2.0)) : 0.0) function var(d::FDist) - (d1, d2) = params(d) - d2 > 4.0 ? 2.0 * d2^2 * (d1 + d2 - 2.0) / (d1 * (d2 - 2.0)^2 * (d2 - 4.0)) : NaN + (ν1, ν2) = params(d) + ν2 > 4.0 ? 2.0 * ν2^2 * (ν1 + ν2 - 2.0) / (ν1 * (ν2 - 2.0)^2 * (ν2 - 4.0)) : NaN end function skewness(d::FDist) - (d1, d2) = params(d) - if d2 > 6.0 - return (2.0 * d1 + d2 - 2.0) * sqrt(8.0 * (d2 - 4.0)) / ((d2 - 6.0) * sqrt(d1 * (d1 + d2 - 2.0))) + (ν1, ν2) = params(d) + if ν2 > 6.0 + return (2.0 * ν1 + ν2 - 2.0) * sqrt(8.0 * (ν2 - 4.0)) / ((ν2 - 6.0) * sqrt(ν1 * (ν1 + ν2 - 2.0))) else return NaN end end function kurtosis(d::FDist) - (d1, d2) = params(d) - if d2 > 8.0 - a = d1 * (5. * d2 - 22.) * (d1 + d2 - 2.) + (d2 - 4.) * (d2 - 2.)^2 - b = d1 * (d2 - 6.) * (d2 - 8.) * (d2 - 2.) + (ν1, ν2) = params(d) + if ν2 > 8.0 + a = ν1 * (5. * ν2 - 22.) * (ν1 + ν2 - 2.) + (ν2 - 4.) * (ν2 - 2.)^2 + b = ν1 * (ν2 - 6.) * (ν2 - 8.) * (ν2 - 2.) return 12. * a / b else return NaN @@ -47,17 +48,17 @@ function kurtosis(d::FDist) end function entropy(d::FDist) - (d1, d2) = params(d) - hd1 = d1 * 0.5 - hd2 = d2 * 0.5 - hs = (d1 + d2) * 0.5 - return log(d2 / d1) + lgamma(hd1) + lgamma(hd2) - lgamma(hs) + - (1.0 - hd1) * digamma(hd1) + (-1.0 - hd2) * digamma(hd2) + + (ν1, ν2) = params(d) + hν1 = ν1 * 0.5 + hν2 = ν2 * 0.5 + hs = (ν1 + ν2) * 0.5 + return log(ν2 / ν1) + lgamma(hν1) + lgamma(hν2) - lgamma(hs) + + (1.0 - hν1) * digamma(hν1) + (-1.0 - hν2) * digamma(hν2) + hs * digamma(hs) end #### Evaluation & Sampling -@_delegate_statsfuns FDist fdist d1 d2 +@_delegate_statsfuns FDist fdist ν1 ν2 -rand(d::FDist) = StatsFuns.Rmath.fdistrand(d.d1, d.d2) +rand(d::FDist) = StatsFuns.Rmath.fdistrand(d.ν1, d.ν2) diff --git a/src/univariate/continuous/frechet.jl b/src/univariate/continuous/frechet.jl index f0faf2e84..4336c0da5 100644 --- a/src/univariate/continuous/frechet.jl +++ b/src/univariate/continuous/frechet.jl @@ -1,12 +1,11 @@ immutable Frechet <: ContinuousUnivariateDistribution α::Float64 - β::Float64 + θ::Float64 - function Frechet(α::Real, β::Real) - α > zero(α) && β > zero(β) || error("Both shape and scale must be positive") - @compat new(Float64(α), Float64(β)) + function Frechet(α::Real, θ::Real) + @check_args(Frechet, α > zero(α) && θ > zero(θ)) + new(α, θ) end - Frechet(α::Real) = Frechet(α, 1.0) Frechet() = new(1.0, 1.0) end @@ -17,23 +16,23 @@ end #### Parameters shape(d::Frechet) = d.α -scale(d::Frechet) = d.β -params(d::Frechet) = (d.α, d.β) +scale(d::Frechet) = d.θ +params(d::Frechet) = (d.α, d.θ) #### Statistics -mean(d::Frechet) = (α = d.α; α > 1.0 ? d.β * gamma(1.0 - 1.0 / α) : Inf) +mean(d::Frechet) = (α = d.α; α > 1.0 ? d.θ * gamma(1.0 - 1.0 / α) : Inf) -median(d::Frechet) = d.β * logtwo^(-1.0 / d.α) +median(d::Frechet) = d.θ * logtwo^(-1.0 / d.α) -mode(d::Frechet) = (iα = -1.0/d.α; d.β * (1.0 - iα) ^ iα) +mode(d::Frechet) = (iα = -1.0/d.α; d.θ * (1.0 - iα) ^ iα) function var(d::Frechet) if d.α > 2.0 - iα = 1.0 / d.α - return d.β^2 * (gamma(1.0 - 2.0 * iα) - gamma(1.0 - iα)^2) - else + iα = 1.0 / d.α + return d.θ^2 * (gamma(1.0 - 2.0 * iα) - gamma(1.0 - iα)^2) + else return Inf end end @@ -60,22 +59,22 @@ function kurtosis(d::Frechet) return (g4 - 4.0 * g3 * g1 + 3 * g2^2) / ((g2 - g1^2)^2) - 6.0 else return Inf - end + end end function entropy(d::Frechet) const γ = 0.57721566490153286060 # γ is the Euler-Mascheroni constant - 1.0 + γ / d.α + γ + log(d.β / d.α) + 1.0 + γ / d.α + γ + log(d.θ / d.α) end #### Evaluation function logpdf(d::Frechet, x::Float64) - (α, β) = params(d) + (α, θ) = params(d) if x > 0.0 - z = β / x - return log(α / β) + (1.0 + α) * log(z) - z^α + z = θ / x + return log(α / θ) + (1.0 + α) * log(z) - z^α else return -Inf end @@ -83,23 +82,21 @@ end pdf(d::Frechet, x::Float64) = exp(logpdf(d, x)) -cdf(d::Frechet, x::Float64) = x > 0.0 ? exp(-((d.β / x) ^ d.α)) : 0.0 -ccdf(d::Frechet, x::Float64) = x > 0.0 ? -expm1(-((d.β / x) ^ d.α)) : 1.0 -logcdf(d::Frechet, x::Float64) = x > 0.0 ? -(d.β / x) ^ d.α : -Inf -logccdf(d::Frechet, x::Float64) = x > 0.0 ? log1mexp(-((d.β / x) ^ d.α)) : 0.0 +cdf(d::Frechet, x::Float64) = x > 0.0 ? exp(-((d.θ / x) ^ d.α)) : 0.0 +ccdf(d::Frechet, x::Float64) = x > 0.0 ? -expm1(-((d.θ / x) ^ d.α)) : 1.0 +logcdf(d::Frechet, x::Float64) = x > 0.0 ? -(d.θ / x) ^ d.α : -Inf +logccdf(d::Frechet, x::Float64) = x > 0.0 ? log1mexp(-((d.θ / x) ^ d.α)) : 0.0 -quantile(d::Frechet, p::Float64) = d.β * (-log(p)) ^ (-1.0 / d.α) -cquantile(d::Frechet, p::Float64) = d.β * (-log1p(-p)) ^ (-1.0 / d.α) -invlogcdf(d::Frechet, lp::Float64) = d.β * (-lp)^(-1.0 / d.α) -invlogccdf(d::Frechet, lp::Float64) = d.β * (-log1mexp(lp))^(-1.0 / d.α) +quantile(d::Frechet, p::Float64) = d.θ * (-log(p)) ^ (-1.0 / d.α) +cquantile(d::Frechet, p::Float64) = d.θ * (-log1p(-p)) ^ (-1.0 / d.α) +invlogcdf(d::Frechet, lp::Float64) = d.θ * (-lp)^(-1.0 / d.α) +invlogccdf(d::Frechet, lp::Float64) = d.θ * (-log1mexp(lp))^(-1.0 / d.α) function gradlogpdf(d::Frechet, x::Float64) - (α, β) = params(d) - insupport(Frechet, x) ? -(α + 1.0) / x + α * (β^α) * x^(-α-1.0) : 0.0 + (α, θ) = params(d) + insupport(Frechet, x) ? -(α + 1.0) / x + α * (θ^α) * x^(-α-1.0) : 0.0 end ## Sampling -rand(d::Frechet) = d.β * randexp() ^ (-1.0 / d.α) - - +rand(d::Frechet) = d.θ * randexp() ^ (-1.0 / d.α) diff --git a/src/univariate/continuous/gamma.jl b/src/univariate/continuous/gamma.jl index 16d96088c..13cd0dda3 100644 --- a/src/univariate/continuous/gamma.jl +++ b/src/univariate/continuous/gamma.jl @@ -1,13 +1,15 @@ immutable Gamma <: ContinuousUnivariateDistribution α::Float64 - β::Float64 + θ::Float64 - function Gamma(α::Real, β::Real) - α > zero(α) && β > zero(β) || error("Gamma: both shape and scale must be positive") - @compat new(Float64(α), Float64(β)) + function Gamma(α::Real, θ::Real) + @check_args(Gamma, α > zero(α) && θ > zero(θ)) + new(α, θ) + end + function Gamma(α::Real) + @check_args(Gamma, α > zero(α)) + new(α, 1.0) end - - Gamma(α::Real) = Gamma(α, 1.0) Gamma() = new(1.0, 1.0) end @@ -17,45 +19,45 @@ end #### Parameters shape(d::Gamma) = d.α -scale(d::Gamma) = d.β -rate(d::Gamma) = 1.0 / d.β +scale(d::Gamma) = d.θ +rate(d::Gamma) = 1.0 / d.θ -params(d::Gamma) = (d.α, d.β) +params(d::Gamma) = (d.α, d.θ) #### Statistics -mean(d::Gamma) = d.α * d.β +mean(d::Gamma) = d.α * d.θ -var(d::Gamma) = d.α * d.β^2 +var(d::Gamma) = d.α * d.θ^2 skewness(d::Gamma) = 2.0 / sqrt(d.α) kurtosis(d::Gamma) = 6.0 / d.α function mode(d::Gamma) - (α, β) = params(d) - α >= 1.0 ? β * (α - 1.0) : error("Gamma has no mode when shape < 1.0") + (α, θ) = params(d) + α >= 1.0 ? θ * (α - 1.0) : error("Gamma has no mode when shape < 1.0") end function entropy(d::Gamma) - (α, β) = params(d) - α + lgamma(α) + (1.0 - α) * digamma(α) + log(β) + (α, θ) = params(d) + α + lgamma(α) + (1.0 - α) * digamma(α) + log(θ) end -mgf(d::Gamma, t::Real) = (1.0 - t * d.β)^(-d.α) +mgf(d::Gamma, t::Real) = (1.0 - t * d.θ)^(-d.α) -cf(d::Gamma, t::Real) = (1.0 - im * t * d.β)^(-d.α) +cf(d::Gamma, t::Real) = (1.0 - im * t * d.θ)^(-d.α) #### Evaluation & Sampling -@_delegate_statsfuns Gamma gamma α β +@_delegate_statsfuns Gamma gamma α θ gradlogpdf(d::Gamma, x::Float64) = - insupport(Gamma, x) ? (d.α - 1.0) / x - 1.0 / d.β : 0.0 + insupport(Gamma, x) ? (d.α - 1.0) / x - 1.0 / d.θ : 0.0 -rand(d::Gamma) = StatsFuns.Rmath.gammarand(d.α, d.β) +rand(d::Gamma) = StatsFuns.Rmath.gammarand(d.α, d.θ) #### Fit model @@ -65,7 +67,7 @@ immutable GammaStats <: SufficientStats slogx::Float64 # (weighted) sum of log(x) tw::Float64 # total sample weight - @compat GammaStats(sx::Real, slogx::Real, tw::Real) = new(Float64(sx), Float64(slogx), Float64(tw)) + GammaStats(sx::Real, slogx::Real, tw::Real) = new(sx, slogx, tw) end function suffstats{T<:Real}(::Type{Gamma}, x::AbstractArray{T}) diff --git a/src/univariate/continuous/gumbel.jl b/src/univariate/continuous/gumbel.jl index 55de73416..cc8ad5105 100644 --- a/src/univariate/continuous/gumbel.jl +++ b/src/univariate/continuous/gumbel.jl @@ -1,13 +1,9 @@ immutable Gumbel <: ContinuousUnivariateDistribution - μ::Float64 # location - β::Float64 # scale + μ::Float64 # location + θ::Float64 # scale - function Gumbel(μ::Real, β::Real) - β > zero(β) || error("The scale of Gumbel must be positive") - @compat new(Float64(μ), Float64(β)) - end - - Gumbel(μ::Real) = Gumbel(μ, 1.0) + Gumbel(μ::Real, θ::Real) = (@check_args(Gumbel, θ > zero(θ)); new(μ, θ)) + Gumbel(μ::Real) = new(μ, 1.0) Gumbel() = new(0.0, 1.0) end @@ -19,52 +15,50 @@ const DoubleExponential = Gumbel #### Parameters location(d::Gumbel) = d.μ -scale(d::Gumbel) = d.β -params(d::Gumbel) = (d.μ, d.β) +scale(d::Gumbel) = d.θ +params(d::Gumbel) = (d.μ, d.θ) #### Statistics -mean(d::Gumbel) = d.μ + d.β * 0.57721566490153286 +mean(d::Gumbel) = d.μ + d.θ * 0.57721566490153286 -median(d::Gumbel) = d.μ + d.β * 0.366512920581664327 +median(d::Gumbel) = d.μ + d.θ * 0.366512920581664327 mode(d::Gumbel) = d.μ -var(d::Gumbel) = 1.6449340668482264 * d.β^2 +var(d::Gumbel) = 1.6449340668482264 * d.θ^2 skewness(d::Gumbel) = 1.13954709940464866 kurtosis(d::Gumbel) = 2.4 -entropy(d::Gumbel) = 1.57721566490153286 + log(d.β) +entropy(d::Gumbel) = 1.57721566490153286 + log(d.θ) #### Evaluation -zval(d::Gumbel, x::Float64) = (x - d.μ) / d.β -xval(d::Gumbel, z::Float64) = x * d.β + d.μ +zval(d::Gumbel, x::Float64) = (x - d.μ) / d.θ +xval(d::Gumbel, z::Float64) = x * d.θ + d.μ function pdf(d::Gumbel, x::Float64) z = zval(d, x) - exp(-z - exp(-z)) / d.β + exp(-z - exp(-z)) / d.θ end function logpdf(d::Gumbel, x::Float64) z = zval(d, x) - - (z + exp(-z) + log(d.β)) + - (z + exp(-z) + log(d.θ)) end cdf(d::Gumbel, x::Float64) = exp(-exp(-zval(d, x))) logcdf(d::Gumbel, x::Float64) = -exp(-zval(d, x)) -quantile(d::Gumbel, p::Float64) = d.μ - d.β * log(-log(p)) +quantile(d::Gumbel, p::Float64) = d.μ - d.θ * log(-log(p)) -gradlogpdf(d::Gumbel, x::Float64) = - (1.0 + exp((d.μ - x) / d.β)) / d.β +gradlogpdf(d::Gumbel, x::Float64) = - (1.0 + exp((d.μ - x) / d.θ)) / d.θ #### Sampling rand(d::Gumbel) = quantile(d, rand()) - - diff --git a/src/univariate/continuous/inversegamma.jl b/src/univariate/continuous/inversegamma.jl index 4aa6b98ee..2d752876f 100644 --- a/src/univariate/continuous/inversegamma.jl +++ b/src/univariate/continuous/inversegamma.jl @@ -1,10 +1,10 @@ immutable InverseGamma <: ContinuousUnivariateDistribution invd::Gamma - β::Float64 + θ::Float64 - function InverseGamma(α::Real, β::Real) - (α > zero(α) && β > zero(β)) || error("Both shape and scale must be positive.") - @compat new(Gamma(α, 1.0 / β), Float64(β)) + function InverseGamma(α::Real, θ::Real) + @check_args(InverseGamma, α > zero(α) && θ > zero(θ)) + new(Gamma(α, 1.0 / θ), θ) end InverseGamma(α::Real) = InverseGamma(α, 1.0) @@ -17,7 +17,7 @@ end #### Parameters shape(d::InverseGamma) = shape(d.invd) -scale(d::InverseGamma) = d.β +scale(d::InverseGamma) = d.θ rate(d::InverseGamma) = scale(d.invd) params(d::InverseGamma) = (shape(d), scale(d)) @@ -25,13 +25,13 @@ params(d::InverseGamma) = (shape(d), scale(d)) #### Parameters -mean(d::InverseGamma) = ((α, β) = params(d); α > 1.0 ? β / (α - 1.0) : Inf) +mean(d::InverseGamma) = ((α, θ) = params(d); α > 1.0 ? θ / (α - 1.0) : Inf) mode(d::InverseGamma) = scale(d) / (shape(d) + 1.0) function var(d::InverseGamma) - (α, β) = params(d) - α > 2.0 ? β^2 / ((α - 1.0)^2 * (α - 2.0)) : Inf + (α, θ) = params(d) + α > 2.0 ? θ^2 / ((α - 1.0)^2 * (α - 2.0)) : Inf end function skewness(d::InverseGamma) @@ -45,8 +45,8 @@ function kurtosis(d::InverseGamma) end function entropy(d::InverseGamma) - (α, β) = params(d) - α + lgamma(α) - (1.0 + α) * digamma(α) + log(β) + (α, θ) = params(d) + α + lgamma(α) - (1.0 + α) * digamma(α) + log(θ) end @@ -55,8 +55,8 @@ end pdf(d::InverseGamma, x::Float64) = exp(logpdf(d, x)) function logpdf(d::InverseGamma, x::Float64) - (α, β) = params(d) - α * log(β) - lgamma(α) - (α + 1.0) * log(x) - β / x + (α, θ) = params(d) + α * log(θ) - lgamma(α) - (α + 1.0) * log(x) - θ / x end cdf(d::InverseGamma, x::Float64) = ccdf(d.invd, 1.0 / x) @@ -71,12 +71,12 @@ invlogccdf(d::InverseGamma, p::Float64) = 1.0 / invlogcdf(d.invd, p) function mgf(d::InverseGamma, t::Real) (a, b) = params(d) - @compat t == zero(t) ? one(Float64(t)) : 2.0*(-b*t)^(0.5a) / gamma(a) * besselk(a, sqrt(-4.0*b*t)) + t == zero(t) ? 1.0 : 2.0*(-b*t)^(0.5a) / gamma(a) * besselk(a, sqrt(-4.0*b*t)) end function cf(d::InverseGamma, t::Real) (a, b) = params(d) - @compat t == zero(t) ? complex(one(Float64(t))) : 2.0*(-im*b*t)^(0.5a) / gamma(a) * besselk(a, sqrt(-4.0*im*b*t)) + t == zero(t) ? 1.0+0.0im : 2.0*(-im*b*t)^(0.5a) / gamma(a) * besselk(a, sqrt(-4.0*im*b*t)) end @@ -92,4 +92,3 @@ function _rand!(d::InverseGamma, A::AbstractArray) end A end - diff --git a/src/univariate/continuous/inversegaussian.jl b/src/univariate/continuous/inversegaussian.jl index ad51b0998..0eed7918b 100644 --- a/src/univariate/continuous/inversegaussian.jl +++ b/src/univariate/continuous/inversegaussian.jl @@ -3,10 +3,9 @@ immutable InverseGaussian <: ContinuousUnivariateDistribution λ::Float64 function InverseGaussian(μ::Real, λ::Real) - (μ > zero(μ) && λ > zero(λ)) || error("InverseGaussian's μ and λ must be positive") - @compat new(Float64(μ), Float64(λ)) + @check_args(InverseGaussian, μ > zero(μ) && λ > zero(λ)) + new(μ, λ) end - InverseGaussian(μ::Real) = InverseGaussian(μ, 1.0) InverseGaussian() = new(1.0, 1.0) end diff --git a/src/univariate/continuous/laplace.jl b/src/univariate/continuous/laplace.jl index 5cd1f63c6..1e0fbbcf4 100644 --- a/src/univariate/continuous/laplace.jl +++ b/src/univariate/continuous/laplace.jl @@ -1,13 +1,9 @@ immutable Laplace <: ContinuousUnivariateDistribution μ::Float64 - β::Float64 + θ::Float64 - function Laplace(μ::Real, β::Real) - β > zero(β) || error("Laplace's scale must be positive") - @compat new(Float64(μ), Float64(β)) - end - - @compat Laplace(μ::Real) = new(Float64(μ), 1.0) + Laplace(μ::Real, θ::Real) = (@check_args(Laplace, θ > zero(θ)); new(μ, θ)) + Laplace(μ::Real) = new(μ, 1.0) Laplace() = new(0.0, 1.0) end @@ -19,8 +15,8 @@ typealias Biexponential Laplace #### Parameters location(d::Laplace) = d.μ -scale(d::Laplace) = d.β -params(d::Laplace) = (d.μ, d.β) +scale(d::Laplace) = d.θ +params(d::Laplace) = (d.μ, d.θ) #### Statistics @@ -29,18 +25,18 @@ mean(d::Laplace) = d.μ median(d::Laplace) = d.μ mode(d::Laplace) = d.μ -var(d::Laplace) = 2.0 * d.β^2 -std(d::Laplace) = sqrt2 * d.β +var(d::Laplace) = 2.0 * d.θ^2 +std(d::Laplace) = sqrt2 * d.θ skewness(d::Laplace) = 0.0 kurtosis(d::Laplace) = 3.0 -entropy(d::Laplace) = log(2.0 * d.β) + 1.0 +entropy(d::Laplace) = log(2.0 * d.θ) + 1.0 #### Evaluations -zval(d::Laplace, x::Float64) = (x - d.μ) / d.β -xval(d::Laplace, z::Float64) = d.μ + z * d.β +zval(d::Laplace, x::Float64) = (x - d.μ) / d.θ +xval(d::Laplace, z::Float64) = d.μ + z * d.θ pdf(d::Laplace, x::Float64) = 0.5 * exp(-abs(zval(d, x))) / scale(d) logpdf(d::Laplace, x::Float64) = - (abs(zval(d, x)) + log(2.0 * scale(d))) @@ -56,25 +52,25 @@ invlogcdf(d::Laplace, lp::Float64) = lp < loghalf ? xval(d, logtwo + lp) : xval( invlogccdf(d::Laplace, lp::Float64) = lp > loghalf ? xval(d, logtwo + log1mexp(lp)) : xval(d, -(logtwo + lp)) function gradlogpdf(d::Laplace, x::Float64) - μ, β = params(d) + μ, θ = params(d) x == μ && error("Gradient is undefined at the location point") - g = 1.0 / β + g = 1.0 / θ x > μ ? -g : g end function mgf(d::Laplace, t::Real) - st = d.β * t + st = d.θ * t exp(t * d.μ) / ((1.0 - st) * (1.0 + st)) end function cf(d::Laplace, t::Real) - st = d.β * t + st = d.θ * t cis(t * d.μ) / (1+st*st) end #### Sampling -rand(d::Laplace) = d.μ + d.β*randexp()*ifelse(rand(Bool), 1, -1) +rand(d::Laplace) = d.μ + d.θ*randexp()*ifelse(rand(Bool), 1, -1) #### Fitting @@ -83,5 +79,3 @@ function fit_mle(::Type{Laplace}, x::Array) a = median(x) Laplace(a, mad(x, a)) end - - diff --git a/src/univariate/continuous/levy.jl b/src/univariate/continuous/levy.jl index 2775cc5fc..dab2519ac 100644 --- a/src/univariate/continuous/levy.jl +++ b/src/univariate/continuous/levy.jl @@ -1,23 +1,19 @@ immutable Levy <: ContinuousUnivariateDistribution μ::Float64 - c::Float64 - - function Levy(μ::Real, c::Real) - c >= zero(c) || error("scale must be non-negative") - @compat new(Float64(μ), Float64(c)) - end + σ::Float64 + Levy(μ::Real, σ::Real) = (@check_args(Levy, σ > zero(σ)); new(μ, σ)) Levy(μ::Real) = new(μ, 1.0) Levy() = new(0.0, 1.0) end -@distr_support Levy location(d) Inf +@distr_support Levy d.μ Inf #### Parameters location(d::Levy) = d.μ -params(d::Levy) = (d.μ, d.c) +params(d::Levy) = (d.μ, d.σ) #### Statistics @@ -27,48 +23,41 @@ var(d::Levy) = Inf skewness(d::Levy) = NaN kurtosis(d::Levy) = NaN -mode(d::Levy) = d.c / 3.0 + d.μ +mode(d::Levy) = d.σ / 3.0 + d.μ -function entropy(d::Levy) - c = scale(d) - (1.0 - 3.0 * digamma(1.0) + log(16.0 * pi * c * c)) / 2.0 -end +entropy(d::Levy) = (1.0 - 3.0 * digamma(1.0) + log(16π * d.σ^2)) / 2.0 -function median(d::Levy) - μ, c = params(d) - μ + c / (2.0 * erfcinv(0.5)^2) -end +median(d::Levy) = d.μ + d.σ / 0.4549364231195728 # 0.454... = (2.0 * erfcinv(0.5)^2) #### Evaluation function pdf(d::Levy, x::Float64) - μ, c = params(d) + μ, σ = params(d) z = x - μ - (sqrt(c) / sqrt2π) * exp((-c) / (2.0 * z)) / z^1.5 + (sqrt(σ) / sqrt2π) * exp((-σ) / (2.0 * z)) / z^1.5 end function logpdf(d::Levy, x::Float64) - μ, c = params(d) + μ, σ = params(d) z = x - μ - 0.5 * (log(c) - log2π - c / z - 3.0 * log(z)) + 0.5 * (log(σ) - log2π - σ / z - 3.0 * log(z)) end -cdf(d::Levy, x::Float64) = erfc(sqrt(d.c / (2.0 * (x - d.μ)))) -ccdf(d::Levy, x::Float64) = erf(sqrt(d.c / (2.0 * (x - d.μ)))) +cdf(d::Levy, x::Float64) = erfc(sqrt(d.σ / (2.0 * (x - d.μ)))) +ccdf(d::Levy, x::Float64) = erf(sqrt(d.σ / (2.0 * (x - d.μ)))) -quantile(d::Levy, p::Float64) = d.μ + d.c / (2.0 * erfcinv(p)^2) -cquantile(d::Levy, p::Float64) = d.μ + d.c / (2.0 * erfinv(p)^2) +quantile(d::Levy, p::Float64) = d.μ + d.σ / (2.0 * erfcinv(p)^2) +cquantile(d::Levy, p::Float64) = d.μ + d.σ / (2.0 * erfinv(p)^2) mgf(d::Levy, t::Real) = t == zero(t) ? 1.0 : NaN function cf(d::Levy, t::Real) - μ, c = params(d) - exp(im * μ * t - sqrt(-2.0 * im * c * t)) + μ, σ = params(d) + exp(im * μ * t - sqrt(-2.0 * im * σ * t)) end #### Sampling -rand(d::Levy) = d.μ + d.c / randn()^2 - +rand(d::Levy) = d.μ + d.σ / randn()^2 diff --git a/src/univariate/continuous/logistic.jl b/src/univariate/continuous/logistic.jl index 56121255e..6b784b0f1 100644 --- a/src/univariate/continuous/logistic.jl +++ b/src/univariate/continuous/logistic.jl @@ -1,13 +1,9 @@ immutable Logistic <: ContinuousUnivariateDistribution μ::Float64 - β::Float64 + θ::Float64 - function Logistic(μ::Float64, β::Float64) - β > zero(β) || error("Logistic: scale must be positive") - @compat new(Float64(μ), Float64(β)) - end - - @compat Logistic(μ::Real) = new(Float64(μ), 1.0) + Logistic(μ::Real, θ::Real) = (@check_args(Logistic, θ > zero(θ)); new(μ, θ)) + Logistic(μ::Real) = new(μ, 1.0) Logistic() = new(0.0, 1.0) end @@ -17,9 +13,9 @@ end #### Parameters location(d::Logistic) = d.μ -scale(d::Logistic) = d.β +scale(d::Logistic) = d.θ -params(d::Logistic) = (d.μ, d.β) +params(d::Logistic) = (d.μ, d.θ) #### Statistics @@ -28,21 +24,21 @@ mean(d::Logistic) = d.μ median(d::Logistic) = d.μ mode(d::Logistic) = d.μ -std(d::Logistic) = π * d.β / sqrt3 -var(d::Logistic) = (π * d.β)^2 / 3.0 +std(d::Logistic) = π * d.θ / sqrt3 +var(d::Logistic) = (π * d.θ)^2 / 3.0 skewness(d::Logistic) = 0.0 kurtosis(d::Logistic) = 1.2 -entropy(d::Logistic) = log(d.β) + 2.0 +entropy(d::Logistic) = log(d.θ) + 2.0 #### Evaluation -zval(d::Logistic, x::Float64) = (x - d.μ) / d.β -xval(d::Logistic, z::Float64) = d.μ + z * d.β +zval(d::Logistic, x::Float64) = (x - d.μ) / d.θ +xval(d::Logistic, z::Float64) = d.μ + z * d.θ -pdf(d::Logistic, x::Float64) = (e = exp(-zval(d, x)); e / (d.β * (1.0 + e)^2)) -logpdf(d::Logistic, x::Float64) = (u = -abs(zval(d, x)); u - 2.0 * log1pexp(u) - log(d.β)) +pdf(d::Logistic, x::Float64) = (e = exp(-zval(d, x)); e / (d.θ * (1.0 + e)^2)) +logpdf(d::Logistic, x::Float64) = (u = -abs(zval(d, x)); u - 2.0 * log1pexp(u) - log(d.θ)) cdf(d::Logistic, x::Float64) = logistic(zval(d, x)) ccdf(d::Logistic, x::Float64) = logistic(-zval(d, x)) @@ -50,19 +46,19 @@ logcdf(d::Logistic, x::Float64) = -log1pexp(-zval(d, x)) logccdf(d::Logistic, x::Float64) = -log1pexp(zval(d, x)) quantile(d::Logistic, p::Float64) = xval(d, logit(p)) -cquantile(d::Logistic, p::Float64) = xval(d, -logit(p)) -invlogcdf(d::Logistic, lp::Float64) = xval(d, -logexpm1(-lp)) -invlogccdf(d::Logistic, lp::Float64) = xval(d, logexpm1(-lp)) +cquantile(d::Logistic, p::Float64) = xval(d, -logit(p)) +invlogcdf(d::Logistic, lp::Float64) = xval(d, -logexpm1(-lp)) +invlogccdf(d::Logistic, lp::Float64) = xval(d, logexpm1(-lp)) function gradlogpdf(d::Logistic, x::Float64) e = exp(-zval(d, x)) - ((2.0 * e) / (1.0 + e) - 1.0) / d.β + ((2.0 * e) / (1.0 + e) - 1.0) / d.θ end -mgf(d::Logistic, t::Real) = exp(t * d.μ) / sinc(d.β * t) +mgf(d::Logistic, t::Real) = exp(t * d.μ) / sinc(d.θ * t) function cf(d::Logistic, t::Real) - a = (π * t) * d.β + a = (π * t) * d.θ a == zero(a) ? complex(one(a)) : cis(t * d.μ) * (a / sinh(a)) end @@ -70,4 +66,3 @@ end #### Sampling rand(d::Logistic) = quantile(d, rand()) - diff --git a/src/univariate/continuous/lognormal.jl b/src/univariate/continuous/lognormal.jl index bb380cdef..91dd0a1e7 100644 --- a/src/univariate/continuous/lognormal.jl +++ b/src/univariate/continuous/lognormal.jl @@ -2,17 +2,14 @@ immutable LogNormal <: ContinuousUnivariateDistribution μ::Float64 σ::Float64 - function LogNormal(μ::Real, σ::Real) - σ > zero(σ) || error("σ must be positive") - @compat new(Float64(μ), Float64(σ)) - end - - LogNormal(μ::Real) = @compat new(Float64(μ), 1.0) + LogNormal(μ::Real, σ::Real) = (@check_args(LogNormal, σ > zero(σ)); new(μ, σ)) + LogNormal(μ::Real) = new(μ, 1.0) LogNormal() = new(0.0, 1.0) end @distr_support LogNormal 0.0 Inf + #### Parameters params(d::LogNormal) = (d.μ, d.σ) diff --git a/src/univariate/continuous/noncentralbeta.jl b/src/univariate/continuous/noncentralbeta.jl index 7aa281c41..7d2b8d415 100644 --- a/src/univariate/continuous/noncentralbeta.jl +++ b/src/univariate/continuous/noncentralbeta.jl @@ -2,15 +2,24 @@ immutable NoncentralBeta <: ContinuousUnivariateDistribution α::Float64 β::Float64 λ::Float64 - function NoncentralBeta(a::Real, b::Real, nc::Real) - a > 0.0 && b > 0.0 && nc >= 0.0 || - error("alpha and beta must be > 0 and ncp >= 0") - @compat new(Float64(a), Float64(b), Float64(nc)) + + function NoncentralBeta(α::Real, β::Real, λ::Real) + @check_args(NoncentralBeta, α > zero(α) && β > zero(β)) + @check_args(NoncentralBeta, λ >= zero(λ)) + new(α, β, λ) end end @distr_support NoncentralBeta 0.0 1.0 + +### Parameters + +params(d::NoncentralBeta) = (d.α, d.β, d.λ) + + +### Evaluation & Sampling + # TODO: add mean and var @_delegate_statsfuns NoncentralBeta nbeta α β λ @@ -18,5 +27,5 @@ end function rand(d::NoncentralBeta) a = rand(NoncentralChisq(2.0 * d.α, d.β)) b = rand(Chisq(2.0 * d.β)) - a / (a+b) + a / (a + b) end diff --git a/src/univariate/continuous/noncentralchisq.jl b/src/univariate/continuous/noncentralchisq.jl index 72669e47d..4373805ad 100644 --- a/src/univariate/continuous/noncentralchisq.jl +++ b/src/univariate/continuous/noncentralchisq.jl @@ -1,28 +1,38 @@ immutable NoncentralChisq <: ContinuousUnivariateDistribution - df::Float64 + ν::Float64 λ::Float64 - function NoncentralChisq(d::Real, nc::Real) - d >= zero(d) && nc >= zero(nc) || error("df and ncp must be non-negative") - @compat new(Float64(d), Float64(nc)) + function NoncentralChisq(ν::Real, λ::Real) + @check_args(NoncentralChisq, ν > zero(ν)) + @check_args(NoncentralChisq, λ >= zero(λ)) + new(ν, λ) end end @distr_support NoncentralChisq 0.0 Inf -mean(d::NoncentralChisq) = d.df + d.λ -var(d::NoncentralChisq) = 2.0*(d.df + 2.0*d.λ) -skewness(d::NoncentralChisq) = 2.0*sqrt2*(d.df + 3.0*d.λ)/sqrt(d.df + 2.0*d.λ)^3 -kurtosis(d::NoncentralChisq) = 12.0*(d.df + 4.0*d.λ)/(d.df + 2.0*d.λ)^2 +### Parameters -function mgf(d::NoncentralChisq, t::Real) - k = d.df - exp(d.λ * t/(1.0 - 2.0 * t))*(1.0 - 2.0 * t)^(-k / 2.0) +params(d::NoncentralChisq) = (d.ν, d.λ) + + +### Statistics + +mean(d::NoncentralChisq) = d.ν + d.λ +var(d::NoncentralChisq) = 2.0*(d.ν + 2.0*d.λ) +skewness(d::NoncentralChisq) = 2.0*sqrt2*(d.ν + 3.0*d.λ)/sqrt(d.ν + 2.0*d.λ)^3 +kurtosis(d::NoncentralChisq) = 12.0*(d.ν + 4.0*d.λ)/(d.ν + 2.0*d.λ)^2 + +function mgf(d::NoncentralChisq, t::Float64) + exp(d.λ * t/(1.0 - 2.0 * t))*(1.0 - 2.0 * t)^(-d.ν / 2.0) end -function cf(d::NoncentralChisq, t::Real) - cis(d.λ * t/(1.0 - 2.0 * im * t))*(1.0 - 2.0 * im * t)^(-d.df / 2.0) +function cf(d::NoncentralChisq, t::Float64) + cis(d.λ * t/(1.0 - 2.0 * im * t))*(1.0 - 2.0 * im * t)^(-d.ν / 2.0) end -@_delegate_statsfuns NoncentralChisq nchisq df λ -rand(d::NoncentralChisq) = StatsFuns.Rmath.nchisqrand(d.df, d.λ) +### Evaluation & Sampling + +@_delegate_statsfuns NoncentralChisq nchisq ν λ + +rand(d::NoncentralChisq) = StatsFuns.Rmath.nchisqrand(d.ν, d.λ) diff --git a/src/univariate/continuous/noncentralf.jl b/src/univariate/continuous/noncentralf.jl index b10768991..5ae30d487 100644 --- a/src/univariate/continuous/noncentralf.jl +++ b/src/univariate/continuous/noncentralf.jl @@ -1,27 +1,38 @@ immutable NoncentralF <: ContinuousUnivariateDistribution - ndf::Float64 - ddf::Float64 + ν1::Float64 + ν2::Float64 λ::Float64 - function NoncentralF(n::Real, d::Real, nc::Real) - n > zero(n) && d > zero(d) && nc >= zero(nc) || - error("ndf and ddf must be > 0 and λ >= 0") - @compat new(Float64(n), Float64(d), Float64(nc)) + function NoncentralF(ν1::Real, ν2::Real, λ::Real) + @check_args(NoncentralF, ν1 > zero(ν1) && ν2 > zero(ν2)) + @check_args(NoncentralF, λ >= zero(λ)) + new(ν1, ν2, λ) end end @distr_support NoncentralF 0.0 Inf -mean(d::NoncentralF) = d.ddf > 2.0 ? d.ddf / (d.ddf - 2.0) * (d.ndf + d.λ) / d.ndf : NaN -var(d::NoncentralF) = d.ddf > 4.0 ? 2.0 * d.ddf^2 * - ((d.ndf+d.λ)^2 + (d.ddf - 2.0)*(d.ndf + 2.0*d.λ)) / - (d.ndf * (d.ddf - 2.0)^2 * (d.ddf - 4.0)) : NaN +### Parameters -@_delegate_statsfuns NoncentralF nfdist ndf ddf λ +params(d::NoncentralF) = (d.ν1, d.ν2, d.λ) + + +### Statistics + +mean(d::NoncentralF) = d.ν2 > 2.0 ? d.ν2 / (d.ν2 - 2.0) * (d.ν1 + d.λ) / d.ν1 : NaN + +var(d::NoncentralF) = d.ν2 > 4.0 ? 2.0 * d.ν2^2 * + ((d.ν1+d.λ)^2 + (d.ν2 - 2.0)*(d.ν1 + 2.0*d.λ)) / + (d.ν1 * (d.ν2 - 2.0)^2 * (d.ν2 - 4.0)) : NaN + + +### Evaluation & Sampling + +@_delegate_statsfuns NoncentralF nfdist ν1 ν2 λ function rand(d::NoncentralF) - rn = rand(NoncentralChisq(d.ndf,d.λ)) / d.ndf - rd = rand(Chisq(d.ddf)) / d.ddf - rn / rd + r1 = rand(NoncentralChisq(d.ν1,d.λ)) / d.ν1 + r2 = rand(Chisq(d.ν2)) / d.ν2 + r1 / r2 end diff --git a/src/univariate/continuous/noncentralt.jl b/src/univariate/continuous/noncentralt.jl index b82e872a5..51b85086f 100644 --- a/src/univariate/continuous/noncentralt.jl +++ b/src/univariate/continuous/noncentralt.jl @@ -1,33 +1,41 @@ immutable NoncentralT <: ContinuousUnivariateDistribution - df::Float64 + ν::Float64 λ::Float64 - function NoncentralT(d::Real, nc::Real) - d >= zero(d) && nc >= zero(nc) || error("df and λ must be non-negative") - @compat new(Float64(d), Float64(nc)) + function NoncentralT(ν::Real, λ::Real) + @check_args(NoncentralT, ν > zero(ν)) + @check_args(NoncentralT, λ >= zero(λ)) + new(ν, λ) end end @distr_support NoncentralT -Inf Inf +### Parameters + +params(d::NoncentralT) = (d.ν, d.λ) + + +### Statistics + function mean(d::NoncentralT) - if d.df > 1.0 - if isinf(d.df) - d.λ - else - sqrt(0.5*d.df)*d.λ*gamma(0.5*(d.df-1))/gamma(0.5*d.df) - end + if d.ν > 1.0 + isinf(d.ν) ? d.λ : + sqrt(0.5*d.ν) * d.λ * gamma(0.5*(d.ν-1)) / gamma(0.5*d.ν) else NaN end end -var(d::NoncentralT) = d.df > 2.0 ? d.df*(1+d.λ^2)/(d.df-2.0) - mean(d)^2 : NaN +var(d::NoncentralT) = d.ν > 2.0 ? d.ν*(1+d.λ^2)/(d.ν-2.0) - mean(d)^2 : NaN + + +### Evaluation & Sampling -@_delegate_statsfuns NoncentralT ntdist df λ +@_delegate_statsfuns NoncentralT ntdist ν λ function rand(d::NoncentralT) z = randn() - v = rand(Chisq(d.df)) - (z+d.λ)/sqrt(v/d.df) + v = rand(Chisq(d.ν)) + (z+d.λ)/sqrt(v/d.ν) end diff --git a/src/univariate/continuous/normal.jl b/src/univariate/continuous/normal.jl index 7da55bcf3..f7158da18 100644 --- a/src/univariate/continuous/normal.jl +++ b/src/univariate/continuous/normal.jl @@ -2,12 +2,8 @@ immutable Normal <: ContinuousUnivariateDistribution μ::Float64 σ::Float64 - function Normal(μ::Real, σ::Real) - σ > zero(σ) || error("std.dev. must be positive") - @compat new(Float64(μ), Float64(σ)) - end - - @compat Normal(μ::Real) = Normal(Float64(μ), 1.0) + Normal(μ::Real, σ::Real) = (@check_args(Normal, σ > zero(σ)); new(μ, σ)) + Normal(μ::Real) = Normal(μ, 1.0) Normal() = Normal(0.0, 1.0) end @@ -57,8 +53,6 @@ immutable NormalStats <: SufficientStats m::Float64 # (weighted) mean of x s2::Float64 # (weighted) sum of (x - μ)^2 tw::Float64 # total sample weight - - @compat NormalStats(s::Real, m::Real, s2::Real, tw::Real) = new(Float64(s), Float64(m), Float64(s2), Float64(tw)) end function suffstats{T<:Real}(::Type{Normal}, x::AbstractArray{T}) @@ -120,7 +114,7 @@ function suffstats{T<:Real}(g::NormalKnownMu, x::AbstractArray{T}) for i = 2:length(x) @inbounds s2 += abs2(x[i] - μ) end - @compat NormalKnownMuStats(g.μ, s2, Float64(length(x))) + NormalKnownMuStats(g.μ, s2, length(x)) end function suffstats{T<:Real}(g::NormalKnownMu, x::AbstractArray{T}, w::AbstractArray{Float64}) diff --git a/src/univariate/continuous/normalcanon.jl b/src/univariate/continuous/normalcanon.jl index 66b1581f8..3714ada28 100644 --- a/src/univariate/continuous/normalcanon.jl +++ b/src/univariate/continuous/normalcanon.jl @@ -1,55 +1,54 @@ ## Canonical Form of Normal distribution immutable NormalCanon <: ContinuousUnivariateDistribution - h::Float64 # σ^(-2) * μ - prec::Float64 # σ^(-2) + η::Float64 # σ^(-2) * μ + λ::Float64 # σ^(-2) μ::Float64 # μ - function NormalCanon(h::Float64, prec::Float64) - prec > 0. || throw(ArgumentError("prec should be positive.")) - new(h, prec, h / prec) + function NormalCanon(η::Real, λ::Real) + @check_args(NormalCanon, λ > zero(λ)) + new(η, λ, η / λ) end - - @compat NormalCanon(h::Real, prec::Real) = NormalCanon(Float64(h), Float64(prec)) NormalCanon() = new(0.0, 1.0, 0.0) end @distr_support NormalCanon -Inf Inf + ## conversion between Normal and NormalCanon -Base.convert(::Type{Normal}, cf::NormalCanon) = Normal(cf.μ, 1.0 / sqrt(cf.prec)) -Base.convert(::Type{NormalCanon}, d::Normal) = (J = 1.0 / abs2(σ); NormalCanon(J * d.μ, J)) +Base.convert(::Type{Normal}, d::NormalCanon) = Normal(d.μ, 1.0 / sqrt(d.λ)) +Base.convert(::Type{NormalCanon}, d::Normal) = (λ = 1.0 / σ^2; NormalCanon(λ * d.μ, λ)) canonform(d::Normal) = convert(NormalCanon, d) #### Parameters -params(d::NormalCanon) = (d.h, d.prec) +params(d::NormalCanon) = (d.η, d.λ) #### Statistics -mean(cf::NormalCanon) = cf.μ -median(cf::NormalCanon) = mean(cf) -mode(cf::NormalCanon) = mean(cf) +mean(d::NormalCanon) = d.μ +median(d::NormalCanon) = mean(d) +mode(d::NormalCanon) = mean(d) -skewness(cf::NormalCanon) = 0.0 -kurtosis(cf::NormalCanon) = 0.0 +skewness(d::NormalCanon) = 0.0 +kurtosis(d::NormalCanon) = 0.0 -var(cf::NormalCanon) = 1.0 / cf.prec -std(cf::NormalCanon) = sqrt(var(cf)) +var(d::NormalCanon) = 1.0 / d.λ +std(d::NormalCanon) = sqrt(var(d)) -entropy(cf::NormalCanon) = 0.5 * (log2π + 1. - log(cf.prec)) +entropy(d::NormalCanon) = 0.5 * (log2π + 1.0 - log(d.λ)) #### Evaluation -pdf(d::NormalCanon, x::Float64) = (sqrt(d.prec) / sqrt2π) * exp(-0.5 * d.prec * abs2(x - d.μ)) -logpdf(d::NormalCanon, x::Float64) = 0.5 * (log(d.prec) - log2π - d.prec * abs2(x - d.μ)) +pdf(d::NormalCanon, x::Float64) = (sqrt(d.λ) / sqrt2π) * exp(-0.5 * d.λ * abs2(x - d.μ)) +logpdf(d::NormalCanon, x::Float64) = 0.5 * (log(d.λ) - log2π - d.λ * abs2(x - d.μ)) -zval(d::NormalCanon, x::Float64) = (x - d.μ) * sqrt(d.prec) -xval(d::NormalCanon, z::Float64) = d.μ + z / sqrt(d.prec) +zval(d::NormalCanon, x::Float64) = (x - d.μ) * sqrt(d.λ) +xval(d::NormalCanon, z::Float64) = d.μ + z / sqrt(d.λ) cdf(d::NormalCanon, x::Float64) = normcdf(zval(d,x)) ccdf(d::NormalCanon, x::Float64) = normccdf(zval(d,x)) @@ -64,5 +63,5 @@ invlogccdf(d::NormalCanon, lp::Float64) = xval(d, norminvlogccdf(lp)) #### Sampling -rand(cf::NormalCanon) = cf.μ + randn() / sqrt(cf.prec) +rand(cf::NormalCanon) = cf.μ + randn() / sqrt(cf.λ) rand!{T<:Real}(cf::NormalCanon, r::AbstractArray{T}) = rand!(convert(Normal, cf), r) diff --git a/src/univariate/continuous/pareto.jl b/src/univariate/continuous/pareto.jl index c426c8c4a..168de7efe 100644 --- a/src/univariate/continuous/pareto.jl +++ b/src/univariate/continuous/pareto.jl @@ -1,36 +1,35 @@ immutable Pareto <: ContinuousUnivariateDistribution α::Float64 - β::Float64 + θ::Float64 - function Pareto(α::Real, β::Real) - (α > zero(α) && β > zero(β)) || error("Pareto: shape and scale must be positive") - @compat new(Float64(α), Float64(β)) + function Pareto(α::Real, θ::Real) + @check_args(Pareto, α > zero(α) && θ > zero(θ)) + new(α, θ) end - Pareto(α::Real) = Pareto(α, 1.0) Pareto() = new(1.0, 1.0) end -@distr_support Pareto d.β Inf +@distr_support Pareto d.θ Inf #### Parameters shape(d::Pareto) = d.α -scale(d::Pareto) = d.β +scale(d::Pareto) = d.θ -params(d::Pareto) = (d.α, d.β) +params(d::Pareto) = (d.α, d.θ) #### Statistics -mean(d::Pareto) = ((α, β) = params(d); α > 1.0 ? α * β / (α - 1.0) : Inf) -median(d::Pareto) = ((α, β) = params(d); β * 2.0 ^ (1.0 / α)) -mode(d::Pareto) = d.β +mean(d::Pareto) = ((α, θ) = params(d); α > 1.0 ? α * θ / (α - 1.0) : Inf) +median(d::Pareto) = ((α, θ) = params(d); θ * 2.0 ^ (1.0 / α)) +mode(d::Pareto) = d.θ function var(d::Pareto) - (α, β) = params(d) - α > 2.0 ? (β^2 * α) / ((α - 1.0)^2 * (α - 2.0)) : Inf + (α, θ) = params(d) + α > 2.0 ? (θ^2 * α) / ((α - 1.0)^2 * (α - 2.0)) : Inf end function skewness(d::Pareto) @@ -43,42 +42,39 @@ function kurtosis(d::Pareto) α > 4.0 ? (6.0 * (α^3 + α^2 - 6.0 * α - 2.0)) / (α * (α - 3.0) * (α - 4.0)) : NaN end -entropy(d::Pareto) = ((α, β) = params(d); log(β / α) + 1.0 / α + 1.0) +entropy(d::Pareto) = ((α, θ) = params(d); log(θ / α) + 1.0 / α + 1.0) #### Evaluation function pdf(d::Pareto, x::Float64) - (α, β) = params(d) - x >= β ? α * (β / x)^α * (1.0 / x) : 0.0 + (α, θ) = params(d) + x >= θ ? α * (θ / x)^α * (1.0 / x) : 0.0 end function logpdf(d::Pareto, x::Float64) - (α, β) = params(d) - x >= β ? log(α) + α * log(β) - (α + 1.0) * log(x) : -Inf + (α, θ) = params(d) + x >= θ ? log(α) + α * log(θ) - (α + 1.0) * log(x) : -Inf end function ccdf(d::Pareto, x::Float64) - (α, β) = params(d) - x >= β ? (β / x)^α : 1.0 + (α, θ) = params(d) + x >= θ ? (θ / x)^α : 1.0 end cdf(d::Pareto, x::Float64) = 1.0 - ccdf(d, x) function logccdf(d::Pareto, x::Float64) - (α, β) = params(d) - x >= β ? α * log(β / x) : 0.0 + (α, θ) = params(d) + x >= θ ? α * log(θ / x) : 0.0 end logcdf(d::Pareto, x::Float64) = log1p(-ccdf(d, x)) -cquantile(d::Pareto, p::Float64) = d.β / p^(1.0 / d.α) +cquantile(d::Pareto, p::Float64) = d.θ / p^(1.0 / d.α) quantile(d::Pareto, p::Float64) = cquantile(d, 1.0 - p) #### Sampling -rand(d::Pareto) = d.β * exp(randexp() / d.α) - - - +rand(d::Pareto) = d.θ * exp(randexp() / d.α) diff --git a/src/univariate/continuous/rayleigh.jl b/src/univariate/continuous/rayleigh.jl index 65add7a37..9f7481e43 100644 --- a/src/univariate/continuous/rayleigh.jl +++ b/src/univariate/continuous/rayleigh.jl @@ -1,11 +1,7 @@ immutable Rayleigh <: ContinuousUnivariateDistribution σ::Float64 - function Rayleigh(σ::Real) - σ > zero(σ) || error("Rayleigh: σ must be positive") - @compat new(Float64(σ)) - end - + Rayleigh(σ::Real) = (@check_args(Rayleigh, σ > zero(σ)); new(σ)) Rayleigh() = new(1.0) end @@ -57,4 +53,3 @@ quantile(d::Rayleigh, p::Float64) = sqrt(-2.0 * d.σ^2 * log1p(-p)) #### Sampling rand(d::Rayleigh) = d.σ * sqrt(2.0 * randexp()) - diff --git a/src/univariate/continuous/symtriangular.jl b/src/univariate/continuous/symtriangular.jl index 51c53d7f8..e58b6e770 100644 --- a/src/univariate/continuous/symtriangular.jl +++ b/src/univariate/continuous/symtriangular.jl @@ -1,25 +1,24 @@ immutable SymTriangularDist <: ContinuousUnivariateDistribution μ::Float64 - s::Float64 + σ::Float64 - function SymTriangularDist(μ::Real, s::Real) - s > zero(s) || error("SymTriangular: scale must be positive") - @compat new(Float64(μ), Float64(s)) + function SymTriangularDist(μ::Real, σ::Real) + @check_args(SymTriangularDist, σ > zero(σ)) + new(μ, σ) end - - @compat SymTriangularDist(μ::Real) = new(Float64(μ), 1.0) + SymTriangularDist(μ::Real) = new(μ, 1.0) SymTriangularDist() = new(0.0, 1.0) end -@distr_support SymTriangularDist d.μ - d.s d.μ + d.s +@distr_support SymTriangularDist d.μ - d.σ d.μ + d.σ #### Parameters location(d::SymTriangularDist) = d.μ -scale(d::SymTriangularDist) = d.s +scale(d::SymTriangularDist) = d.σ -params(d::SymTriangularDist) = (d.μ, d.s) +params(d::SymTriangularDist) = (d.μ, d.σ) #### Statistics @@ -28,17 +27,17 @@ mean(d::SymTriangularDist) = d.μ median(d::SymTriangularDist) = d.μ mode(d::SymTriangularDist) = d.μ -var(d::SymTriangularDist) = d.s^2 / 6.0 +var(d::SymTriangularDist) = d.σ^2 / 6.0 skewness(d::SymTriangularDist) = 0.0 kurtosis(d::SymTriangularDist) = -0.6 -entropy(d::SymTriangularDist) = 0.5 + log(d.s) +entropy(d::SymTriangularDist) = 0.5 + log(d.σ) #### Evaluation -zval(d::SymTriangularDist, x::Float64) = (x - d.μ) / d.s -xval(d::SymTriangularDist, z::Float64) = d.μ + z * d.s +zval(d::SymTriangularDist, x::Float64) = (x - d.μ) / d.σ +xval(d::SymTriangularDist, z::Float64) = d.μ + z * d.σ pdf(d::SymTriangularDist, x::Float64) = insupport(d, x) ? (1.0 - abs(zval(d, x))) / scale(d) : 0.0 @@ -46,38 +45,38 @@ pdf(d::SymTriangularDist, x::Float64) = insupport(d, x) ? (1.0 - abs(zval(d, x)) logpdf(d::SymTriangularDist, x::Float64) = insupport(d, x) ? log((1.0 - abs(zval(d, x))) / scale(d)) : -Inf function cdf(d::SymTriangularDist, x::Float64) - (μ, s) = params(d) - x <= μ - s ? 0.0 : + (μ, σ) = params(d) + x <= μ - σ ? 0.0 : x <= μ ? 0.5 * (1.0 + zval(d, x))^2 : - x < μ + s ? 1.0 - 0.5 * (1.0 - zval(d, x))^2 : 1.0 + x < μ + σ ? 1.0 - 0.5 * (1.0 - zval(d, x))^2 : 1.0 end function ccdf(d::SymTriangularDist, x::Float64) - (μ, s) = params(d) - x <= μ - s ? 1.0 : + (μ, σ) = params(d) + x <= μ - σ ? 1.0 : x <= μ ? 1.0 - 0.5 * (1.0 + zval(d, x))^2 : - x < μ + s ? 0.5 * (1.0 - zval(d, x))^2 : 0.0 + x < μ + σ ? 0.5 * (1.0 - zval(d, x))^2 : 0.0 end function logcdf(d::SymTriangularDist, x::Float64) - (μ, s) = params(d) - x <= μ - s ? -Inf : + (μ, σ) = params(d) + x <= μ - σ ? -Inf : x <= μ ? loghalf + 2.0 * log1p(zval(d, x)) : - x < μ + s ? log1p(-0.5 * (1.0 - zval(d, x))^2) : 0.0 + x < μ + σ ? log1p(-0.5 * (1.0 - zval(d, x))^2) : 0.0 end function logccdf(d::SymTriangularDist, x::Float64) - (μ, s) = params(d) - x <= μ - s ? 0.0 : + (μ, σ) = params(d) + x <= μ - σ ? 0.0 : x <= μ ? log1p(-0.5 * (1.0 + zval(d, x))^2) : - x < μ + s ? loghalf + 2.0 * log1p(-zval(d, x)) : -Inf + x < μ + σ ? loghalf + 2.0 * log1p(-zval(d, x)) : -Inf end -quantile(d::SymTriangularDist, p::Float64) = p < 0.5 ? xval(d, sqrt(2.0 * p) - 1.0) : +quantile(d::SymTriangularDist, p::Float64) = p < 0.5 ? xval(d, sqrt(2.0 * p) - 1.0) : xval(d, 1.0 - sqrt(2.0 * (1.0 - p))) cquantile(d::SymTriangularDist, p::Float64) = p > 0.5 ? xval(d, sqrt(2.0 * (1.0-p)) - 1.0) : - xval(d, 1.0 - sqrt(2.0 * p)) + xval(d, 1.0 - sqrt(2.0 * p)) invlogcdf(d::SymTriangularDist, lp::Float64) = lp < loghalf ? xval(d, expm1(0.5*(lp - loghalf))) : xval(d, 1.0 - sqrt(-2.0 * expm1(lp))) @@ -86,16 +85,16 @@ invlogccdf(d::SymTriangularDist, lp::Float64) = lp > loghalf ? xval(d, sqrt(-2.0 xval(d, -(expm1(0.5 * (lp - loghalf)))) -function mgf(d::SymTriangularDist, t::Real) - (μ, s) = params(d) - a = s * t +function mgf(d::SymTriangularDist, t::Float64) + (μ, σ) = params(d) + a = σ * t a == zero(a) && return one(a) 4.0 * exp(μ * t) * (sinh(0.5 * a) / a)^2 end -function cf(d::SymTriangularDist, t::Real) - (μ, s) = params(d) - a = s * t +function cf(d::SymTriangularDist, t::Float64) + (μ, σ) = params(d) + a = σ * t a == zero(a) && return complex(one(a)) 4.0 * cis(μ * t) * (sin(0.5 * a) / a)^2 end @@ -104,5 +103,3 @@ end #### Sampling rand(d::SymTriangularDist) = xval(d, rand() - rand()) - - diff --git a/src/univariate/continuous/tdist.jl b/src/univariate/continuous/tdist.jl index aa3c01ab5..e3b29efac 100644 --- a/src/univariate/continuous/tdist.jl +++ b/src/univariate/continuous/tdist.jl @@ -1,9 +1,7 @@ immutable TDist <: ContinuousUnivariateDistribution - df::Float64 - function TDist(d::Real) - d > zero(d) || error("TDist: df must be positive") - @compat new(Float64(d)) - end + ν::Float64 + + TDist(ν::Real) = (@check_args(TDist, ν > zero(ν)); new(ν)) end @distr_support TDist -Inf Inf @@ -11,49 +9,49 @@ end #### Parameters -dof(d::TDist) = d.df -params(d::TDist) = (d.df,) +dof(d::TDist) = d.ν +params(d::TDist) = (d.ν,) #### Statistics -mean(d::TDist) = d.df > 1.0 ? 0.0 : NaN +mean(d::TDist) = d.ν > 1.0 ? 0.0 : NaN median(d::TDist) = 0.0 mode(d::TDist) = 0.0 function var(d::TDist) - df = d.df - df > 2.0 ? df / (df - 2.0) : - df > 1.0 ? Inf : NaN + ν = d.ν + ν > 2.0 ? ν / (ν - 2.0) : + ν > 1.0 ? Inf : NaN end -skewness(d::TDist) = d.df > 3.0 ? 0.0 : NaN +skewness(d::TDist) = d.ν > 3.0 ? 0.0 : NaN function kurtosis(d::TDist) - df = d.df - df > 4.0 ? 6.0 / (df - 4.0) : - df > 2.0 ? Inf : NaN + ν = d.ν + ν > 4.0 ? 6.0 / (ν - 4.0) : + ν > 2.0 ? Inf : NaN end function entropy(d::TDist) - hdf = 0.5 * d.df - hdfph = hdf + 0.5 - hdfph * (digamma(hdfph) - digamma(hdf)) + 0.5 * log(d.df) + lbeta(hdf,0.5) + h = 0.5 * d.ν + h1 = h + 0.5 + h1 * (digamma(h1) - digamma(h)) + 0.5 * log(d.ν) + lbeta(h, 0.5) end #### Evaluation & Sampling -@_delegate_statsfuns TDist tdist df +@_delegate_statsfuns TDist tdist ν -rand(d::TDist) = StatsFuns.Rmath.tdistrand(d.df) +rand(d::TDist) = StatsFuns.Rmath.tdistrand(d.ν) function cf(d::TDist, t::Real) t == 0 && return complex(1.0) - h = d.df/2 - q = d.df/4 + h = d.ν * 0.5 + q = d.ν * 0.25 t2 = t*t - complex(2*(q*t2)^q*besselk(h,sqrt(d.df)*abs(t))/gamma(h)) + complex(2*(q*t2)^q*besselk(h,sqrt(d.ν)*abs(t))/gamma(h)) end -gradlogpdf(d::TDist, x::Float64) = -((d.df + 1.0) * x) / (x^2 + d.df) +gradlogpdf(d::TDist, x::Float64) = -((d.ν + 1.0) * x) / (x^2 + d.ν) diff --git a/src/univariate/continuous/triangular.jl b/src/univariate/continuous/triangular.jl index 2dbe620bb..f0c3e081d 100644 --- a/src/univariate/continuous/triangular.jl +++ b/src/univariate/continuous/triangular.jl @@ -4,17 +4,13 @@ immutable TriangularDist <: ContinuousUnivariateDistribution c::Float64 function TriangularDist(a::Real, b::Real, c::Real) - a < b || error("TriangularDist: a < b must be true") - a <= c <= b || error("a <= c <= b must be true") - @compat new(Float64(a), Float64(b), Float64(c)) + @check_args(TriangularDist, a < b) + @check_args(TriangularDist, a <= c <= b) + new(a, b, c) end - function TriangularDist(a::Real, b::Real) - a < b || error("TriangularDist: a < b must be true") - @compat a_ = Float64(a) - @compat b_ = Float64(b) - c_ = middle(a_, b_) - new(a_, b_, c_) + @check_args(TriangularDist, a < b) + new(a, b, middle(a, b)) end end @@ -113,7 +109,6 @@ function rand(d::TriangularDist) (a, b, c) = params(d) b_m_a = b - a u = rand() - b_m_a * u < (c - a) ? d.a + sqrt(u * b_m_a * (c - a)) : + b_m_a * u < (c - a) ? d.a + sqrt(u * b_m_a * (c - a)) : d.b - sqrt((1.0 - u) * b_m_a * (b - c)) end - diff --git a/src/univariate/continuous/triweight.jl b/src/univariate/continuous/triweight.jl index 1c2c688e9..82edbf7fd 100644 --- a/src/univariate/continuous/triweight.jl +++ b/src/univariate/continuous/triweight.jl @@ -1,55 +1,57 @@ immutable Triweight <: ContinuousUnivariateDistribution - location::Float64 - scale::Float64 - function Triweight(l::Real, s::Real) - s > zero(s) || error("scale must be positive") - @compat new(Float64(l), Float64(s)) - end + μ::Float64 + σ::Float64 + + Triweight(μ::Real, σ::Real) = (@check_args(Triweight, σ > zero(σ)); new(μ, σ)) + Triweight(μ::Real) = new(μ, 1.0) + Triweight() = new(0.0, 1.0) end -Triweight(location::Real) = Triweight(location, 1.0) -Triweight() = Triweight(0.0, 1.0) +@distr_support Triweight d.μ - d.σ d.μ + d.σ -@distr_support Triweight d.location-d.scale d.location+d.scale ## Parameters -params(d::Triweight) = (d.location, d.scale) + +location(d::Triweight) = d.μ +scale(d::Triweight) = d.σ +params(d::Triweight) = (d.μ, d.σ) + ## Properties -mean(d::Triweight) = d.location -median(d::Triweight) = d.location -mode(d::Triweight) = d.location +mean(d::Triweight) = d.μ +median(d::Triweight) = d.μ +mode(d::Triweight) = d.μ -var(d::Triweight) = d.scale*d.scale/9.0 +var(d::Triweight) = d.σ^2 / 9.0 skewness(d::Triweight) = 0.0 -kurtosis(d::Triweight) = 1/33-3 +kurtosis(d::Triweight) = -2.9696969696969697 # 1/33-3 ## Functions function pdf(d::Triweight, x::Real) - u = abs(x - d.location)/d.scale - u >= 1 ? 0.0 : 1.09375*(1-u*u)^3/d.scale + u = abs(x - d.μ)/d.σ + u >= 1 ? 0.0 : 1.09375*(1-u*u)^3/d.σ end function cdf(d::Triweight, x::Real) - u = (x - d.location)/d.scale + u = (x - d.μ)/d.σ u <= -1 ? 0.0 : u >= 1 ? 1.0 : 0.03125*(1+u)^4*@horner(u,16.0,-29.0,20.0,-5.0) end + function ccdf(d::Triweight, x::Real) - u = (d.location - x)/d.scale + u = (d.μ - x)/d.σ u <= -1 ? 1.0 : u >= 1 ? 0.0 : 0.03125*(1+u)^4*@horner(u,16.0,-29.0,20.0,-5.0) end @quantile_newton Triweight - -function mgf(d::Triweight, t::Real) - a = d.scale*t +function mgf(d::Triweight, t::Float64) + a = d.σ*t a2 = a*a - a == 0 ? one(a) : 105.0*exp(d.location*t)*((15.0/a2+1.0)*cosh(a)-(15.0/a2-6.0)/a*sinh(a))/(a2*a2) + a == 0 ? one(a) : 105.0*exp(d.μ*t)*((15.0/a2+1.0)*cosh(a)-(15.0/a2-6.0)/a*sinh(a))/(a2*a2) end -function cf(d::Triweight, t::Real) - a = d.scale*t +function cf(d::Triweight, t::Float64) + a = d.σ*t a2 = a*a - a == 0 ? complex(one(a)) : 105.0*cis(d.location*t)*((1.0-15.0/a2)*cos(a)+(15.0/a2-6.0)/a*sin(a))/(a2*a2) + a == 0 ? complex(one(a)) : 105.0*cis(d.μ*t)*((1.0-15.0/a2)*cos(a)+(15.0/a2-6.0)/a*sin(a))/(a2*a2) end diff --git a/src/univariate/continuous/uniform.jl b/src/univariate/continuous/uniform.jl index 9af5007d0..d31a1f6d1 100644 --- a/src/univariate/continuous/uniform.jl +++ b/src/univariate/continuous/uniform.jl @@ -2,11 +2,7 @@ immutable Uniform <: ContinuousUnivariateDistribution a::Float64 b::Float64 - function Uniform(a::Real, b::Real) - a < b || error("Uniform: a must be less than b") - @compat new(Float64(a), Float64(b)) - end - + Uniform(a::Real, b::Real) = (@check_args(Uniform, a < b); new(a, b)) Uniform() = new(0.0, 1.0) end diff --git a/src/univariate/continuous/vonmises.jl b/src/univariate/continuous/vonmises.jl index b963d9824..30f5ef23b 100644 --- a/src/univariate/continuous/vonmises.jl +++ b/src/univariate/continuous/vonmises.jl @@ -9,12 +9,11 @@ immutable VonMises <: ContinuousUnivariateDistribution I0κ::Float64 # I0(κ), where I0 is the modified Bessel function of order 0 function VonMises(μ::Real, κ::Real) - κ > zero(κ) || error("kappa must be positive") - @compat new(Float64(μ), Float64(κ), besseli(0, κ)) + @check_args(VonMises, κ > zero(κ)) + new(μ, κ, besseli(0, κ)) end - - @compat VonMises(κ::Real) = VonMises(0.0, Float64(κ)) - VonMises() = VonMises(0.0, 1.0) + VonMises(κ::Real) = VonMises(0.0, κ) + VonMises() = new(0.0, 1.0) end show(io::IO, d::VonMises) = show(io, d, (:μ, :κ)) diff --git a/src/univariate/continuous/weibull.jl b/src/univariate/continuous/weibull.jl index f146e9a9b..739525602 100644 --- a/src/univariate/continuous/weibull.jl +++ b/src/univariate/continuous/weibull.jl @@ -1,12 +1,11 @@ immutable Weibull <: ContinuousUnivariateDistribution α::Float64 # shape - β::Float64 # scale + θ::Float64 # scale - function Weibull(α::Real, β::Real) - (zero(α) < α && zero(β) < β) || error("Weibull: both shape and scale must be positive") - @compat new(Float64(α), Float64(β)) + function Weibull(α::Real, θ::Real) + @check_args(Weibull, α > zero(α) && θ > zero(θ)) + new(α, θ) end - Weibull(α::Real) = Weibull(α, 1.0) Weibull() = new(1.0, 1.0) end @@ -17,41 +16,41 @@ end #### Parameters shape(d::Weibull) = d.α -scale(d::Weibull) = d.β +scale(d::Weibull) = d.θ -params(d::Weibull) = (d.α, d.β) +params(d::Weibull) = (d.α, d.θ) #### Statistics -mean(d::Weibull) = d.β * gamma(1.0 + 1.0 / d.α) -median(d::Weibull) = d.β * logtwo ^ (1.0 / d.α) -mode(d::Weibull) = d.α > 1.0 ? (iα = 1.0 / d.α; d.β * (1.0 - iα) ^ iα) : 0.0 +mean(d::Weibull) = d.θ * gamma(1.0 + 1.0 / d.α) +median(d::Weibull) = d.θ * logtwo ^ (1.0 / d.α) +mode(d::Weibull) = d.α > 1.0 ? (iα = 1.0 / d.α; d.θ * (1.0 - iα) ^ iα) : 0.0 -var(d::Weibull) = d.β^2 * gamma(1.0 + 2.0 / d.α) - mean(d)^2 +var(d::Weibull) = d.θ^2 * gamma(1.0 + 2.0 / d.α) - mean(d)^2 function skewness(d::Weibull) μ = mean(d) σ2 = var(d) σ = sqrt(σ2) r = μ / σ - gamma(1.0 + 3.0 / d.α) * (d.β / σ)^3 - 3.0 * r - r^3 + gamma(1.0 + 3.0 / d.α) * (d.θ / σ)^3 - 3.0 * r - r^3 end function kurtosis(d::Weibull) - α, β = params(d) + α, θ = params(d) μ = mean(d) σ = std(d) γ = skewness(d) r = μ / σ r2 = r^2 r4 = r2^2 - (β / σ)^4 * gamma(1.0 + 4.0 / α) - 4.0 * γ * r - 6.0 * r2 - r4 - 3.0 + (θ / σ)^4 * gamma(1.0 + 4.0 / α) - 4.0 * γ * r - 6.0 * r2 - r4 - 3.0 end function entropy(d::Weibull) - α, β = params(d) - 0.5772156649015328606 * (1.0 - 1.0 / α) + log(β / α) + 1.0 + α, θ = params(d) + 0.5772156649015328606 * (1.0 - 1.0 / α) + log(θ / α) + 1.0 end @@ -59,9 +58,9 @@ end function pdf(d::Weibull, x::Float64) if x >= 0.0 - α, β = params(d) - z = x / β - (α / β) * z^(α - 1.0) * exp(-z^α) + α, θ = params(d) + z = x / θ + (α / θ) * z^(α - 1.0) * exp(-z^α) else 0.0 end @@ -69,16 +68,16 @@ end function logpdf(d::Weibull, x::Float64) if x >= 0.0 - α, β = params(d) - z = x / β - log(α / β) + (α - 1.0) * log(z) - z^α + α, θ = params(d) + z = x / θ + log(α / θ) + (α - 1.0) * log(z) - z^α else -Inf end end -zv(d::Weibull, x::Float64) = (x / d.β) ^ d.α -xv(d::Weibull, z::Float64) = d.β * z ^ (1.0 / d.α) +zv(d::Weibull, x::Float64) = (x / d.θ) ^ d.α +xv(d::Weibull, z::Float64) = d.θ * z ^ (1.0 / d.α) cdf(d::Weibull, x::Float64) = x > 0.0 ? -expm1(-zv(d, x)) : 0.0 ccdf(d::Weibull, x::Float64) = x > 0.0 ? exp(-zv(d, x)) : 1.0 @@ -92,9 +91,9 @@ invlogccdf(d::Weibull, lp::Float64) = xv(d, -lp) function gradlogpdf(d::Weibull, x::Float64) if insupport(Weibull, x) - α, β = params(d) - (α - 1.0) / x - α * x^(α - 1.0) / (β^α) - else + α, θ = params(d) + (α - 1.0) / x - α * x^(α - 1.0) / (θ^α) + else 0.0 end end @@ -103,5 +102,3 @@ end #### Sampling rand(d::Weibull) = xv(d, randexp()) - - diff --git a/src/univariate/discrete/bernoulli.jl b/src/univariate/discrete/bernoulli.jl index ed89b2e35..e784452f8 100644 --- a/src/univariate/discrete/bernoulli.jl +++ b/src/univariate/discrete/bernoulli.jl @@ -1,17 +1,10 @@ -# Bernoulli distribution - - -#### Type and Constructor - immutable Bernoulli <: DiscreteUnivariateDistribution p::Float64 - function Bernoulli(p::Float64) - 0.0 <= p <= 1.0 || error("p must be in [0, 1].") + function Bernoulli(p::Real) + @check_args(Bernoulli, zero(p) <= p <= one(p)) new(p) end - - @compat Bernoulli(p::Real) = Bernoulli(Float64(p)) Bernoulli() = new(0.5) end @@ -48,7 +41,7 @@ function median(d::Bernoulli) p > 0.5 ? 1.0 : 0.5 end -function entropy(d::Bernoulli) +function entropy(d::Bernoulli) p0 = failprob(d) p1 = succprob(d) (p0 == 0.0 || p0 == 1.0) ? 0.0 : -(p0 * log(p0) + p1 * log(p1)) @@ -57,7 +50,7 @@ end #### Evaluation pdf(d::Bernoulli, x::Bool) = x ? succprob(d) : failprob(d) -pdf(d::Bernoulli, x::Int) = x == 0 ? failprob(d) : +pdf(d::Bernoulli, x::Int) = x == 0 ? failprob(d) : x == 1 ? succprob(d) : 0.0 pdf(d::Bernoulli) = Float64[failprob(d), succprob(d)] @@ -126,5 +119,3 @@ function suffstats{T<:Integer}(::Type{Bernoulli}, x::AbstractArray{T}, w::Abstra end BernoulliStats(c0, c1) end - - diff --git a/src/univariate/discrete/binomial.jl b/src/univariate/discrete/binomial.jl index 2ed142c52..0f2f808d2 100644 --- a/src/univariate/discrete/binomial.jl +++ b/src/univariate/discrete/binomial.jl @@ -1,21 +1,22 @@ - immutable Binomial <: DiscreteUnivariateDistribution n::Int p::Float64 - function Binomial(n::Int, p::Float64) - n >= 0 || error("n must be non-negative but is $n.") - 0.0 <= p <= 1.0 || error("p must be in [0, 1] but is $p") + function Binomial(n::Real, p::Real) + @check_args(Binomial, n >= zero(n)) + @check_args(Binomial, zero(p) <= p <= one(p)) new(n, p) end - - @compat Binomial(n::Integer, p::Real) = Binomial(round(Int, n), Float64(p)) - @compat Binomial(n::Integer) = Binomial(round(Int, n), 0.5) + function Binomial(n::Real) + @check_args(Binomial, n >= zero(n)) + new(n, 0.5) + end Binomial() = new(1, 0.5) end @distr_support Binomial 0 d.n + #### Parameters ntrials(d::Binomial) = d.n @@ -127,9 +128,7 @@ immutable BinomialStats <: SufficientStats ne::Float64 # the number of experiments n::Int # the number of trials in each experiment - function BinomialStats(ns::Real, ne::Real, n::Integer) - @compat new(Float64(ns), Float64(ne), round(Int, n)) - end + BinomialStats(ns::Real, ne::Real, n::Integer) = new(ns, ne, n) end function suffstats{T<:Integer}(::Type{Binomial}, n::Integer, x::AbstractArray{T}) diff --git a/src/univariate/discrete/categorical.jl b/src/univariate/discrete/categorical.jl index 6731aac28..4a52163da 100644 --- a/src/univariate/discrete/categorical.jl +++ b/src/univariate/discrete/categorical.jl @@ -5,12 +5,12 @@ immutable Categorical <: DiscreteUnivariateDistribution Categorical(p::Vector{Float64}, ::NoArgCheck) = new(length(p), p) function Categorical(p::Vector{Float64}) - isprobvec(p) || error("p is not a valid probability vector.") + @check_args(Categorical, isprobvec(p)) new(length(p), p) end function Categorical(k::Integer) - k >= 1 || error("k must be a positive integer.") + @check_args(Categorical, k >= 1) new(k, fill(1.0/k, k)) end end diff --git a/src/univariate/discrete/discreteuniform.jl b/src/univariate/discrete/discreteuniform.jl index 110c7e0c4..45aabd3fb 100644 --- a/src/univariate/discrete/discreteuniform.jl +++ b/src/univariate/discrete/discreteuniform.jl @@ -3,13 +3,11 @@ immutable DiscreteUniform <: DiscreteUnivariateDistribution b::Int pv::Float64 - function DiscreteUniform(a::Int, b::Int) - a <= b || error("a and b must satisfy a <= b") + function DiscreteUniform(a::Real, b::Real) + @check_args(DiscreteUniform, a <= b) new(a, b, 1.0 / (b - a + 1)) end - - @compat DiscreteUniform(a::Real, b::Real) = DiscreteUniform(round(Int, a), round(Int, b)) - DiscreteUniform(b::Real) = DiscreteUniform(0, round(Int, b)) + DiscreteUniform(b::Real) = DiscreteUniform(0, b) DiscreteUniform() = new(0, 1, 0.5) end @@ -33,16 +31,16 @@ mean(d::DiscreteUniform) = middle(d.a, d.b) median(d::DiscreteUniform) = middle(d.a, d.b) -@compat var(d::DiscreteUniform) = (abs2(Float64(span(d))) - 1.0) / 12.0 +var(d::DiscreteUniform) = (span(d)^2 - 1.0) / 12.0 skewness(d::DiscreteUniform) = 0.0 function kurtosis(d::DiscreteUniform) - @compat n2 = abs2(Float64(span(d))) - return -1.2 * (n2 + 1.0) / (n2 - 1.0) + n2 = span(d)^2 + -1.2 * (n2 + 1.0) / (n2 - 1.0) end -@compat entropy(d::DiscreteUniform) = log(Float64(span(d))) +entropy(d::DiscreteUniform) = log(span(d)) mode(d::DiscreteUniform) = d.a modes(d::DiscreteUniform) = [d.a:d.b] diff --git a/src/univariate/discrete/geometric.jl b/src/univariate/discrete/geometric.jl index d77901ac7..728949f0b 100644 --- a/src/univariate/discrete/geometric.jl +++ b/src/univariate/discrete/geometric.jl @@ -2,11 +2,10 @@ immutable Geometric <: DiscreteUnivariateDistribution p::Float64 function Geometric(p::Real) - zero(p) < p < one(p) || error("prob must be in (0, 1)") - @compat new(Float64(p)) + @check_args(Geometric, zero(p) < p < one(p)) + new(p) end - - Geometric() = Geometric(0.5) # Flips of a fair coin + Geometric() = new(0.5) end @distr_support Geometric 0 Inf @@ -116,7 +115,7 @@ immutable GeometricStats <: SufficientStats sx::Float64 tw::Float64 - @compat GeometricStats(sx::Real, tw::Real) = new(Float64(sx), Float64(tw)) + GeometricStats(sx::Real, tw::Real) = new(sx, tw) end suffstats{T<:Integer}(::Type{Geometric}, x::AbstractArray{T}) = GeometricStats(sum(x), length(x)) diff --git a/src/univariate/discrete/hypergeometric.jl b/src/univariate/discrete/hypergeometric.jl index 8845eb54f..5c58c89aa 100644 --- a/src/univariate/discrete/hypergeometric.jl +++ b/src/univariate/discrete/hypergeometric.jl @@ -11,14 +11,11 @@ immutable Hypergeometric <: DiscreteUnivariateDistribution nf::Int # number of failures in population n::Int # sample size - function Hypergeometric(ns::Int, nf::Int, n::Int) - ns >= 0 || error("ns must be non-negative.") - nf >= 0 || error("nf must be non-negative.") - 0 < n < ns + nf || error("n must have 0 < n < ns + nf") + function Hypergeometric(ns::Real, nf::Real, n::Real) + @check_args(Hypergeometric, ns >= zero(ns) && nf >= zero(nf)) + @check_args(Hypergeometric, zero(n) < n < ns + nf) new(ns, nf, n) end - - @compat Hypergeometric(ns::Real, nf::Real, n::Real) = Hypergeometric(round(Int, ns), round(Int, nf), round(Int, n)) end @distr_support Hypergeometric max(d.n - d.nf, 0) min(d.ns, d.n) diff --git a/src/univariate/discrete/negativebinomial.jl b/src/univariate/discrete/negativebinomial.jl index 5cbde211a..4990662e4 100644 --- a/src/univariate/discrete/negativebinomial.jl +++ b/src/univariate/discrete/negativebinomial.jl @@ -8,15 +8,13 @@ immutable NegativeBinomial <: DiscreteUnivariateDistribution r::Int p::Float64 - function NegativeBinomial(r::Int, p::Float64) - r > 0 || error("r must be positive.") - 0.0 < p <= 1.0 || error("prob must be in (0, 1].") + function NegativeBinomial(r::Real, p::Real) + @check_args(NegativeBinomial, r > zero(r)) + @check_args(NegativeBinomial, zero(p) < p <= one(p)) new(r, p) end - - @compat NegativeBinomial(r::Real, p::Real) = NegativeBinomial(round(Int, r), Float64(p)) NegativeBinomial(r::Real) = NegativeBinomial(r, 0.5) - NegativeBinomial() = new(1.0, 0.5) + NegativeBinomial() = new(1, 0.5) end @distr_support NegativeBinomial 0 Inf diff --git a/src/univariate/discrete/noncentralhypergeometric.jl b/src/univariate/discrete/noncentralhypergeometric.jl index c04cdf863..bfee2f00e 100644 --- a/src/univariate/discrete/noncentralhypergeometric.jl +++ b/src/univariate/discrete/noncentralhypergeometric.jl @@ -1,4 +1,5 @@ # Noncentral hypergeometric distribution +# TODO: this distribution needs clean-up and testing abstract NoncentralHypergeometric <: DiscreteUnivariateDistribution @@ -38,13 +39,12 @@ immutable FisherNoncentralHypergeometric <: NoncentralHypergeometric nf::Int # number of failures in population n::Int # sample size ω::Float64 # odds ratio - function FisherNoncentralHypergeometric(s::Real, f::Real, n::Real, ω::Float64) - isinteger(s) && zero(s) <= s || error("ns must be a non-negative integer") - isinteger(f) && zero(f) <= f || error("nf must be a non-negative integer") - isinteger(n) && zero(n) < n < s + f || - error("n must be a positive integer <= (ns + nf)") - zero(ω) < ω || error("ω must be a positive real value") - @compat new(Float64(s), Float64(f), Float64(n), Float64(ω)) + + function FisherNoncentralHypergeometric(ns::Real, nf::Real, n::Real, ω::Float64) + @check_args(FisherNoncentralHypergeometric, ns >= zero(ns) && nf >= zero(nf)) + @check_args(FisherNoncentralHypergeometric, zero(n) < n < ns + nf) + @check_args(FisherNoncentralHypergeometric, ω > zero(ω)) + new(ns, nf, n, ω) end end @@ -81,13 +81,11 @@ immutable WalleniusNoncentralHypergeometric <: NoncentralHypergeometric nf::Int # number of failures in population n::Int # sample size ω::Float64 # odds ratio - function WalleniusNoncentralHypergeometric(s::Real, f::Real, n::Real, ω::Float64) - isinteger(s) && zero(s) <= s || error("ns must be a non-negative integer") - isinteger(f) && zero(f) <= f || error("nf must be a non-negative integer") - isinteger(n) && zero(n) < n < s + f || - error("n must be a positive integer <= (ns + nf)") - zero(ω) < ω || error("ω must be a positive real value") - @compat new(Float64(s), Float64(f), Float64(n), Float64(ω)) + function WalleniusNoncentralHypergeometric(ns::Real, nf::Real, n::Real, ω::Float64) + @check_args(WalleniusNoncentralHypergeometric, ns >= zero(ns) && nf >= zero(nf)) + @check_args(WalleniusNoncentralHypergeometric, zero(n) < n < ns + nf) + @check_args(WalleniusNoncentralHypergeometric, ω > zero(ω)) + new(ns, nf, n, ω) end end @@ -104,6 +102,3 @@ function pdf(d::WalleniusNoncentralHypergeometric, k::Int) end logpdf(d::WalleniusNoncentralHypergeometric, k::Int) = log(pdf(d, k)) - - - diff --git a/src/univariate/discrete/poisson.jl b/src/univariate/discrete/poisson.jl index c0edc5867..f6f3bb886 100644 --- a/src/univariate/discrete/poisson.jl +++ b/src/univariate/discrete/poisson.jl @@ -1,12 +1,7 @@ immutable Poisson <: DiscreteUnivariateDistribution λ::Float64 - function Poisson(λ::Float64) - λ > 0.0 || error("λ must be positive.") - new(λ) - end - - @compat Poisson(λ::Real) = Poisson(Float64(λ)) + Poisson(λ::Real) = (@check_args(Poisson, λ > zero(λ)); new(λ)) Poisson() = new(1.0) end @@ -88,7 +83,7 @@ immutable PoissonStats <: SufficientStats tw::Float64 # total sample weight end -@compat suffstats{T<:Integer}(::Type{Poisson}, x::AbstractArray{T}) = PoissonStats(Float64(sum(x)), Float64(length(x))) +suffstats{T<:Integer}(::Type{Poisson}, x::AbstractArray{T}) = PoissonStats(sum(x), length(x)) function suffstats{T<:Integer}(::Type{Poisson}, x::AbstractArray{T}, w::AbstractArray{Float64}) n = length(x) diff --git a/src/univariate/discrete/poissonbinomial.jl b/src/univariate/discrete/poissonbinomial.jl index 1c8e7d73d..afa19394f 100644 --- a/src/univariate/discrete/poissonbinomial.jl +++ b/src/univariate/discrete/poissonbinomial.jl @@ -1,3 +1,4 @@ +# TODO: this distribution may need clean-up # Computes the pdf of a poisson-binomial random variable using # fast fourier transform diff --git a/src/univariate/discrete/skellam.jl b/src/univariate/discrete/skellam.jl index d5fdd3a5c..afb744721 100644 --- a/src/univariate/discrete/skellam.jl +++ b/src/univariate/discrete/skellam.jl @@ -2,21 +2,16 @@ immutable Skellam <: DiscreteUnivariateDistribution μ1::Float64 μ2::Float64 - function Skellam(μ1::Float64, μ2::Float64) - μ1 > 0.0 && μ2 > 0.0 || error("μ1 and μ2 must be positive.") + function Skellam(μ1::Real, μ2::Real) + @check_args(Skellam, μ1 > zero(μ1) && μ2 > zero(μ2)) new(μ1, μ2) end - - @compat Skellam(μ1::Real, μ2::Real) = Skellam(Float64(μ1), Float64(μ2)) - Skellam(μ::Real) = Skellam(μ, μ) - Skellam() = new(1.0, 1.0) end @distr_support Skellam -Inf Inf - ### Parameters params(d::Skellam) = (d.μ1, d.μ2) @@ -55,4 +50,3 @@ end ### Sampling rand(d::Skellam) = rand(Poisson(d.μ1)) - rand(Poisson(d.μ2)) - diff --git a/src/utils.jl b/src/utils.jl index c29a33906..3c24d778c 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,3 +1,13 @@ +## macro for argument checking + +macro check_args(D, cond) + quote + if !($cond) + throw(ArgumentError(string( + $(string(D)), ": the condition ", $(string(cond)), " is not satisfied."))) + end + end +end ## a type to indicate zero vector @@ -117,4 +127,3 @@ function trycholfact(a::Matrix{Float64}) return e end end -