Adding unnormalized logpdf functions #185
Description
JuliaStats/Distributions.jl#1978 will add logupdf and logulikelihood functions to the API: logupdf includes all terms in logpdf that depend on the argument (parameter-only terms are excluded), while logulikelihood includes all terms in logpdf that depend on the parameters (argument-only terms are excluded). This more efficient function evaluations and is especially useful when the omitted constants are expensive.
I propose StatsFuns gets a matching function like logupdf. e.g. if a function <distr>logpdf is implemented, so is a function <distr>logupdf. In some cases <distr>logpdf could just call the logupdf version and add a normalization factor.
Here are a few example implementations:
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 = -xθ
# xlogy(k - 1, xθ) - xθ -> -∞ for xθ -> ∞ so we only add the first term
# when it's safe
if isfinite(xθ)
val += xlogy(k - 1, xθ)
end
return x < 0 ? oftype(val, -Inf) : val
end
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)
return x < 0 || x > 1 ? oftype(val, -Inf) : val
end
normlogupdf(z::Number) = abs2(z) / -2
function normlogupdf(μ::Real, σ::Real, x::Number)
if iszero(σ) && x == μ
z = zval(μ, one(σ), x)
else
z = zval(μ, σ, x)
end
return normlogupdf(z)
end
function poislogupdf(λ::T, x::T) where {T <: Real}
val = xlogy(x, λ) - loggamma(x + 1)
return x >= 0 && isinteger(x) ? val : oftype(val, -Inf)
end
function tdistlogupdf(ν::T, x::T) where {T <: Real}
isinf(ν) && return normlogupdf(x)
νp12 = (ν + 1) / 2
return -νp12 * log1p(x^2 / ν)
end