diff --git a/ext/StatsFunsChainRulesCoreExt.jl b/ext/StatsFunsChainRulesCoreExt.jl
index e94e549..b007090 100644
--- a/ext/StatsFunsChainRulesCoreExt.jl
+++ b/ext/StatsFunsChainRulesCoreExt.jl
@@ -13,6 +13,14 @@ 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),
@@ -22,15 +30,30 @@ ChainRulesCore.@scalar_rule(
         ChainRulesCore.NoTangent(),
     ),
 )
+ChainRulesCore.@scalar_rule(
+    binomlogupdf(n::Real, p::Real, k::Real),
+    (
+        ChainRulesCore.NoTangent(),
+        (k / p - n) / (1 - p),
+        ChainRulesCore.NoTangent(),
+    ),
+)
 
 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),
+    (
+        log(x) / 2,
+        (k - x - 2) / (2 * x),
+    ),
+)
 
 ChainRulesCore.@scalar_rule(
     fdistlogpdf(ν1::Real, ν2::Real, x::Number),
@@ -47,6 +70,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,
+        xotmp = x / tmp,
+        a = xotmp * (ν1 + ν2),
+    ),
+    (
+        (log(xotmp) - a) / 2,
+        -(log(tmp) + (ν1 + ν2) / tmp) / 2,
+        (ν1 * (1 - a) - 2) / (2 * x),
+    ),
+)
 
 ChainRulesCore.@scalar_rule(
     gammalogpdf(k::Real, θ::Real, x::Number),
@@ -61,11 +97,27 @@ ChainRulesCore.@scalar_rule(
         - (1 + z) / x,
     ),
 )
+ChainRulesCore.@scalar_rule(
+    gammalogupdf(k::Real, θ::Real, x::Number),
+    (
+        log(x),
+        x / θ^2,
+        (k - 1) / x - inv(θ),
+    ),
+)
 
 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),
@@ -81,5 +133,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
diff --git a/src/StatsFuns.jl b/src/StatsFuns.jl
index 3a39252..4e43ec8 100644
--- a/src/StatsFuns.jl
+++ b/src/StatsFuns.jl
@@ -56,6 +56,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
@@ -68,6 +70,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
@@ -80,6 +84,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
@@ -92,6 +98,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
@@ -104,6 +112,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
@@ -116,6 +126,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
@@ -128,6 +140,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
@@ -140,6 +154,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
@@ -152,6 +168,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
@@ -164,6 +182,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
@@ -176,6 +196,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
@@ -188,6 +210,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
@@ -200,6 +224,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
@@ -212,6 +238,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
@@ -224,6 +252,8 @@ export
     # distrs/signrank
     signrankpdf,
     signranklogpdf,
+    signranklogupdf,
+    signranklogulikelihood,
     signrankcdf,
     signranklogcdf,
     signrankccdf,
@@ -246,6 +276,8 @@ export
     # distrs/wilcox
     wilcoxpdf,
     wilcoxlogpdf,
+    wilcoxlogupdf,
+    wilcoxlogulikelihood,
     wilcoxcdf,
     wilcoxlogcdf,
     wilcoxccdf,
diff --git a/src/distrs/beta.jl b/src/distrs/beta.jl
index 79c6045..f4b9dcc 100644
--- a/src/distrs/beta.jl
+++ b/src/distrs/beta.jl
@@ -7,12 +7,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
diff --git a/src/distrs/binom.jl b/src/distrs/binom.jl
index 7d16c48..8809a6a 100644
--- a/src/distrs/binom.jl
+++ b/src/distrs/binom.jl
@@ -5,11 +5,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"))
diff --git a/src/distrs/chisq.jl b/src/distrs/chisq.jl
index 2204e94..577a88d 100644
--- a/src/distrs/chisq.jl
+++ b/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", "cdf", "ccdf", "logcdf", "logccdf", "invcdf", "invccdf", "invlogcdf", "invlogccdf")
     _chisqf = Symbol("chisq" * f)
     _gammaf = Symbol("gamma" * f)
     @eval begin
@@ -9,3 +9,11 @@ for f in ("pdf", "logpdf", "cdf", "ccdf", "logcdf", "logccdf", "invcdf", "invccd
         $(_chisqf)(k::T, x::T) where {T <: Real} = $(_gammaf)(k / 2, 2, x)
     end
 end
+
+chisqlogulikelihood(k::Real, x::Real) = chisqlogulikelihood(promote(k, x)...)
+function chisqlogulikelihood(k::T, x::T) where {T <: Real}
+    y = max(x, 0)
+    k2 = k / 2
+    val = xlogy(k2, x / 2) - loggamma(k2)
+    return x < 0 ? oftype(val, -Inf) : val
+end
diff --git a/src/distrs/fdist.jl b/src/distrs/fdist.jl
index 856d40f..aad3e31 100644
--- a/src/distrs/fdist.jl
+++ b/src/distrs/fdist.jl
@@ -12,6 +12,27 @@ 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`
+    y = max(x, 0)
+    val = (xlogy(ν1 - 2, y) - xlogy(ν1 + ν2, ν1 * y + ν2)) / 2
+    return x < 0 ? oftype(val, -Inf) : val
+end
+
+fdistlogulikelihood(ν1::Real, ν2::Real, x::Real) = fdistlogulikelihood(promote(ν1, ν2, x)...)
+function fdistlogulikelihood(ν1::T, ν2::T, x::T) where {T <: Real}
+    # we ensure that `log(x)` does not error if `x < 0`
+    y = max(x, 0)
+    tmp = ν1 * y + ν2
+    a = ν1 / tmp
+    b = ν2 / tmp
+    halfν1 = ν1 / 2
+    halfν2 = ν2 / 2
+    val = (xlogy(halfν1, a) + xlogy(halfν2, b)) - logbeta(halfν1, halfν2)
+    return x < 0 ? oftype(val, -Inf) : val
+end
+
 for f in ("cdf", "ccdf", "logcdf", "logccdf")
     ff = Symbol("fdist" * f)
     bf = Symbol("beta" * f)
diff --git a/src/distrs/gamma.jl b/src/distrs/gamma.jl
index cee2fca..064dc33 100644
--- a/src/distrs/gamma.jl
+++ b/src/distrs/gamma.jl
@@ -7,13 +7,32 @@ 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) - xlogy(k, θ) : 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θ
+    y = max(x, 0)
+    val = -float(y / θ)
     # 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θ)
+    if isfinite(val)
+        val += xlogy(k - 1, y)
+    end
+    return x < 0 ? oftype(val, -Inf) : val
+end
+
+gammalogulikelihood(k::Real, θ::Real, x::Real) = gammalogulikelihood(promote(k, θ, x)...)
+function gammalogulikelihood(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θ - loggamma(k)
+    # xlogy(k, xθ) - xθ -> -∞ for xθ -> ∞ so we only add the first term
+    # when it's safe
+    if isfinite(val)
+        val += xlogy(k, xθ)
     end
     return x < 0 ? oftype(val, -Inf) : val
 end
diff --git a/src/distrs/hyper.jl b/src/distrs/hyper.jl
index e7cf1c7..a6b0f70 100644
--- a/src/distrs/hyper.jl
+++ b/src/distrs/hyper.jl
@@ -1,5 +1,8 @@
 # functions related to hyper-geometric distribution
 
+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)
+
 # Rmath implementations
 function hyperpdf(ms::Real, mf::Real, n::Real, x::Real)
     T = float(Base.promote_typeof(ms, mf, n, x))
diff --git a/src/distrs/nbeta.jl b/src/distrs/nbeta.jl
index 9ed99a7..0d7d82e 100644
--- a/src/distrs/nbeta.jl
+++ b/src/distrs/nbeta.jl
@@ -1,5 +1,8 @@
 # functions related to noncentral beta distribution
 
+nbetalogupdf(α::Real, β::Real, λ::Real, x::Real) = nbetalogpdf(α, β, λ, x)
+nbetalogulikelihood(α::Real, β::Real, λ::Real, x::Real) = nbetalogpdf(α, β, λ, x)
+
 # Rmath implementations
 function nbetapdf(α::Real, β::Real, λ::Real, x::Real)
     T = float(Base.promote_typeof(α, β, λ, x))
diff --git a/src/distrs/nbinom.jl b/src/distrs/nbinom.jl
index 580e5db..7b785b8 100644
--- a/src/distrs/nbinom.jl
+++ b/src/distrs/nbinom.jl
@@ -12,10 +12,16 @@ nbinompdf(r::Real, p::Real, k::Real) = exp(nbinomlogpdf(r, p, k))
 
 nbinomlogpdf(r::Real, p::Real, k::Real) = nbinomlogpdf(promote(r, p, k)...)
 function nbinomlogpdf(r::T, p::T, k::T) where {T <: Real}
+    logupdf = nbinomlogupdf(r, p, k)
+    return isfinite(logupdf) ? logupdf + xlogy(r, p) : logupdf
+end
+
+nbinomlogupdf(r::Real, p::Real, k::Real) = nbinomlogupdf(promote(r, p, k)...)
+function nbinomlogupdf(r::T, p::T, k::T) where {T <: Real}
     if !(0 <= p <= 1) || r <= 0
         return float(T)(NaN)
     elseif k >= 0 && isinteger(k)
-        z = xlog1py(k, -p) + xlogy(r, p)
+        z = xlog1py(k, -p)
         iszero(k) && return z
         res = z - (logbeta(k + 1, r) + log(r + k))
         if !isnan(res)
@@ -27,6 +33,9 @@ function nbinomlogpdf(r::T, p::T, k::T) where {T <: Real}
         return float(T)(isnan(k) ? NaN : -Inf)
     end
 end
+
+nbinomlogulikelihood(r::Real, p::Real, k::Real) = nbinomlogpdf(r, p, k)
+
 function nbinomcdf(r::Real, p::Real, k::Real)
     if k < 0
         return zero(float(first(promote(r, p, k))))
diff --git a/src/distrs/nchisq.jl b/src/distrs/nchisq.jl
index c1f8aea..b6045a0 100644
--- a/src/distrs/nchisq.jl
+++ b/src/distrs/nchisq.jl
@@ -1,5 +1,8 @@
 # functions related to noncentral chi-square distribution
 
+nchisqlogupdf(k::Real, λ::Real, x::Real) = nchisqlogpdf(k, λ, x)
+nchisqlogulikelihood(k::Real, λ::Real, x::Real) = nchisqlogpdf(k, λ, x)
+
 # Rmath implementations
 function nchisqpdf(k::Real, λ::Real, x::Real)
     T = float(Base.promote_typeof(k, λ, x))
diff --git a/src/distrs/nfdist.jl b/src/distrs/nfdist.jl
index a6c564e..ad0adb6 100644
--- a/src/distrs/nfdist.jl
+++ b/src/distrs/nfdist.jl
@@ -1,5 +1,8 @@
 # functions related to noncentral F distribution
 
+nfdistlogupdf(k1::Real, k2::Real, λ::Real, x::Real) = nfdistlogpdf(k1, k2, λ, x)
+nfdistlogulikelihood(k1::Real, k2::Real, λ::Real, x::Real) = nfdistlogpdf(k1, k2, λ, x)
+
 # Rmath implementations
 function nfdistpdf(k1::Real, k2::Real, λ::Real, x::Real)
     T = float(Base.promote_typeof(k1, k2, λ, x))
diff --git a/src/distrs/norm.jl b/src/distrs/norm.jl
index e79713c..204dd65 100644
--- a/src/distrs/norm.jl
+++ b/src/distrs/norm.jl
@@ -41,6 +41,32 @@ function normlogpdf(μ::Real, σ::Real, x::Number)
     return normlogpdf(z) - log(σ)
 end
 
+# logupdf
+normlogupdf(z::Number) = -abs2(z) / 2
+function normlogupdf(μ::Real, σ::Real, x::Number)
+    return if iszero(σ) && x == μ
+        oftype(normlogupdf(zval(μ, oneunit(σ), x)), Inf)
+    else
+        normlogupdf(zval(μ, σ, x))
+    end
+end
+
+# logulikelihood
+normlogulikelihood(z::Number) = normlogupdf(z)
+function normlogulikelihood(μ::Real, σ::Real, x::Number)
+    if iszero(σ)
+        if x == μ
+            z = zval(μ, oneunit(σ), x)
+        else
+            z = zval(μ, σ, x)
+            σ = oneunit(σ)
+        end
+    else
+        z = zval(μ, σ, x)
+    end
+    return normlogulikelihood(z) - log(σ)
+end
+
 # cdf
 normcdf(z::Number) = erfc(-z * invsqrt2) / 2
 function normcdf(μ::Real, σ::Real, x::Number)
diff --git a/src/distrs/ntdist.jl b/src/distrs/ntdist.jl
index 8f50540..7e690b8 100644
--- a/src/distrs/ntdist.jl
+++ b/src/distrs/ntdist.jl
@@ -1,5 +1,8 @@
 # functions related to noncentral T distribution
 
+ntdistlogupdf(k::Real, λ::Real, x::Real) = ntdistlogpdf(k, λ, x)
+ntdistlogulikelihood(k::Real, λ::Real, x::Real) = ntdistlogpdf(k, λ, x)
+
 # Rmath implementations
 function ntdistpdf(k::Real, λ::Real, x::Real)
     T = float(Base.promote_typeof(k, λ, x))
diff --git a/src/distrs/pois.jl b/src/distrs/pois.jl
index 1c9b565..a6bdb7a 100644
--- a/src/distrs/pois.jl
+++ b/src/distrs/pois.jl
@@ -3,9 +3,19 @@
 # Julia implementations
 poispdf(λ::Real, x::Real) = exp(poislogpdf(λ, x))
 
-poislogpdf(λ::Real, x::Real) = poislogpdf(promote(λ, x)...)
-function poislogpdf(λ::T, x::T) where {T <: Real}
-    val = xlogy(x, λ) - λ - loggamma(x + 1)
+function poislogpdf(λ::Real, x::Real)
+    logupdf = poislogupdf(λ, x)
+    return isfinite(logupdf) ? logupdf - λ : logupdf
+end
+
+poislogupdf(λ::Real, x::Real) = poislogupdf(promote(λ, x)...)
+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 poislogulikelihood(λ::Real, x::Real)
+    val = xlogy(x, λ) - λ
     return x >= 0 && isinteger(x) ? val : oftype(val, -Inf)
 end
 
diff --git a/src/distrs/signrank.jl b/src/distrs/signrank.jl
index 96936db..64b3224 100644
--- a/src/distrs/signrank.jl
+++ b/src/distrs/signrank.jl
@@ -25,25 +25,36 @@ the number of ways {1,2,...,j} can sum to W-i+1.
     return DP
 end
 
-function signrankpdf(n::Int, W::Float64)
-    return isinteger(W) ? signrankpdf(n, Int(W)) : 0.0
+function signrankpdf(n::Int, W::Union{Int, Float64})
+    numsets = signrank_numsets(n, W)
+    return iszero(numsets) ? 0.0 : ldexp(float(numsets), -n)
 end
-function signrankpdf(n::Int, W::Int)
+
+function signrank_numsets(n::Int, W::Float64)
+    return isinteger(W) ? signrank_numsets(n, Int(W)) : 0
+end
+function signrank_numsets(n::Int, W::Int)
     if W < 0
-        return 0.0
+        return 0
     end
     max_W = (n * (n + 1)) >> 1
     W2 = max_W - W
     if W2 < W
-        return signrankpdf(n, W2)
+        return signrank_numsets(n, W2)
     end
     DP = signrankDP(n, W)
-    return ldexp(float(DP[1]), -n)
+    return DP[1]
 end
 
 function signranklogpdf(n::Int, W::Union{Float64, Int})
     return log(signrankpdf(n, W))
 end
+function signranklogupdf(n::Int, W::Union{Float64, Int})
+    return log(signrank_numsets(n, W))
+end
+function signranklogulikelihood(n::Int, W::Union{Float64, Int})
+    return signranklogpdf(n, W)
+end
 
 function signrankcdf(n::Int, W::Float64)
     return signrankcdf(n, round(Int, W, RoundNearestTiesUp))
diff --git a/src/distrs/tdist.jl b/src/distrs/tdist.jl
index 8b1fef9..9f2e180 100644
--- a/src/distrs/tdist.jl
+++ b/src/distrs/tdist.jl
@@ -9,6 +9,20 @@ function tdistlogpdf(ν::T, x::T) where {T <: Real}
     return loggamma(νp12) - (logπ + log(ν)) / 2 - loggamma(ν / 2) - νp12 * log1p(x^2 / ν)
 end
 
+tdistlogupdf(ν::Real, x::Real) = tdistlogupdf(promote(ν, x)...)
+function tdistlogupdf(ν::T, x::T) where {T <: Real}
+    isinf(ν) && return normlogupdf(x)
+    νp12 = (ν + 1) / 2
+    return - νp12 * log1p(x^2 / ν)
+end
+
+tdistlogulikelihood(ν::Real, x::Real) = tdistlogulikelihood(promote(ν, x)...)
+function tdistlogulikelihood(ν::T, x::T) where {T <: Real}
+    isinf(ν) && return normlogulikelihood(x)
+    νp12 = (ν + 1) / 2
+    return loggamma(νp12) - log(ν) / 2 - loggamma(ν / 2) - νp12 * log1p(x^2 / ν)
+end
+
 function tdistcdf(ν::T, x::T) where {T <: Real}
     if isinf(ν)
         return normcdf(x)
diff --git a/src/distrs/wilcox.jl b/src/distrs/wilcox.jl
index 64511cb..d54d3cc 100644
--- a/src/distrs/wilcox.jl
+++ b/src/distrs/wilcox.jl
@@ -74,22 +74,33 @@ A. Löffler: "Über eine Partition der nat. Zahlen und ihre Anwendung beim U-Tes
     return partitions
 end
 
-function wilcoxpdf(nx::Int, ny::Int, U::Float64)
-    return isinteger(U) ? wilcoxpdf(nx, ny, Int(U)) : 0.0
+function wilcoxpdf(nx::Int, ny::Int, U::Union{Float64, Int})
+    numseqs = wilcox_numseqs(nx, ny, U)
+    return iszero(numseqs) ? 0.0 : numseqs / binomial(nx + ny, nx)
 end
-function wilcoxpdf(nx::Int, ny::Int, U::Int)
+
+function wilcox_numseqs(nx::Int, ny::Int, U::Float64)
+    return isinteger(U) ? wilcox_numseqs(nx, ny, Int(U)) : 0
+end
+function wilcox_numseqs(nx::Int, ny::Int, U::Int)
     max_U = nx * ny
     if !(0 <= U <= max_U)
-        return 0.0
+        return 0
     end
     U = min(U, max_U - U)
     partitions = wilcox_partitions(nx, ny, U)
-    return partitions[end] / binomial(nx + ny, nx)
+    return partitions[end]
 end
 
 function wilcoxlogpdf(nx::Int, ny::Int, U::Union{Float64, Int})
     return log(wilcoxpdf(nx, ny, U))
 end
+function wilcoxlogupdf(nx::Int, ny::Int, U::Union{Float64, Int})
+    return log(wilcox_numseqs(nx, ny, U))
+end
+function wilcoxlogulikelihood(nx::Int, ny::Int, U::Union{Float64, Int})
+    return wilcoxlogpdf(nx, ny, U)
+end
 
 function wilcoxcdf(nx::Int, ny::Int, U::Float64)
     return wilcoxcdf(nx, ny, round(Int, U, RoundNearestTiesUp))
diff --git a/test/chainrules.jl b/test/chainrules.jl
index d59b5ea..adf618c 100644
--- a/test/chainrules.jl
+++ b/test/chainrules.jl
@@ -9,45 +9,63 @@ using Random
     z = logistic(randn())
     test_frule(betalogpdf, x, y, z)
     test_rrule(betalogpdf, x, y, z)
+    test_frule(betalogupdf, x, y, z)
+    test_rrule(betalogupdf, x, y, z)
 
     x = exp(randn())
     y = exp(randn())
     z = exp(randn())
     test_frule(gammalogpdf, x, y, z)
     test_rrule(gammalogpdf, x, y, z)
+    test_frule(gammalogupdf, x, y, z)
+    test_rrule(gammalogupdf, x, y, z)
 
     x = exp(randn())
     y = exp(randn())
     test_frule(chisqlogpdf, x, y)
     test_rrule(chisqlogpdf, x, y)
+    test_frule(chisqlogupdf, x, y)
+    test_rrule(chisqlogupdf, x, y)
 
     x = exp(randn())
     y = exp(randn())
     z = exp(randn())
     test_frule(fdistlogpdf, x, y, z)
     test_rrule(fdistlogpdf, x, y, z)
+    test_frule(fdistlogupdf, x, y, z)
+    test_rrule(fdistlogupdf, x, y, z)
 
     x = exp(randn())
     y = randn()
     test_frule(tdistlogpdf, x, y)
     test_rrule(tdistlogpdf, x, y)
+    test_frule(tdistlogupdf, x, y)
+    test_rrule(tdistlogupdf, x, y)
 
     x = rand(1:100)
     y = logistic(randn())
     z = rand(1:x)
     test_frule(binomlogpdf, x, y, z)
     test_rrule(binomlogpdf, x, y, z)
+    test_frule(binomlogupdf, x, y, z)
+    test_rrule(binomlogupdf, x, y, z)
 
     x = exp(randn())
     y = rand(1:100)
     test_frule(poislogpdf, x, y)
     test_rrule(poislogpdf, x, y)
+    test_frule(poislogupdf, x, y)
+    test_rrule(poislogupdf, x, y)
+    test_frule(poislogulikelihood, x, y)
+    test_rrule(poislogulikelihood, x, y)
 
     # test special case λ = 0
-    _, pb = rrule(poislogpdf, 0.0, 0)
-    _, x̄1, _ = pb(1)
-    @test x̄1 == -1
-    _, pb = rrule(poislogpdf, 0.0, 1)
-    _, x̄1, _ = pb(1)
-    @test x̄1 == Inf
+    for f in (poislogpdf, poislogupdf, poislogulikelihood)
+        _, pb = rrule(f, 0.0, 0)
+        _, x̄1, _ = pb(1)
+        @test x̄1 == (f === poislogupdf ? 0 : -1)
+        _, pb = rrule(f, 0.0, 1)
+        _, x̄1, _ = pb(1)
+        @test x̄1 == Inf
+    end
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index 0cd38a1..aeed0cd 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,4 +1,4 @@
-tests = ["rmath", "generic", "misc", "chainrules", "inverse", "tvpack", "qa"]
+tests = ["rmath", "generic", "misc", "chainrules", "inverse", "tvpack", "unnormalized", "qa"]
 
 for t in tests
     fp = "$t.jl"
diff --git a/test/unnormalized.jl b/test/unnormalized.jl
new file mode 100644
index 0000000..4bb3fad
--- /dev/null
+++ b/test/unnormalized.jl
@@ -0,0 +1,115 @@
+using StatsFuns
+using SpecialFunctions
+using Test
+
+@testset "logupdf + logulikelihood" begin
+    @testset "optimized" begin
+        # Beta distribution
+        for ((α, β, x), T) in (((0.1, 2.5, 0.7), Float64), ((1, 2.5, 0.7), Float64), ((0.1f0, 2.5f0, 0.7f0), Float32), ((1, 2.5f0, 0.7f0), Float32))
+            @test (@inferred betalogupdf(α, β, x))::T ≈ betalogpdf(α, β, x) + logbeta(α, β)
+            @test (@inferred betalogulikelihood(α, β, x))::T == betalogpdf(α, β, x)
+        end
+
+        # Binomial distribution
+        for ((n, p, k), T) in (((9, 0.4, 6), Float64), ((9.0f0, 0.4f0, 6.0f0), Float32), ((9, 0.4f0, 6), Float32))
+            @test (@inferred binomlogupdf(n, p, k))::T ≈ binomlogpdf(n, p, k) + log(T(n) + 1)
+            @test (@inferred binomlogulikelihood(n, p, k))::T == binomlogpdf(n, p, k)
+        end
+
+        # Chi-squared distribution
+        for ((k, x), T) in (((4.2, 3.1), Float64), ((4, 3.1), Float64), ((4.2f0, 3.1f0), Float32), ((4, 3.1f0), Float32))
+            @test (@inferred chisqlogupdf(k, x))::T ≈ chisqlogpdf(k, x) + T(k) / 2 * log(T(2)) + loggamma(T(k) / 2)
+            @test (@inferred chisqlogulikelihood(k, x))::T ≈ chisqlogpdf(k, x) + log(T(x)) + T(x) / 2
+        end
+
+        # F distribution
+        for ((ν1, ν2, x), T) in (((0.9, 1.5, 2.1), Float64), ((1, 2, 2.1), Float64), ((0.9f0, 1.5f0, 2.1f0), Float32), ((1, 2, 2.1f0), Float32))
+            @test (@inferred fdistlogupdf(ν1, ν2, x))::T ≈ fdistlogpdf(ν1, ν2, x) + logbeta(T(ν1) / 2, T(ν2) / 2) - (T(ν1) * log(T(ν1)) + T(ν2) * log(T(ν2))) / 2
+            @test (@inferred fdistlogulikelihood(ν1, ν2, x))::T ≈ fdistlogpdf(ν1, ν2, x) - (T(ν1) / 2 - 1) * log(T(x))
+        end
+
+        # Gamma distribution
+        for ((k, θ, x), T) in (((1.4, 2.3, 1.9), Float64), ((2, 2.3, 1.9), Float64), ((1.4f0, 2.3f0, 1.9f0), Float32), ((2, 2.3f0, 1.9f0), Float32))
+            @test (@inferred gammalogupdf(k, θ, x))::T ≈ gammalogpdf(k, θ, x) + loggamma(T(k)) + T(k) * log(T(θ))
+            @test (@inferred gammalogulikelihood(k, θ, x))::T ≈ gammalogpdf(k, θ, x) + log(T(x))
+        end
+
+        # Negative binomial distribution
+        for ((r, p, x), T) in (((3, 0.7, 2), Float64), ((3.0f0, 0.7f0, 2.0f0), Float32), ((3, 0.7f0, 2), Float32))
+            @test (@inferred nbinomlogupdf(r, p, x))::T ≈ nbinomlogpdf(r, p, x) - xlogy(T(r), T(p))
+            @test (@inferred nbinomlogulikelihood(r, p, x))::T == nbinomlogpdf(r, p, x)
+        end
+        for ((r, p, x), T) in (((3, 1.5, 2), Float64), ((-1, 0.5, 2), Float64), ((3.0f0, 1.5f0, 2.0f0), Float32))
+            @test isnan((@inferred nbinomlogupdf(r, p, x))::T)
+        end
+
+        # Normal distribution
+        for (z, T) in ((3.1, Float64), (3, Float64), (3.1f0, Float32))
+            @test (@inferred normlogupdf(z))::T ≈ normlogpdf(z) + log(T(2π)) / 2
+            @test (@inferred normlogulikelihood(z))::T ≈ normlogpdf(z) + log(T(2π)) / 2
+        end
+        for ((μ, σ, z), T) in (
+                ((-0.1, 2.4, 3.1), Float64), ((-0.1, 0.0, 3.1), Float64), ((3.1, 0.0, 3.1), Float64),
+                ((-0.1f0, 2.4f0, 3.1f0), Float32), ((-0.1f0, 0.0f0, 3.1f0), Float32), ((3.1f0, 0.0f0, 3.1f0), Float32), ((0, 2.4f0, 3.1f0), Float32),
+            )
+            @test (@inferred normlogupdf(μ, σ, z))::T ≈ (iszero(σ) ? normlogpdf(μ, σ, z) : normlogpdf(μ, σ, z) + log(T(2π)) / 2 + log(σ))
+            @test (@inferred normlogulikelihood(μ, σ, z))::T ≈ normlogpdf(μ, σ, z) + log(T(2π)) / 2
+        end
+
+        # Poisson distribution
+        for ((λ, x), T) in (((3.5, 2), Float64), ((3.5f0, 2.0f0), Float32), ((3.5f0, 2), Float32))
+            @test (@inferred poislogupdf(λ, x))::T ≈ poislogpdf(λ, x) + T(λ)
+            @test (@inferred poislogulikelihood(λ, x))::T ≈ poislogpdf(λ, x) + loggamma(T(x) + 1)
+        end
+
+        # Student's t distribution
+        for ((ν, x), T) in (((4.5, 1.3), Float64), ((5, 1.3), Float64), ((4.5f0, 1.3f0), Float32), ((5, 1.3f0), Float32))
+            @test (@inferred tdistlogupdf(ν, x))::T ≈ tdistlogpdf(ν, x) - loggamma((T(ν) + 1) / 2) + (log(T(π)) + log(T(ν))) / 2 + loggamma(T(ν) / 2)
+            @test (@inferred tdistlogulikelihood(ν, x))::T ≈ tdistlogpdf(ν, x) + log(T(π)) / 2
+        end
+
+        # Wilcoxon signed rank distribution
+        for (n, W) in ((5, 7), (5, 7.0))
+            @test (@inferred signranklogupdf(n, W))::Float64 ≈ signranklogpdf(n, W) + n * log(2)
+            @test (@inferred signranklogulikelihood(n, W))::Float64 == signranklogpdf(n, W)
+        end
+
+        # Wilcoxon rank sum distribution
+        for (nx, ny, U) in ((3, 4, 5), (3, 4, 5.0))
+            @test (@inferred wilcoxlogupdf(nx, ny, U))::Float64 ≈ wilcoxlogpdf(nx, ny, U) + first(logabsbinomial(nx + ny, nx))
+            @test (@inferred wilcoxlogulikelihood(nx, ny, U))::Float64 == wilcoxlogpdf(nx, ny, U)
+        end
+    end
+
+    @testset "fallback" begin
+        # Hyper-geometric distribution
+        for (ms, mf, n, x) in ((2, 3, 4, 2), (2, 3, 4, 2.0))
+            @test (@inferred hyperlogupdf(ms, mf, n, x))::Float64 == hyperlogpdf(ms, mf, n, x)
+            @test (@inferred hyperlogulikelihood(ms, mf, n, x))::Float64 == hyperlogpdf(ms, mf, n, x)
+        end
+
+        # Non-central beta distribution
+        for ((α, β, λ, x), T) in (((0.8, 2.1, 1.1, 0.8), Float64), ((1, 2.1, 1.1, 0.8), Float64), ((0.8f0, 2.1f0, 1.1f0, 0.8f0), Float32), ((1, 2.1f0, 1.1f0, 0.8f0), Float32))
+            @test (@inferred nbetalogupdf(α, β, λ, x))::T == nbetalogpdf(α, β, λ, x)
+            @test (@inferred nbetalogulikelihood(α, β, λ, x))::T == nbetalogpdf(α, β, λ, x)
+        end
+
+        # Non-central chi-squared distribution
+        for ((k, λ, x), T) in (((3.0, 1.5, 4.2), Float64), ((3, 1.5, 4.2), Float64), ((3.0f0, 1.5f0, 4.2f0), Float32), ((3, 1.5f0, 4.2f0), Float32))
+            @test (@inferred nchisqlogupdf(k, λ, x))::T == nchisqlogpdf(k, λ, x)
+            @test (@inferred nchisqlogulikelihood(k, λ, x))::T == nchisqlogpdf(k, λ, x)
+        end
+
+        # Non-central F distribution
+        for ((k1, k2, λ, x), T) in (((2.0, 3.0, 1.0, 1.5), Float64), ((2, 3, 1.0, 1.5), Float64), ((2.0f0, 3.0f0, 1.0f0, 1.5f0), Float32), ((2, 3, 1.0f0, 1.5f0), Float32))
+            @test (@inferred nfdistlogupdf(k1, k2, λ, x))::T == nfdistlogpdf(k1, k2, λ, x)
+            @test (@inferred nfdistlogulikelihood(k1, k2, λ, x))::T == nfdistlogpdf(k1, k2, λ, x)
+        end
+
+        # Non-central t distribution
+        for ((k, λ, x), T) in (((5.0, 1.2, 2.0), Float64), ((5, 1.2, 2.0), Float64), ((5.0f0, 1.2f0, 2.0f0), Float32), ((5, 1.2f0, 2.0f0), Float32))
+            @test (@inferred ntdistlogupdf(k, λ, x))::T == ntdistlogpdf(k, λ, x)
+            @test (@inferred ntdistlogulikelihood(k, λ, x))::T == ntdistlogpdf(k, λ, x)
+        end
+    end
+end