Add unnormalized log pdf functions

Author: devmotion

ext/StatsFunsChainRulesCoreExt.jl

@@ -13,10 +13,25 @@ ChainRulesCore.@scalar_rule(
         (α - 1) / x + (1 - β) / (1 - x),
     ),
 )
+ChainRulesCore.@scalar_rule(
+    betalogupdf(α::Real, β::Real, x::Number),
+    (
+        log(x),
+        log1p(-x),
+        (α - 1) / x + (1 - β) / (1 - x),
+    ),
+)
 
 ChainRulesCore.@scalar_rule(
     binomlogpdf(n::Real, p::Real, k::Real),
-    @setup(z = digamma(n - k + 1)),
+    (
+        ChainRulesCore.NoTangent(),
+        (k / p - n) / (1 - p),
+        ChainRulesCore.NoTangent(),
+    ),
+)
+ChainRulesCore.@scalar_rule(
+    binomlogupdf(n::Real, p::Real, k::Real),
     (
         ChainRulesCore.NoTangent(),
         (k / p - n) / (1 - p),
@@ -28,7 +43,15 @@ ChainRulesCore.@scalar_rule(
     chisqlogpdf(k::Real, x::Number),
     @setup(hk = k / 2),
     (
-        (log(x) - logtwo - digamma(hk)) / 2,
+        (log(x / 2) - digamma(hk)) / 2,
+        (hk - 1) / x - one(hk) / 2,
+    ),
+)
+ChainRulesCore.@scalar_rule(
+    chisqlogupdf(k::Real, x::Number),
+    @setup(hk = k / 2),
+    (
+        log(x / 2) / 2,
         (hk - 1) / x - one(hk) / 2,
     ),
 )
@@ -48,6 +71,19 @@ ChainRulesCore.@scalar_rule(
         ((ν1 - 2) / x - ν1 * νsum / temp1) / 2,
     ),
 )
+ChainRulesCore.@scalar_rule(
+    fdistlogupdf(ν1::Real, ν2::Real, x::Number),
+    @setup(
+        tmp = x * ν1 + ν2,
+        a = x * (ν1 + ν2) / tmp,
+        b = ν2 / tmp,
+    ),
+    (
+        (log(x * b) - a) / 2,
+        (log(b) + (ν1 / ν2) * a) / 2,
+        (ν1 * (1 - a) - 2) / (2 * x),
+    ),
+)
 
 ChainRulesCore.@scalar_rule(
     gammalogpdf(k::Real, θ::Real, x::Number),
@@ -62,11 +98,32 @@ ChainRulesCore.@scalar_rule(
         - (1 + z) / x,
     ),
 )
+ChainRulesCore.@scalar_rule(
+    gammalogupdf(k::Real, θ::Real, x::Number),
+    @setup(
+        invθ = inv(θ),
+        xoθ = invθ * x,
+        z = xoθ - (k - 1),
+    ),
+    (
+        log(xoθ),
+        invθ * z,
+        - z / x,
+    ),
+)
 
 ChainRulesCore.@scalar_rule(
     poislogpdf(λ::Number, x::Number),
     ((iszero(x) && iszero(λ) ? zero(x / λ) : x / λ) - 1, ChainRulesCore.NoTangent()),
 )
+ChainRulesCore.@scalar_rule(
+    poislogupdf(λ::Number, x::Number),
+    ((iszero(x) && iszero(λ) ? zero(x / λ) : x / λ), ChainRulesCore.NoTangent()),
+)
+ChainRulesCore.@scalar_rule(
+    poislogulikelihood(λ::Number, x::Number),
+    ((iszero(x) && iszero(λ) ? zero(x / λ) : x / λ) - 1, ChainRulesCore.NoTangent()),
+)
 
 ChainRulesCore.@scalar_rule(
     tdistlogpdf(ν::Real, x::Number),
@@ -82,5 +139,17 @@ ChainRulesCore.@scalar_rule(
         - x * b,
     ),
 )
+ChainRulesCore.@scalar_rule(
+    tdistlogupdf(ν::Real, x::Number),
+    @setup(
+        xsq = x^2,
+        a = xsq / ν,
+        b = (ν + 1) / (ν + xsq),
+    ),
+    (
+        (a * b - log1p(a)) / 2,
+        - x * b,
+    ),
+)
 
 end # module

src/StatsFuns.jl

@@ -54,6 +54,8 @@ export
     # distrs/beta
     betapdf, # pdf of beta distribution
     betalogpdf, # logpdf of beta distribution
+    betalogupdf, # unnormalized logpdf of beta distribution (parameters constant)
+    betalogulikelihood, # unnormalized logpdf of beta distribution (data constant)
     betacdf, # cdf of beta distribution
     betaccdf, # ccdf of beta distribution
     betalogcdf, # logcdf of beta distribution
@@ -66,6 +68,8 @@ export
     # distrs/binom
     binompdf, # pdf of binomial distribution
     binomlogpdf, # logpdf of binomial distribution
+    binomlogupdf, # unnormalized logpdf of binomial distribution (parameters constant)
+    binomlogulikelihood, # unnormalized logpdf of binomial distribution (data constant)
     binomcdf, # cdf of binomial distribution
     binomccdf, # ccdf of binomial distribution
     binomlogcdf, # logcdf of binomial distribution
@@ -78,6 +82,8 @@ export
     # distrs/chisq
     chisqpdf, # pdf of chi-square distribution
     chisqlogpdf, # logpdf of chi-square distribution
+    chisqlogupdf, # unnormalized logpdf of chi-square distribution (parameters constant)
+    chisqlogulikelihood, # unnormalized logpdf of chi-square distribution (data constant)
     chisqcdf, # cdf of chi-square distribution
     chisqccdf, # ccdf of chi-square distribution
     chisqlogcdf, # logcdf of chi-square distribution
@@ -90,6 +96,8 @@ export
     # distrs/fdist
     fdistpdf, # pdf of F distribution
     fdistlogpdf, # logpdf of F distribution
+    fdistlogupdf, # unnormalized logpdf of F distribution (parameters constant)
+    fdistlogulikelihood, # unnormalized logpdf of F distribution (data constant)
     fdistcdf, # cdf of F distribution
     fdistccdf, # ccdf of F distribution
     fdistlogcdf, # logcdf of F distribution
@@ -102,6 +110,8 @@ export
     # distrs/gamma
     gammapdf, # pdf of gamma distribution
     gammalogpdf, # logpdf of gamma distribution
+    gammalogupdf, # unnormalized logpdf of gamma distribution (parameters constant)
+    gammalogulikelihood, # unnormalized logpdf of gamma distribution (data constant)
     gammacdf, # cdf of gamma distribution
     gammaccdf, # ccdf of gamma distribution
     gammalogcdf, # logcdf of gamma distribution
@@ -114,6 +124,8 @@ export
     # distrs/hyper
     hyperpdf, # pdf of hypergeometric distribution
     hyperlogpdf, # logpdf of hypergeometric distribution
+    hyperlogupdf, # unnormalized logpdf of hypergeometric distribution (parameters constant)
+    hyperlogulikelihood, # unnormalized logpdf of hypergeometric distribution (data constant)
     hypercdf, # cdf of hypergeometric distribution
     hyperccdf, # ccdf of hypergeometric distribution
     hyperlogcdf, # logcdf of hypergeometric distribution
@@ -126,6 +138,8 @@ export
     # distrs/nbeta
     nbetapdf, # pdf of noncentral beta distribution
     nbetalogpdf, # logpdf of noncentral beta distribution
+    nbetalogupdf, # unnormalized logpdf of noncentral beta distribution (parameters constant)
+    nbetalogulikelihood, # unnormalized logpdf of noncentral beta distribution (data constant)
     nbetacdf, # cdf of noncentral beta distribution
     nbetaccdf, # ccdf of noncentral beta distribution
     nbetalogcdf, # logcdf of noncentral beta distribution
@@ -138,6 +152,8 @@ export
     # distrs/nbinom
     nbinompdf, # pdf of negative nbinomial distribution
     nbinomlogpdf, # logpdf of negative nbinomial distribution
+    nbinomlogupdf, # unnormalized logpdf of negative nbinomial distribution (parameters constant)
+    nbinomlogulikelihood, # unnormalized logpdf of negative nbinomial distribution (data constant)
     nbinomcdf, # cdf of negative nbinomial distribution
     nbinomccdf, # ccdf of negative nbinomial distribution
     nbinomlogcdf, # logcdf of negative nbinomial distribution
@@ -150,6 +166,8 @@ export
     # distrs/nchisq
     nchisqpdf, # pdf of noncentral chi-square distribution
     nchisqlogpdf, # logpdf of noncentral chi-square distribution
+    nchisqlogupdf, # unnormalized logpdf of noncentral chi-square distribution (parameters constant)
+    nchisqlogulikelihood, # unnormalized logpdf of noncentral chi-square distribution (data constant)
     nchisqcdf, # cdf of noncentral chi-square distribution
     nchisqccdf, # ccdf of noncentral chi-square distribution
     nchisqlogcdf, # logcdf of noncentral chi-square distribution
@@ -162,6 +180,8 @@ export
     # distrs/nfdist
     nfdistpdf, # pdf of noncentral F distribution
     nfdistlogpdf, # logpdf of noncentral F distribution
+    nfdistlogupdf, # unnormalized logpdf of noncentral F distribution (parameters constant)
+    nfdistlogulikelihood, # unnormalized logpdf of noncentral F distribution (data constant)
     nfdistcdf, # cdf of noncentral F distribution
     nfdistccdf, # ccdf of noncentral F distribution
     nfdistlogcdf, # logcdf of noncentral F distribution
@@ -174,6 +194,8 @@ export
     # distrs/norm
     normpdf, # pdf of normal distribution
     normlogpdf, # logpdf of normal distribution
+    normlogupdf, # unnormalized logpdf of normal distribution (parameters constant)
+    normlogulikelihood, # unnormalized logpdf of normal distribution (data constant)
     normcdf, # cdf of normal distribution
     normccdf, # ccdf of normal distribution
     normlogcdf, # logcdf of normal distribution
@@ -186,6 +208,8 @@ export
     # distrs/ntdist
     ntdistpdf, # pdf of noncentral t distribution
     ntdistlogpdf, # logpdf of noncentral t distribution
+    ntdistlogupdf, # unnormalized logpdf of noncentral t distribution (parameters constant)
+    ntdistlogulikelihood, # unnormalized logpdf of noncentral t distribution (data constant)
     ntdistcdf, # cdf of noncentral t distribution
     ntdistccdf, # ccdf of noncentral t distribution
     ntdistlogcdf, # logcdf of noncentral t distribution
@@ -198,6 +222,8 @@ export
     # distrs/pois
     poispdf, # pdf of Poisson distribution
     poislogpdf, # logpdf of Poisson distribution
+    poislogupdf, # unnormalized logpdf of Poisson distribution (parameters constant)
+    poislogulikelihood, # unnormalized logpdf of Poisson distribution (data constant)
     poiscdf, # cdf of Poisson distribution
     poisccdf, # ccdf of Poisson distribution
     poislogcdf, # logcdf of Poisson distribution
@@ -210,6 +236,8 @@ export
     # distrs/tdist
     tdistpdf, # pdf of student's t distribution
     tdistlogpdf, # logpdf of student's t distribution
+    tdistlogupdf, # unnormalized logpdf of student's t distribution (parameters constant)
+    tdistlogulikelihood, # unnormalized logpdf of student's t distribution (data constant)
     tdistcdf, # cdf of student's t distribution
     tdistccdf, # ccdf of student's t distribution
     tdistlogcdf, # logcdf of student's t distribution
@@ -222,6 +250,8 @@ export
     # distrs/signrank
     signrankpdf,
     signranklogpdf,
+    signranklogupdf,
+    signranklogulikelihood,
     signrankcdf,
     signranklogcdf,
     signrankccdf,
@@ -244,6 +274,8 @@ export
     # distrs/wilcox
     wilcoxpdf,
     wilcoxlogpdf,
+    wilcoxlogupdf,
+    wilcoxlogulikelihood,
     wilcoxcdf,
     wilcoxlogcdf,
     wilcoxccdf,

src/distrs/beta.jl

@@ -20,12 +20,20 @@ betapdf(α::Real, β::Real, x::Real) = exp(betalogpdf(α, β, x))
 
 betalogpdf(α::Real, β::Real, x::Real) = betalogpdf(promote(α, β, x)...)
 function betalogpdf(α::T, β::T, x::T) where {T <: Real}
+    logupdf = betalogupdf(α, β, x)
+    return isfinite(logupdf) ? logupdf - logbeta(α, β) : logupdf
+end
+
+betalogupdf(α::Real, β::Real, x::Real) = betalogupdf(promote(α, β, x)...)
+function betalogupdf(α::T, β::T, x::T) where {T <: Real}
     # we ensure that `log(x)` and `log1p(-x)` do not error
     y = clamp(x, 0, 1)
-    val = xlogy(α - 1, y) + xlog1py(β - 1, -y) - logbeta(α, β)
+    val = xlogy(α - 1, y) + xlog1py(β - 1, -y)
     return x < 0 || x > 1 ? oftype(val, -Inf) : val
 end
 
+betalogulikelihood(α::Real, β::Real, x::Real) = betalogpdf(α, β, x)
+
 function betacdf(α::Real, β::Real, x::Real)
     # Handle degenerate cases
     if iszero(α) && β > 0

src/distrs/binom.jl

@@ -18,11 +18,23 @@ binompdf(n::Real, p::Real, k::Real) = exp(binomlogpdf(n, p, k))
 
 binomlogpdf(n::Real, p::Real, k::Real) = binomlogpdf(promote(n, p, k)...)
 function binomlogpdf(n::T, p::T, k::T) where {T <: Real}
+    logupdf = binomlogupdf(n, p, k)
+    if isfinite(logupdf)
+        return min(0, logupdf - log(n + 1))
+    else
+        return logupdf
+    end
+end
+
+binomlogupdf(n::Real, p::Real, k::Real) = binomlogupdf(promote(n, p, k)...)
+function binomlogupdf(n::T, p::T, k::T) where {T <: Real}
     m = clamp(k, 0, n)
-    val = min(0, betalogpdf(m + 1, n - m + 1, p) - log(n + 1))
+    val = betalogpdf(m + 1, n - m + 1, p)
     return 0 <= k <= n && isinteger(k) ? val : oftype(val, -Inf)
 end
 
+binomlogulikelihood(n::Real, p::Real, k::Real) = binomlogpdf(n, p, k)
+
 for l in ("", "log"), compl in (false, true)
     fbinom = Symbol(string("binom", l, ifelse(compl, "c", ""), "cdf"))
     fbeta = Symbol(string("beta", l, ifelse(compl, "", "c"), "cdf"))

src/distrs/chisq.jl

@@ -1,7 +1,7 @@
 # functions related to chi-square distribution
 
 # Just use the Gamma definitions
-for f in ("pdf", "logpdf", "cdf", "ccdf", "logcdf", "logccdf", "invcdf", "invccdf", "invlogcdf", "invlogccdf")
+for f in ("pdf", "logpdf", "logupdf", "logulikelihood", "cdf", "ccdf", "logcdf", "logccdf", "invcdf", "invccdf", "invlogcdf", "invlogccdf")
     _chisqf = Symbol("chisq" * f)
     _gammaf = Symbol("gamma" * f)
     @eval begin

src/distrs/fdist.jl

@@ -12,6 +12,17 @@ function fdistlogpdf(ν1::T, ν2::T, x::T) where {T <: Real}
     return x < 0 ? oftype(val, -Inf) : val
 end
 
+fdistlogupdf(ν1::Real, ν2::Real, x::Real) = fdistlogupdf(promote(ν1, ν2, x)...)
+function fdistlogupdf(ν1::T, ν2::T, x::T) where {T <: Real}
+    # we ensure that `log(x)` does not error if `x < 0`
+    ν1ν2 = ν1 / ν2
+    y = max(x, 0)
+    val = (xlogy(ν1 - 2, y) - xlogy(ν1 + ν2, 1 + ν1ν2 * y)) / 2
+    return x < 0 ? oftype(val, -Inf) : val
+end
+
+fdistlogulikelihood(ν1::Real, ν2::Real, x::Real) = fdistlogpdf(ν1, ν2, x)
+
 for f in ("cdf", "ccdf", "logcdf", "logccdf")
     ff = Symbol("fdist" * f)
     bf = Symbol("beta" * f)

src/distrs/gamma.jl

@@ -20,9 +20,15 @@ gammapdf(k::Real, θ::Real, x::Real) = exp(gammalogpdf(k, θ, x))
 
 gammalogpdf(k::Real, θ::Real, x::Real) = gammalogpdf(promote(k, θ, x)...)
 function gammalogpdf(k::T, θ::T, x::T) where {T <: Real}
+    logupdf = gammalogupdf(k, θ, x)
+    return isfinite(logupdf) ? logupdf - loggamma(k) - log(θ) : logupdf
+end
+
+gammalogupdf(k::Real, θ::Real, x::Real) = gammalogupdf(promote(k, θ, x)...)
+function gammalogupdf(k::T, θ::T, x::T) where {T <: Real}
     # we ensure that `log(x)` does not error if `x < 0`
     xθ = max(x, 0) / θ
-    val = -loggamma(k) - log(θ) - xθ
+    val = -float(xθ)
     # xlogy(k - 1, xθ) - xθ -> -∞ for xθ -> ∞ so we only add the first term
     # when it's safe
     if isfinite(xθ)
@@ -31,6 +37,8 @@ function gammalogpdf(k::T, θ::T, x::T) where {T <: Real}
     return x < 0 ? oftype(val, -Inf) : val
 end
 
+gammalogulikelihood(k::Real, θ::Real, x::Real) = gammalogpdf(k, θ, x)
+
 function gammacdf(k::T, θ::T, x::T) where {T <: Real}
     # Handle the degenerate case
     if iszero(k)

src/distrs/hyper.jl

@@ -12,3 +12,7 @@ using .RFunctions:
     hyperinvccdf,
     hyperinvlogcdf,
     hyperinvlogccdf
+
+
+hyperlogupdf(ms::Real, mf::Real, n::Real, x::Real) = hyperlogpdf(ms, mf, n, x)
+hyperlogulikelihood(ms::Real, mf::Real, n::Real, x::Real) = hyperlogpdf(ms, mf, n, x)

src/distrs/nbeta.jl

@@ -12,3 +12,6 @@ using .RFunctions:
     nbetainvccdf,
     nbetainvlogcdf,
     nbetainvlogccdf
+
+nbetalogupdf(α::Real, β::Real, λ::Real, x::Real) = nbetalogpdf(α, β, λ, x)
+nbetalogulikelihood(α::Real, β::Real, λ::Real, x::Real) = nbetalogpdf(α, β, λ, x)

src/distrs/nbinom.jl

@@ -12,3 +12,6 @@ using .RFunctions:
     nbinominvccdf,
     nbinominvlogcdf,
     nbinominvlogccdf
+
+nbinomlogupdf(r::Real, p::Real, x::Real) = nbinomlogpdf(r, p, x)
+nbinomlogulikelihood(r::Real, p::Real, x::Real) = nbinomlogpdf(r, p, x)