Use Number instead of Real for mgf and cf

src/Distributions.jl

@@ -1,7 +1,7 @@
 module Distributions
 
 using StatsBase, PDMats, StatsFuns, Statistics
-using StatsFuns: logtwo, invsqrt2, invsqrt2π
+using StatsFuns: logtwo, invsqrt2, invsqrt2π, sqrt2
 
 import QuadGK: quadgk
 import Base: size, length, convert, show, getindex, rand, vec, inv

src/univariate/continuous/biweight.jl

@@ -54,16 +54,18 @@ end
 
 @quantile_newton Biweight
 
-function mgf(d::Biweight{T}, t::Real) where T<:Real
-    a = d.σ*t
+function mgf(d::Biweight, t::Number)
+    a = d.σ * t
+    abs(real(a)) < 1 || throw(DomainError(t, "|σ ⋅ real(t)| has to be smaller than 1"))
     a2 = a^2
-    a == 0 ? one(T) :
-    15exp(d.μ * t) * (-3cosh(a) + (a + 3/a) * sinh(a)) / (a2^2)
+    result = 15exp(d.μ * t) * (-3cosh(a) + (a + 3/a) * sinh(a)) / (a2^2)
+    return iszero(a) ? one(result) : result
+
 end
 
-function cf(d::Biweight{T}, t::Real) where T<:Real
+function cf(d::Biweight, t::Number)
     a = d.σ * t
     a2 = a^2
-    a == 0 ? one(T)+zero(T)*im :
-    -15cis(d.μ * t) * (3cos(a) + (a - 3/a) * sin(a)) / (a2^2)
+    result = -15cis(d.μ * t) * (3cos(a) + (a - 3/a) * sin(a)) / (a2^2)
+    return iszero(a) ? complex(one(result)) : result
 end

src/univariate/continuous/cauchy.jl

@@ -97,8 +97,11 @@ function cquantile(d::Cauchy, p::Real)
     μ + σ * tan(π * (1//2 - p))
 end
 
-mgf(d::Cauchy{T}, t::Real) where {T<:Real} = t == zero(t) ? one(T) : T(NaN)
-cf(d::Cauchy, t::Real) = exp(im * (t * d.μ) - d.σ * abs(t))
+function mgf(d::Cauchy, t::Number)
+    T = float(partype(d))
+    return iszero(t) ? one(T) : T(NaN)
+end
+cf(d::Cauchy, t::Number) = exp(im * (t * d.μ) - d.σ * abs(t))
 
 #### Affine transformations
 

src/univariate/continuous/chisq.jl

@@ -75,9 +75,12 @@ end
 
 @_delegate_statsfuns Chisq chisq ν
 
-mgf(d::Chisq, t::Real) = (1 - 2 * t)^(-d.ν/2)
+function mgf(d::Chisq, t::Number)
+    real(t) < 0.5 || throw(DomainError(t, "the real part of t has to be smaller than 0.5"))
+    return (1 - 2 * t)^(-d.ν/2)
+end
 
-cf(d::Chisq, t::Real) = (1 - 2 * im * t)^(-d.ν/2)
+cf(d::Chisq, t::Number) = (1 - 2 * im * t)^(-d.ν/2)
 
 gradlogpdf(d::Chisq{T}, x::Real) where {T<:Real} =  x > 0 ? (d.ν/2 - 1) / x - 1//2 : zero(T)
 

src/univariate/continuous/epanechnikov.jl

@@ -65,14 +65,15 @@ end
 
 @quantile_newton Epanechnikov
 
-function mgf(d::Epanechnikov{T}, t::Real) where T<:Real
+function mgf(d::Epanechnikov{T}, t::Number) where T<:Real
     a = d.σ * t
-    a == 0 ? one(T) :
-    3exp(d.μ * t) * (cosh(a) - sinh(a) / a) / a^2
+    abs(real(a)) < 1 || throw(DomainError(t, "|σ ⋅ real(t)| has to be smaller than 1"))
+    result = 3exp(d.μ * t) * (cosh(a) - sinh(a) / a) / a^2
+    return iszero(a) ? one(result) : result
 end
 
-function cf(d::Epanechnikov{T}, t::Real) where T<:Real
+function cf(d::Epanechnikov{T}, t::Number) where T<:Real
     a = d.σ * t
-    a == 0 ? one(T)+zero(T)*im :
-    -3exp(im * d.μ * t) * (cos(a) - sin(a) / a) / a^2
+    result = -3exp(im * d.μ * t) * (cos(a) - sin(a) / a) / a^2
+    return iszero(a) ? complex(one(result)) : result
 end

src/univariate/continuous/erlang.jl

@@ -75,8 +75,11 @@ function entropy(d::Erlang)
     α + loggamma(α) + (1 - α) * digamma(α) + log(θ)
 end
 
-mgf(d::Erlang, t::Real) = (1 - t * d.θ)^(-d.α)
-cf(d::Erlang, t::Real)  = (1 - im * t * d.θ)^(-d.α)
+function mgf(d::Erlang, t::Number)
+    real(t) < inv(d.θ) || throw(DomainError(t, "the real part of t should be smaller than θ⁻¹")) 
+    return (1 - t * d.θ)^(-d.α)
+end
+cf(d::Erlang, t::Number)  = (1 - im * t * d.θ)^(-d.α)
 
 
 #### Evaluation & Sampling

src/univariate/continuous/exponential.jl

@@ -93,8 +93,11 @@ invlogccdf(d::Exponential, lp::Real) = -xval(d, lp)
 
 gradlogpdf(d::Exponential{T}, x::Real) where {T<:Real} = x > 0 ? -rate(d) : zero(T)
 
-mgf(d::Exponential, t::Real) = 1/(1 - t * scale(d))
-cf(d::Exponential, t::Real) = 1/(1 - t * im * scale(d))
+function mgf(d::Exponential, t::Number)
+    real(t) < rate(d) || throw(DomainError(t, "the real part of t should be smaller than λ"))
+    return 1/(1 - t * scale(d))
+end
+cf(d::Exponential, t::Number) = 1/(1 - t * im * scale(d))
 
 
 #### Sampling

src/univariate/continuous/gamma.jl

@@ -75,9 +75,12 @@ function entropy(d::Gamma)
     α + loggamma(α) + (1 - α) * digamma(α) + log(θ)
 end
 
-mgf(d::Gamma, t::Real) = (1 - t * d.θ)^(-d.α)
+function mgf(d::Gamma, t::Number)
+    real(t) < rate(d) || throw(DomainError(t, "the real part of t should smaller than θ⁻¹"))
+    return (1 - t * d.θ)^(-d.α)
+end
 
-cf(d::Gamma, t::Real) = (1 - im * t * d.θ)^(-d.α)
+cf(d::Gamma, t::Number) = (1 - im * t * d.θ)^(-d.α)
 
 function kldivergence(p::Gamma, q::Gamma)
     # We use the parametrization with the scale θ

src/univariate/continuous/inversegamma.jl

@@ -108,14 +108,18 @@ cquantile(d::InverseGamma, p::Real) = inv(quantile(d.invd, p))
 invlogcdf(d::InverseGamma, p::Real) = inv(invlogccdf(d.invd, p))
 invlogccdf(d::InverseGamma, p::Real) = inv(invlogcdf(d.invd, p))
 
-function mgf(d::InverseGamma{T}, t::Real) where T<:Real
+function mgf(d::InverseGamma, t::Number)
     (a, b) = params(d)
-    t == zero(t) ? one(T) : 2(-b*t)^(0.5a) / gamma(a) * besselk(a, sqrt(-4*b*t))
+    # Where does this formula comes from?
+    # Wikipedia says (wrongly I think) that it is not defined
+    result = 2(-b*t)^(0.5a) / gamma(a) * besselk(a, sqrt(-4*b*t))
+    return iszero(t) ? one(result) : result
 end
 
-function cf(d::InverseGamma{T}, t::Real) where T<:Real
+function cf(d::InverseGamma, t::Number)
     (a, b) = params(d)
-    t == zero(t) ? one(T)+zero(T)*im : 2(-im*b*t)^(0.5a) / gamma(a) * besselk(a, sqrt(-4*im*b*t))
+    result = 2(-im*b*t)^(0.5a) / gamma(a) * besselk(a, sqrt(-4*im*b*t))
+    iszero(t) ? complex(one(result)) : result
 end
 
 

src/univariate/continuous/laplace.jl

@@ -95,11 +95,12 @@ function gradlogpdf(d::Laplace, x::Real)
     x > μ ? -g : g
 end
 
-function mgf(d::Laplace, t::Real)
+function mgf(d::Laplace, t::Number)
     st = d.θ * t
+    abs(real(t)) < 1 / d.θ || throw(DomainError(t, "the absolute value of the real part of t should be smaller than θ⁻¹"))
     exp(t * d.μ) / ((1 - st) * (1 + st))
 end
-function cf(d::Laplace, t::Real)
+function cf(d::Laplace, t::Number)
     st = d.θ * t
     cis(t * d.μ) / (1+st*st)
 end

src/univariate/continuous/levy.jl

@@ -90,9 +90,12 @@ ccdf(d::Levy{T}, x::Real) where {T<:Real} =  x <= d.μ ? one(T) : erf(sqrt(d.σ
 quantile(d::Levy, p::Real) = d.μ + d.σ / (2*erfcinv(p)^2)
 cquantile(d::Levy, p::Real) = d.μ + d.σ / (2*erfinv(p)^2)
 
-mgf(d::Levy{T}, t::Real) where {T<:Real} = t == zero(t) ? one(T) : T(NaN)
+function mgf(d::Levy, t::Number)
+    T = float(partype(d))
+    iszero(t) ? one(T) : T(NaN)
+end
 
-function cf(d::Levy, t::Real)
+function cf(d::Levy, t::Number)
     μ, σ = params(d)
     exp(im * μ * t - sqrt(-2im * σ * t))
 end

src/univariate/continuous/logistic.jl

@@ -95,9 +95,12 @@ function gradlogpdf(d::Logistic, x::Real)
     ((2e) / (1 + e) - 1) / d.θ
 end
 
-mgf(d::Logistic, t::Real) = exp(t * d.μ) / sinc(d.θ * t)
+function mgf(d::Logistic, t::Number)
+    abs(real(t)) < inv(d.θ) || throw(DomainError(t, "the absolute value of the real part of t should be smaller than θ⁻¹"))
+    return exp(t * d.μ) / sinc(d.θ * t)
+end
 
-function cf(d::Logistic, t::Real)
+function cf(d::Logistic, t::Number)
     a = (π * t) * d.θ
-    a == zero(a) ? complex(one(a)) : cis(t * d.μ) * (a / sinh(a))
+    return iszero(a) ? complex(one(a)) : cis(t * d.μ) * (a / sinh(a))
 end

src/univariate/continuous/noncentralchisq.jl

@@ -61,11 +61,12 @@ var(d::NoncentralChisq) = 2(d.ν + 2d.λ)
 skewness(d::NoncentralChisq) = 2sqrt2*(d.ν + 3d.λ)/sqrt(d.ν + 2d.λ)^3
 kurtosis(d::NoncentralChisq) = 12(d.ν + 4d.λ)/(d.ν + 2d.λ)^2
 
-function mgf(d::NoncentralChisq, t::Real)
-    exp(d.λ * t/(1 - 2t))*(1 - 2t)^(-d.ν/2)
+function mgf(d::NoncentralChisq, t::Number)
+    real(t) < 0.5 || throw(DomainError(t, "the real part of t should be smaller than 1/2"))
+    return exp(d.λ * t/(1 - 2t))*(1 - 2t)^(-d.ν/2)
 end
 
-function cf(d::NoncentralChisq, t::Real)
+function cf(d::NoncentralChisq, t::Number)
     cis(d.λ * t/(1 - 2im*t))*(1 - 2im*t)^(-d.ν/2)
 end
 

src/univariate/continuous/normal.jl

@@ -90,8 +90,8 @@ end
 
 gradlogpdf(d::Normal, x::Real) = (d.μ - x) / d.σ^2
 
-mgf(d::Normal, t::Real) = exp(t * d.μ + d.σ^2 / 2 * t^2)
-cf(d::Normal, t::Real) = exp(im * t * d.μ - d.σ^2 / 2 * t^2)
+mgf(d::Normal, t::Number) = exp(t * d.μ + d.σ^2 / 2 * t^2)
+cf(d::Normal, t::Number) = exp(im * t * d.μ - d.σ^2 / 2 * t^2)
 
 #### Affine transformations
 

src/univariate/continuous/skewnormal.jl

@@ -54,9 +54,9 @@ logpdf(d::SkewNormal, x::Real) = log(2) - log(d.ω) + normlogpdf((x-d.ξ) / d.ω
 #cdf requires Owen's T function.
 #cdf/quantile etc 
 
-mgf(d::SkewNormal, t::Real) = 2 * exp(d.ξ * t + (d.ω^2 * t^2)/2 ) * normcdf(d.ω * delta(d) * t)
+mgf(d::SkewNormal, t::Number) = 2 * exp(d.ξ * t + (d.ω^2 * t^2)/2) * normcdf(d.ω * delta(d) * t)
 
-cf(d::SkewNormal, t::Real) = exp(im * t * d.ξ - (d.ω^2 * t^2)/2) * (1 + im * erfi((d.ω * delta(d) * t)/(sqrt(2))) )
+cf(d::SkewNormal, t::Number) = exp(im * t * d.ξ - (d.ω^2 * t^2)/2) * (1 + im * erfi((d.ω * delta(d) * t)/sqrt2))
 
 #### Sampling
 function rand(rng::AbstractRNG, d::SkewNormal)

src/univariate/continuous/tdist.jl

@@ -80,12 +80,12 @@ end
 
 rand(rng::AbstractRNG, d::TDist) = randn(rng) / ( isinf(d.ν) ? 1 : sqrt(rand(rng, Chisq(d.ν))/d.ν) )
 
-function cf(d::TDist{T}, t::Real) where T <: Real
+function cf(d::TDist{T}, t::Real) where {T<:Real}
     isinf(d.ν) && return cf(Normal(zero(T), one(T)), t)
-    t == 0 && return complex(1)
-    h = d.ν/2
-    q = d.ν/4
-    complex(2(q*t^2)^q * besselk(h, sqrt(d.ν) * abs(t)) / gamma(h))
+    h = d.ν / 2 
+    q = d.ν / 4
+    c = complex(2(q * t^2)^q * besselk(h, sqrt(d.ν) * abs(t)) / gamma(h))
+    iszero(t) ? complex(one(c)) : c
 end
 
 gradlogpdf(d::TDist{T}, x::Real) where {T<:Real} = isinf(d.ν) ? gradlogpdf(Normal(zero(T), one(T)), x) : -((d.ν + 1) * x) / (x^2 + d.ν)

src/univariate/continuous/triangular.jl

@@ -118,27 +118,18 @@ function quantile(d::TriangularDist, p::Real)
               b - sqrt(b_m_a * (b - c) * (1 - p))
 end
 
-function mgf(d::TriangularDist{T}, t::Real) where T<:Real
-    if t == zero(t)
-        return one(T)
-    else
-        (a, b, c) = params(d)
-        u = (b - c) * exp(a * t) - (b - a) * exp(c * t) + (c - a) * exp(b * t)
-        v = (b - a) * (c - a) * (b - c) * t^2
-        return 2u / v
-    end
+function mgf(d::TriangularDist, t::Number)
+    (a, b, c) = params(d)
+    u = (b - c) * exp(a * t) - (b - a) * exp(c * t) + (c - a) * exp(b * t)
+    v = (b - a) * (c - a) * (b - c) * t^2
+    return iszero(t) ? one(u) : 2u / v
 end
 
-function cf(d::TriangularDist{T}, t::Real) where T<:Real
-    # Is this correct?
-    if t == zero(t)
-        return one(Complex{T})
-    else
-        (a, b, c) = params(d)
-        u = (b - c) * cis(a * t) - (b - a) * cis(c * t) + (c - a) * cis(b * t)
-        v = (b - a) * (c - a) * (b - c) * t^2
-        return -2u / v
-    end
+function cf(d::TriangularDist, t::Number)
+    (a, b, c) = params(d)
+    u = (b - c) * cis(a * t) - (b - a) * cis(c * t) + (c - a) * cis(b * t)
+    v = (b - a) * (c - a) * (b - c) * t^2
+    return iszero(t) ? complex(one(u)) : -2u / v
 end
 
 

src/univariate/continuous/triweight.jl

@@ -59,14 +59,17 @@ end
 
 @quantile_newton Triweight
 
-function mgf(d::Triweight{T}, t::Float64) where T<:Real
-    a = d.σ*t
-    a2 = a*a
-    a == 0 ? one(T) : 105*exp(d.μ*t)*((15/a2+1)*cosh(a)-(15/a2-6)/a*sinh(a))/(a2*a2)
+function mgf(d::Triweight, t::Number)
+    a = d.σ * t
+    a2 = a * a
+    abs(real(a)) < 1 || throw(DomainError(t, "|σ ⋅ real(t)| should be smaller than 1"))
+    result = 105 * exp(d.μ * t) * ((15/a2 + 1) * cosh(a) - (15/a2 - 6) / a * sinh(a)) / (a2 * a2)
+    return iszero(t) ? one(result) : result
 end
 
-function cf(d::Triweight{T}, t::Float64) where T<:Real
-    a = d.σ*t
-    a2 = a*a
-    a == 0 ? complex(one(T)) : 105*cis(d.μ*t)*((1-15/a2)*cos(a)+(15/a2-6)/a*sin(a))/(a2*a2)
+function cf(d::Triweight, t::Number)
+    a = d.σ * t
+    a2 = a * a
+    result = 105 * cis(d.μ * t) * ((1 - 15/a2) * cos(a) + (15/a2 - 6) / a * sin(a)) / (a2 * a2)
+    return iszero(t) ? complex(one(result)) : result
 end

src/univariate/continuous/uniform.jl

@@ -93,20 +93,21 @@ quantile(d::Uniform, p::Real) = d.a + p * (d.b - d.a)
 cquantile(d::Uniform, p::Real) = d.b + p * (d.a - d.b)
 
 
-function mgf(d::Uniform, t::Real)
+function mgf(d::Uniform, t::Number)
     (a, b) = params(d)
     u = (b - a) * t / 2
-    u == zero(u) && return one(u)
+    iszero(u) && return one(float(typeof(u)))
+    abs(real(u)) < 1 || throw(DomainError(t, "|(b - a) ⋅ real(t) / 2| should be smaller than 1"))
     v = (a + b) * t / 2
-    exp(v) * (sinh(u) / u)
+    return exp(v) * (sinh(u) / u)
 end
 
-function cf(d::Uniform, t::Real)
+function cf(d::Uniform, t::Number)
     (a, b) = params(d)
     u = (b - a) * t / 2
-    u == zero(u) && return complex(one(u))
+    iszero(u) && return complex(one(float(typeof(u))))
     v = (a + b) * t / 2
-    cis(v) * (sin(u) / u)
+    return cis(v) * (sin(u) / u)
 end
 
 #### Affine transformations

src/univariate/discrete/bernoulli.jl

@@ -102,8 +102,8 @@ function cquantile(d::Bernoulli{T}, p::Real) where T<:Real
     0 <= p <= 1 ? (p >= succprob(d) ? zero(T) : one(T)) : T(NaN)
 end
 
-mgf(d::Bernoulli, t::Real) = failprob(d) + succprob(d) * exp(t)
-cf(d::Bernoulli, t::Real) = failprob(d) + succprob(d) * cis(t)
+mgf(d::Bernoulli, t::Number) = failprob(d) + succprob(d) * exp(t)
+cf(d::Bernoulli, t::Number) = failprob(d) + succprob(d) * cis(t)
 
 
 #### Sampling

src/univariate/discrete/binomial.jl

@@ -125,12 +125,12 @@ function rand(rng::AbstractRNG, d::Binomial)
     p <= 0.5 ? y : n-y
 end
 
-function mgf(d::Binomial, t::Real)
+function mgf(d::Binomial, t::Number)
     n, p = params(d)
     (one(p) - p + p * exp(t)) ^ n
 end
 
-function cf(d::Binomial, t::Real)
+function cf(d::Binomial, t::Number)
     n, p = params(d)
     (one(p) - p + p * cis(t)) ^ n
 end

src/univariate/discrete/dirac.jl

@@ -49,8 +49,8 @@ logccdf(d::Dirac, x::Real) = x < d.value ? 0.0 : isnan(x) ? NaN : -Inf
 
 quantile(d::Dirac{T}, p::Real) where {T} = 0 <= p <= 1 ? d.value : T(NaN)
 
-mgf(d::Dirac, t) = exp(t * d.value)
-cf(d::Dirac, t) = cis(t * d.value)
+mgf(d::Dirac, t::Number) = exp(t * d.value)
+cf(d::Dirac, t::Number) = cis(t * d.value)
 
 #### Sampling