Add mean et al. for truncated log normal

Fixes 709

Author: ararslan

src/truncate.jl

@@ -261,6 +261,7 @@ include(joinpath("truncated", "exponential.jl"))
 include(joinpath("truncated", "uniform.jl"))
 include(joinpath("truncated", "loguniform.jl"))
 include(joinpath("truncated", "discrete_uniform.jl"))
+include(joinpath("truncated", "lognormal.jl"))
 
 #### Utilities
 

src/truncated/lognormal.jl

@@ -0,0 +1,52 @@
+# Moments of the truncated log-normal can be computed directly from the moment generating
+# function of the truncated normal:
+#   Let Y ~ LogNormal(μ, σ) truncated to (a, b). Then log(Y) ~ Normal(μ, σ) truncated
+#   to (log(a), log(b)), and E[Y^n] = E[(e^log(Y))^n] = E[e^(nlog(Y))] = mgf(log(Y), n).
+
+# Given `truncate(LogNormal(μ, σ), a, b)`, return `truncate(Normal(μ, σ), log(a), log(b))`
+function _truncnorm(d::Truncated{<:LogNormal})
+    μ, σ = params(d.untruncated)
+    a = d.lower === nothing ? nothing : log(minimum(d))
+    b = d.upper === nothing ? nothing : log(maximum(d))
+    return truncated(Normal(μ, σ), a, b)
+end
+
+mean(d::Truncated{<:LogNormal}) = mgf(_truncnorm(d), 1)
+
+function var(d::Truncated{<:LogNormal})
+    tn = _truncnorm(d)
+    m1 = mgf(tn, 1)
+    m2 = sqrt(mgf(tn, 2))
+    return (m2 - m1) * (m2 + m1)
+end
+
+function skewness(d::Truncated{<:LogNormal})
+    tn = _truncnorm(d)
+    m1 = mgf(tn, 1)
+    m2 = sqrt(mgf(tn, 2))
+    m3 = mgf(tn, 3)
+    v = (m2 - m1) * (m2 + m1)
+    return (m3 - 3 * m1 * v - m1^3) / (v * sqrt(v))
+end
+
+function kurtosis(d::Truncated{<:LogNormal})
+    tn = _truncnorm(d)
+    m1 = mgf(tn, 1)
+    m2 = mgf(tn, 2)
+    m3 = mgf(tn, 3)
+    m4 = mgf(tn, 4)
+    sm2 = sqrt(m2)
+    v = (sm2 - m1) * (sm2 + m1)
+    return evalpoly(m1, (m4, -4m3, 6m2, 0, -3)) / v^2 - 3
+end
+
+median(d::Truncated{<:LogNormal}) = exp(median(_truncnorm(d)))
+
+# TODO: The entropy can be written "directly" as well, according to Mathematica, but
+# the expression for it fills me with regret. There are some recognizable components,
+# so a sufficiently motivated person could try to manually simplify it into something
+# comprehensible. For reference, you can obtain the entropy with Mathematica like so:
+#
+#   d = TruncatedDistribution[{a, b}, LogNormalDistribution[m, s]];
+#   Expectation[-LogLikelihood[d, {x}], Distributed[x, d],
+#               Assumptions -> Element[x | m | s | a | b, Reals] && s > 0 && 0 < a < x < b]

src/truncated/normal.jl

@@ -118,6 +118,16 @@ function entropy(d::Truncated{<:Normal{<:Real},Continuous})
     0.5 * (log2π + 1.) + log(σ * z) + (aφa - bφb) / (2.0 * z)
 end
 
+function mgf(d::Truncated{Normal{T},Continuous}, t::Real) where {T}
+    d0 = d.untruncated
+    μ = mean(d0)
+    σ = std(d0)
+    σt = σ * t
+    a = (minimum(d) - μ) / σ - σt
+    b = (maximum(d) - μ) / σ - σt
+    stdnorm = Normal{T}(zero(T), one(T))
+    return exp(μ * t + σt^2 / 2 + logdiffcdf(d0, b, a) - d.logtp)
+end
 
 ### sampling