Addition of heslogpdf of cont. univariate distributions

Addition of the heslogpdf, i.e., second derivative of the logpdf, and associated tests. Included distributions: arcsine, beta, chi, chisq, exponential, frechet, gamma, gumbel, laplace, locationscale, logistic, logitnormal, lognormal, normal, tdist and weibull.

docs/src/univariate.md

@@ -71,6 +71,7 @@ cf(::UnivariateDistribution, ::Any)
 insupport(::UnivariateDistribution, x::Any)
 pdf(::UnivariateDistribution, ::Real)
 logpdf(::UnivariateDistribution, ::Real)
+heslogpdf(::ContinuousUnivariateDistribution, ::Real)
 loglikelihood(::UnivariateDistribution, ::Union{Real,AbstractArray})
 cdf(::UnivariateDistribution, ::Real)
 logcdf(::UnivariateDistribution, ::Real)

src/Distributions.jl

@@ -250,6 +250,7 @@ export
     varlogx,            # variance of log(x)
     expected_logdet,    # expected logarithm of random matrix determinant
     gradlogpdf,         # gradient (or derivative) of logpdf(d,x) wrt x
+    heslogpdf,          # hessian (or second derivative) of logpdf(d,x) wrt x
 
     # reexport from StatsBase
     sample, sample!,        # sample from a source array

src/univariate/continuous/arcsine.jl

@@ -99,3 +99,13 @@ function gradlogpdf(d::Arcsine{T}, x::R) where {T, R <: Real}
     x == b && return TP(Inf)
     return TP(0.5 * (inv(b - x) - inv(x - a)))
 end
+function heslogpdf(d::Arcsine{T}, x::R) where {T, R <: Real}
+    TP = promote_type(T, R)
+    (a, b) = extrema(d)
+    # on the bounds, we consider the gradient limit inside the domain
+    # right side for the left bound,
+    # left side for the right bound
+    a < x <= b || return TP(-Inf)
+    x == b && return TP(Inf)
+    return TP(0.5 * (inv((x - a)^2) + inv((x - b)^2)))
+end

src/univariate/continuous/beta.jl

@@ -114,7 +114,8 @@ end
 
 gradlogpdf(d::Beta{T}, x::Real) where {T<:Real} =
     ((α, β) = params(d); 0 <= x <= 1 ? (α - 1) / x - (β - 1) / (1 - x) : zero(T))
-
+heslogpdf(d::Beta{T}, x::Real) where {T<:Real} =
+    ((α, β) = params(d); 0 <= x <= 1 ? - (α - 1) / x^2 - (β - 1) / (1 - x)^2 : zero(T))
 
 #### Sampling
 

src/univariate/continuous/chi.jl

@@ -82,6 +82,7 @@ logpdf(d::Chi, x::Real) = (ν = d.ν;
 )
 
 gradlogpdf(d::Chi{T}, x::Real) where {T<:Real} = x >= 0 ? (d.ν - 1) / x - x : zero(T)
+heslogpdf(d::Chi{T}, x::Real) where {T<:Real} = x >= 0 ? - (d.ν - 1) / x^2 - 1 : zero(T)
 
 cdf(d::Chi, x::Real) = chisqcdf(d.ν, x^2)
 ccdf(d::Chi, x::Real) = chisqccdf(d.ν, x^2)

src/univariate/continuous/chisq.jl

@@ -79,8 +79,8 @@ mgf(d::Chisq, t::Real) = (1 - 2 * t)^(-d.ν/2)
 
 cf(d::Chisq, t::Real) = (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)
-
+gradlogpdf(d::Chisq{T}, x::Real) where {T<:Real} = x > 0 ? (d.ν/2 - 1) / x - 1//2 : zero(T)
+heslogpdf(d::Chisq{T}, x::Real) where {T<:Real} = x > 0 ? - (d.ν/2 - 1) / x^2 : zero(T)
 
 #### Sampling
 

src/univariate/continuous/exponential.jl

@@ -82,6 +82,7 @@ invlogcdf(d::Exponential, lp::Real) = -xval(d, log1mexp(lp))
 invlogccdf(d::Exponential, lp::Real) = -xval(d, lp)
 
 gradlogpdf(d::Exponential{T}, x::Real) where {T<:Real} = x > 0 ? -rate(d) : zero(T)
+heslogpdf(d::Exponential{T}, x::Real) where {T<:Real} = zero(T)
 
 mgf(d::Exponential, t::Real) = 1/(1 - t * scale(d))
 cf(d::Exponential, t::Real) = 1/(1 - t * im * scale(d))

src/univariate/continuous/frechet.jl

@@ -131,7 +131,11 @@ invlogccdf(d::Frechet, lp::Real) = d.θ * (-log1mexp(lp))^(-1 / d.α)
 
 function gradlogpdf(d::Frechet{T}, x::Real) where T<:Real
     (α, θ) = params(d)
-    insupport(Frechet, x) ? -(α + 1) / x + α * (θ^α) * x^(-α-1)  : zero(T)
+    insupport(Frechet, x) ? -(α + 1) / x + α * (θ^α) * x^(-α-1) : zero(T)
+end
+function heslogpdf(d::Frechet{T}, x::Real) where T<:Real
+    (α, θ) = params(d)
+    insupport(Frechet, x) ? -(α + 1) / x^2 + α * (θ^α) * (-α-1) * x^(-α-2) : zero(T)
 end
 
 ## Sampling

src/univariate/continuous/gamma.jl

@@ -86,6 +86,8 @@ cf(d::Gamma, t::Real) = (1 - im * t * d.θ)^(-d.α)
 
 gradlogpdf(d::Gamma{T}, x::Real) where {T<:Real} =
     insupport(Gamma, x) ? (d.α - 1) / x - 1 / d.θ : zero(T)
+heslogpdf(d::Gamma{T}, x::Real) where {T<:Real} =
+    insupport(Gamma, x) ? - (d.α - 1) / x^2 : zero(T)
 
 function rand(rng::AbstractRNG, d::Gamma)
     if shape(d) < 1.0

src/univariate/continuous/gumbel.jl

@@ -93,3 +93,4 @@ logcdf(d::Gumbel, x::Real) = -exp(-zval(d, x))
 quantile(d::Gumbel, p::Real) = d.μ - d.θ * log(-log(p))
 
 gradlogpdf(d::Gumbel, x::Real) = - (1 + exp((d.μ - x) / d.θ)) / d.θ
+heslogpdf(d::Gumbel, x::Real) = exp((d.μ - x) / d.θ) / d.θ^2

src/univariate/continuous/laplace.jl

@@ -95,6 +95,7 @@ function gradlogpdf(d::Laplace, x::Real)
     g = 1 / θ
     x > μ ? -g : g
 end
+heslogpdf(d::Laplace{T}, x::Real) where {T<:Real} = zero(T)
 
 function mgf(d::Laplace, t::Real)
     st = d.θ * t

src/univariate/continuous/locationscale.jl

@@ -81,3 +81,4 @@ quantile(d::LocationScale,q::Real) = d.μ + d.σ * quantile(d.ρ,q)
 rand(rng::AbstractRNG, d::LocationScale) = d.μ + d.σ * rand(rng, d.ρ)
 cf(d::LocationScale, t::Real) = cf(d.ρ,t*d.σ) * exp(1im*t*d.μ)
 gradlogpdf(d::LocationScale, x::Real) = gradlogpdf(d.ρ,(x-d.μ)/d.σ) / d.σ
+heslogpdf(d::LocationScale, x::Real) = heslogpdf(d.ρ,(x-d.μ)/d.σ) / d.σ^2

src/univariate/continuous/logistic.jl

@@ -95,6 +95,10 @@ function gradlogpdf(d::Logistic, x::Real)
     e = exp(-zval(d, x))
     ((2e) / (1 + e) - 1) / d.θ
 end
+function heslogpdf(d::Logistic, x::Real)
+    e = exp(-zval(d,x))
+    - ((2e) / (1 + e)^2) / d.θ^2
+end
 
 mgf(d::Logistic, t::Real) = exp(t * d.μ) / sinc(d.θ * t)
 

src/univariate/continuous/logitnormal.jl

@@ -141,6 +141,11 @@ function gradlogpdf(d::LogitNormal{T}, x::Real) where T<:Real
     (μ, σ) = params(d)
     0 < x < 1 ? - ((log(x) - μ) / ((σ^2) + 1) * x * (1-x)) : zero(T)
 end
+function heslogpdf(d::LogitNormal{T}, x::Real) where T<:Real
+    #TODO check
+    (μ, σ) = params(d)
+    0 < x < 1 ? ((2 * x - 1) * log(x) + (1 - 2 * μ) * x + μ - 1) / (σ^2 + 1) : zero(T)
+end
 
 # mgf(d::LogitNormal)
 # cf(d::LogitNormal)

src/univariate/continuous/lognormal.jl

@@ -144,6 +144,8 @@ function gradlogpdf(d::LogNormal, x::Real)
     z = (gradlogpdf(Normal(d.μ, d.σ), log(y)) - 1) / y
     return outofsupport ? zero(z) : z
 end
+heslogpdf(d::LogNormal, x::Real) =
+    x > zero(x) ? (log(x) + d.σ^2 - d.μ - 1) / (d.σ^2 * x^2) : zero(x)
 
 # mgf(d::LogNormal)
 # cf(d::LogNormal)

src/univariate/continuous/normal.jl

@@ -98,6 +98,7 @@ Computes the z-value based on a Normal distribution and a x-value.
 zval(d::Normal, x::Real) = (x - d.μ) / d.σ
 
 gradlogpdf(d::Normal, x::Real) = -zval(d, x) / d.σ
+heslogpdf(d::Normal, x::Real) = - 1 / d.σ^2
 
 # logpdf
 _normlogpdf(z::Real) = -(abs2(z) + log2π)/2

src/univariate/continuous/tdist.jl

@@ -88,4 +88,7 @@ function cf(d::TDist{T}, t::Real) where T <: Real
     complex(2(q*t^2)^q * besselk(h, sqrt(d.ν) * abs(t)) / gamma(h))
 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.ν)
+gradlogpdf(d::TDist{T}, x::Real) where {T<:Real} = 
+    isinf(d.ν) ? gradlogpdf(Normal(zero(T), one(T)), x) : -((d.ν + 1) * x) / (x^2 + d.ν)
+heslogpdf(d::TDist{T}, x::Real) where {T<:Real} = 
+    isinf(d.ν) ? heslogpdf(Normal(zero(T), one(T)), x) : (d.ν + 1) * (x^2 - d.ν) / (x^2 + d.ν)^2

src/univariate/continuous/weibull.jl

@@ -134,7 +134,14 @@ function gradlogpdf(d::Weibull{T}, x::Real) where T<:Real
         zero(T)
     end
 end
-
+function heslogpdf(d::Weibull{T}, x::Real) where T<:Real
+    if insupport(Weibull, x)
+        α, θ = params(d)
+        - (α - 1) / x^2 - α * (α - 1) * x^(α - 2) / (θ^α)
+    else
+        zero(T)
+    end
+end
 
 #### Sampling
 

src/univariates.jl

@@ -440,6 +440,13 @@ invlogccdf(d::UnivariateDistribution, lp::Real) = quantile(d, -expm1(lp))
 
 gradlogpdf(d::ContinuousUnivariateDistribution, x::Real) = throw(MethodError(gradlogpdf, (d, x)))
 
+"""
+    heslogpdf(d::ContinuousUnivariateDistribution, x::Real)
+
+The second derivative of the logarithm of probability density (mass) w.r.t `x` evaluated at `x`.
+"""
+heslogpdf(d::ContinuousUnivariateDistribution, x::Real) = throw(MethodError(heslogpdf, (d, x)))
+
 
 function _pdf_fill_outside!(r::AbstractArray, d::DiscreteUnivariateDistribution, X::UnitRange)
     vl = vfirst = first(X)

test/arcsine.jl

@@ -24,4 +24,11 @@ using ForwardDiff
         glog = gradlogpdf(d, v)
         @test fgrad ≈ glog
     end
+
+    # second derivative
+    for v in 3.01:0.1:4.99
+        fhes = ForwardDiff.derivative(x -> gradlogpdf(d, x), v)
+        glog = heslogpdf(d, v)
+        @test fhes ≈ glog
+    end
 end

test/edgecases.jl

@@ -17,6 +17,7 @@ using Test, Distributions, StatsBase
     @test pdf(T, x) ≈ pdf(N, x)
     @test logpdf(T, x) ≈ logpdf(N, x)
     @test gradlogpdf(T, x) ≈ gradlogpdf(N, x)
+    @test heslogpdf(T, x) ≈ heslogpdf(N, x)
     @test cf(T, x) ≈ cf(N, x)
 
     fnecdf = ecdf(z)

test/heslogpdf.jl

@@ -0,0 +1,17 @@
+using Distributions
+using Test
+
+# Test for heslogpdf on univariate distributions
+
+@test isapprox(heslogpdf(Beta(1.5, 3.0), 0.7)   , -23.242630385487523 , atol=1.0e-8)
+@test isapprox(heslogpdf(Chi(5.0), 5.5)         , -1.1322314049586777 , atol=1.0e-8)
+@test isapprox(heslogpdf(Chisq(7.0), 12.0)      , -0.01736111111111111, atol=1.0e-8)
+@test isapprox(heslogpdf(Exponential(2.0), 7.0) ,  0.0                , atol=1.0e-8)
+@test isapprox(heslogpdf(Gamma(9.0, 0.5), 11.0) , -0.06611570247933884, atol=1.0e-8)
+@test isapprox(heslogpdf(Gumbel(3.5, 1.0), 4.0) ,  0.6065306597126334 , atol=1.0e-8)
+@test isapprox(heslogpdf(Laplace(7.0), 34.0)    ,  0.0                , atol=1.0e-8)
+@test isapprox(heslogpdf(Logistic(-6.0), 1.0)   , -0.00182044236024365, atol=1.0e-8)
+@test isapprox(heslogpdf(LogNormal(5.5), 2.0)   , -1.2017132048600137 , atol=1.0e-8)
+@test isapprox(heslogpdf(Normal(-4.5, 2.0), 1.6), -0.25               , atol=1.0e-8)
+@test isapprox(heslogpdf(TDist(8.0), 9.1)       ,  0.0816459812355031 , atol=1.0e-8)
+@test isapprox(heslogpdf(Weibull(2.0), 3.5)     , -2.0816326530612246 , atol=1.0e-8)

test/locationscale.jl

@@ -72,6 +72,7 @@ function test_location_scale_normal(μ::Real, σ::Real, μD::Real, σD::Real,
     @test std(r) ≈ std(dref) atol=0.01
     @test cf(d, -0.1) ≈ cf(dref,-0.1)
     @test gradlogpdf(d, 0.1) ≈ gradlogpdf(dref, 0.1)
+    @test heslogpdf(d, 0.1) ≈ heslogpdf(dref, 0.1)
 
 end
 

test/lognormal.jl

@@ -285,6 +285,12 @@ end
     @test @inferred(gradlogpdf(LogNormal(0.0f0, 1.0f0), 1.0)) === -1.0
     @test @inferred(gradlogpdf(LogNormal(0.0f0, 1.0f0), 1.0f0)) === -1.0f0
 
+    # heslogpdf
+    @test @inferred(heslogpdf(LogNormal(0.0, 1.0), 1.0)) === 0.0
+    @test @inferred(heslogpdf(LogNormal(0.0, 1.0), exp(-1))) ≈ -7.3890560989306495
+    @test @inferred(heslogpdf(LogNormal(0.0, 1.0), 0.0)) === 0.0
+    @test @inferred(heslogpdf(LogNormal(0.0, 1.0), -0.5)) === 0.0
+
     @test isnan_type(Float64, @inferred(logccdf(LogNormal(0.0, 1.0), NaN)))
     @test isnan_type(Float64, @inferred(logccdf(LogNormal(NaN, 1.0), 1.0f0)))
     @test isnan_type(Float64, @inferred(logccdf(LogNormal(NaN, 1.0), 0.0f0)))

test/runtests.jl

@@ -43,6 +43,7 @@ const tests = [
     "convolution",
     "mixture",
     "gradlogpdf",
+    "heslogpdf",
     "noncentralt",
     "locationscale",
     "quantile_newton",