Use Number instead of Real for mgf and cf

src/univariate/continuous/biweight.jl

@@ -54,16 +54,16 @@ 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
     a2 = a^2
-    a == 0 ? one(T) :
+    iszero(a) ? one(Base.promote_typeof(a, d.μ)) :
     15exp(d.μ * t) * (-3cosh(a) + (a + 3/a) * sinh(a)) / (a2^2)
 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 :
+    iszero(a) ? complex(one(Base.promote_typeof(a, d.μ))) :
     -15cis(d.μ * t) * (3cos(a) + (a - 3/a) * sin(a)) / (a2^2)
 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,9 @@ end
 
 @_delegate_statsfuns Chisq chisq ν
 
-mgf(d::Chisq, t::Real) = (1 - 2 * t)^(-d.ν/2)
+mgf(d::Chisq, t::Number) = (1 - 2 * t)^(-d.ν/2)
 
-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,14 @@ 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) :
+    iszero(a) ? one(Base.promote_typeof(a, d.μ)) :
     3exp(d.μ * t) * (cosh(a) - sinh(a) / a) / a^2
 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 :
+    iszero(a) ? complex(one(Base.promote_typeof(a, d.μ))) :
     -3exp(im * d.μ * t) * (cos(a) - sin(a) / a) / a^2
 end

src/univariate/continuous/erlang.jl

@@ -75,8 +75,8 @@ 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.α)
+mgf(d::Erlang, t::Number) = (1 - t * d.θ)^(-d.α)
+cf(d::Erlang, t::Number)  = (1 - im * t * d.θ)^(-d.α)
 
 
 #### Evaluation & Sampling

src/univariate/continuous/exponential.jl

@@ -93,8 +93,8 @@ 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))
+mgf(d::Exponential, t::Number) = 1/(1 - t * scale(d))
+cf(d::Exponential, t::Number) = 1/(1 - t * im * scale(d))
 
 
 #### Sampling

src/univariate/continuous/gamma.jl

@@ -75,9 +75,9 @@ function entropy(d::Gamma)
     α + loggamma(α) + (1 - α) * digamma(α) + log(θ)
 end
 
-mgf(d::Gamma, t::Real) = (1 - t * d.θ)^(-d.α)
+mgf(d::Gamma, t::Number) = (1 - t * d.θ)^(-d.α)
 
-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,15 @@ 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))
+    iszero(t) ? one(Base.promote_typeof(a, b, t)) : 2(-b*t)^(0.5a) / gamma(a) * besselk(a, sqrt(-4*b*t))
 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))
+    T = promote_rule(typeof(a), typeof(b))
+    iszero(t) ? complex(one(Base.promote_typeof(a, b, t))) : 2(-im*b*t)^(0.5a) / gamma(a) * besselk(a, sqrt(-4*im*b*t))
 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)
+    abs(t) < 1 / d.θ || return Base.promote_typeof(t, d.θ, d.μ)(NaN)
     st = d.θ * t
     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,9 @@ function gradlogpdf(d::Logistic, x::Real)
     ((2e) / (1 + e) - 1) / d.θ
 end
 
-mgf(d::Logistic, t::Real) = exp(t * d.μ) / sinc(d.θ * t)
+mgf(d::Logistic, t::Number) = exp(t * d.μ) / sinc(d.θ * t)
 
-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))
+    iszero(a) ? complex(one(a)) : cis(t * d.μ) * (a / sinh(a))
 end

src/univariate/continuous/noncentralchisq.jl

@@ -61,11 +61,11 @@ 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)
+function mgf(d::NoncentralChisq, t::Number)
     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)/sqrttwo) )
 
 #### Sampling
 function rand(rng::AbstractRNG, d::SkewNormal)

src/univariate/continuous/triangular.jl

@@ -118,23 +118,22 @@ 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)
+function mgf(d::TriangularDist{T}, t::Number) where T<:Real
+    (a, b, c) = params(d)
+    if iszero(t)
+        return one(Base.promote_typeof(a, b, c, 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
 end
 
-function cf(d::TriangularDist{T}, t::Real) where T<:Real
-    # Is this correct?
-    if t == zero(t)
-        return one(Complex{T})
+function cf(d::TriangularDist, t::Number)
+    (a, b, c) = params(d)
+    if iszero(t)
+        return complex(one(Base.promote_typeof(a, b, c, 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

src/univariate/continuous/triweight.jl

@@ -59,14 +59,22 @@ 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
+    if iszero(t)
+        return one(Base.promote_typeof(a, d.μ))
+    else
+        return 105 * exp(d.μ * t) * ((15/a2 + 1) * cosh(a) - (15/a2 - 6) / a * sinh(a)) / (a2 * a2)
+    end
 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
+    if iszero(t)
+        return complex(one(Base.promote_typeof(a, d.μ)))
+    else
+        return 105 * cis(d.μ * t) * ((1 - 15/a2) * cos(a) + (15/a2 - 6) / a * sin(a)) / (a2 * a2)
+    end
 end

src/univariate/continuous/uniform.jl

@@ -93,20 +93,20 @@ 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)))
     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
 

src/univariate/discrete/discretenonparametric.jl

@@ -213,22 +213,20 @@ function modes(d::DiscreteNonParametric)
     return modes(x, Weights(p, one(eltype(p))))
 end
 
-function mgf(d::DiscreteNonParametric, t::Real)
-    x, p = params(d)
-    s = zero(Float64)
-    for i in 1:length(x)
-        s += p[i] * exp(t*x[i])
-    end
-    s
+function mgf(d::DiscreteNonParametric, t::Number)
+    xs, ps = params(d)
+    s = sum(Broadcast.instantiate(Broadcast.broadcasted(ps, xs) do p, x
+        return p * exp(t * x)
+    end))
+    return s
 end
 
-function cf(d::DiscreteNonParametric, t::Real)
-    x, p = params(d)
-    s = zero(Complex{Float64})
-    for i in 1:length(x)
-       s += p[i] * cis(t*x[i])
-    end
-    s
+function cf(d::DiscreteNonParametric, t::Number)
+    xs, ps = params(d)
+    s = sum(Broadcast.instantiate(Broadcast.broadcasted(ps, xs) do p, x
+        return p * cis(t * x)
+    end))
+    return s
 end
 
 # Sufficient statistics

src/univariate/discrete/discreteuniform.jl

@@ -87,14 +87,14 @@ end
 
 quantile(d::DiscreteUniform, p::Real) = iszero(p) ? d.a : d.a - 1 + ceil(Int, p * span(d))
 
-function mgf(d::DiscreteUniform, t::Real)
+function mgf(d::DiscreteUniform, t::Number)
     a, b = d.a, d.b
     u = b - a + 1
     result = (exp(t*a) * expm1(t*u)) / (u*expm1(t))
     return iszero(t) ? one(result) : result
 end
 
-function cf(d::DiscreteUniform, t::Real)
+function cf(d::DiscreteUniform, t::Number)
     a, b = d.a, d.b
     u = b - a + 1
     result = (im*cos(t*(a+b)/2) + sin(t*(a-b-1)/2)) / (u*sin(t/2))

src/univariate/discrete/geometric.jl

@@ -109,12 +109,12 @@ function invlogccdf(d::Geometric{T}, lp::Real) where T<:Real
     max(ceil(lp/log1p(-d.p)) - 1, zero(T))
 end
 
-function mgf(d::Geometric, t::Real)
+function mgf(d::Geometric, t::Number)
     p = succprob(d)
     p / (expm1(-t) + p)
 end
 
-function cf(d::Geometric, t::Real)
+function cf(d::Geometric, t::Number)
     p = succprob(d)
     # replace with expm1 when complex version available
     p / (exp(-t*im) - 1 + p)