diff --git a/src/matrixvariates.jl b/src/matrixvariates.jl
index c7ac4ed4d2..a5bdd719dd 100644
--- a/src/matrixvariates.jl
+++ b/src/matrixvariates.jl
@@ -226,11 +226,17 @@ rows and `size(d, 2)` columns, or an array of matrices of size `size(d)`.
 """
 loglikelihood(d::MatrixDistribution, X::AbstractMatrix{<:Real}) = logpdf(d, X)
 function loglikelihood(d::MatrixDistribution, X::AbstractArray{<:Real,3})
-    (size(X, 1), size(X, 2)) == size(d) || throw(DimensionMismatch("Inconsistent array dimensions."))
-    return sum(i -> _logpdf(d, view(X, :, :, i)), axes(X, 3))
+    (size(X, 1), size(X, 2)) == size(d) || throw(DimensionMismatch("inconsistent array dimensions"))
+    # we use pairwise summation (https://github.com/JuliaLang/julia/pull/31020)
+    # to compute `sum(Base.Fix1(_logpdf, d), eachslice(X; dims=3))`
+    broadcasted = Broadcast.broadcasted(_logpdf, (d,), (view(X, :, :, i) for i in axes(X, 3)))
+    return sum(Broadcast.instantiate(broadcasted))
 end
 function loglikelihood(d::MatrixDistribution, X::AbstractArray{<:AbstractMatrix{<:Real}})
-    return sum(x -> logpdf(d, x), X)
+    # we use pairwise summation (https://github.com/JuliaLang/julia/pull/31020)
+    # to compute `sum(Base.Fix1(logpdf, d), X)`
+    broadcasted = Broadcast.broadcasted(logpdf, (d,), X)
+    return sum(Broadcast.instantiate(broadcasted))
 end
 
 #  for testing
diff --git a/src/multivariates.jl b/src/multivariates.jl
index 0d2f214c9d..8209291d3e 100644
--- a/src/multivariates.jl
+++ b/src/multivariates.jl
@@ -263,11 +263,15 @@ vectors of length `dim(d)`.
 """
 loglikelihood(d::MultivariateDistribution, X::AbstractVector{<:Real}) = logpdf(d, X)
 function loglikelihood(d::MultivariateDistribution, X::AbstractMatrix{<:Real})
-    size(X, 1) == length(d) || throw(DimensionMismatch("Inconsistent array dimensions."))
-    return sum(i -> _logpdf(d, view(X, :, i)), 1:size(X, 2))
+    size(X, 1) == length(d) || throw(DimensionMismatch("inconsistent array dimensions"))
+    # we use pairwise summation (https://github.com/JuliaLang/julia/pull/31020)
+    broadcasted = Broadcast.broadcasted(_logpdf, (d,), (view(X, :, i) for i in axes(X, 2)))
+    return sum(Broadcast.instantiate(broadcasted))
 end
 function loglikelihood(d::MultivariateDistribution, X::AbstractArray{<:AbstractVector})
-    return sum(x -> logpdf(d, x), X)
+    # we use pairwise summation (https://github.com/JuliaLang/julia/pull/31020)
+    broadcasted = Broadcast.broadcasted(logpdf, (d,), X)
+    return sum(Broadcast.instantiate(broadcasted))
 end
 
 ##### Specific distributions #####
diff --git a/src/univariates.jl b/src/univariates.jl
index 52073863b1..f2235fe3c9 100644
--- a/src/univariates.jl
+++ b/src/univariates.jl
@@ -475,7 +475,11 @@ The log-likelihood of distribution `d` with respect to all samples contained in
 
 Here `x` can be a single scalar sample or an array of samples.
 """
-loglikelihood(d::UnivariateDistribution, X::AbstractArray) = sum(x -> logpdf(d, x), X)
+function loglikelihood(d::UnivariateDistribution, X::AbstractArray)
+    # we use pairwise summation (https://github.com/JuliaLang/julia/pull/31020)
+    broadcasted = Broadcast.broadcasted(logpdf, (d,), X)
+    return sum(Broadcast.instantiate(broadcasted))
+end
 loglikelihood(d::UnivariateDistribution, x::Real) = logpdf(d, x)
 
 ### special definitions for distributions with integer-valued support
@@ -550,11 +554,12 @@ function integerunitrange_cdf(d::DiscreteUnivariateDistribution, x::Integer)
     minimum_d, maximum_d = extrema(d)
     isfinite(minimum_d) || isfinite(maximum_d) || error("support is unbounded")
 
+    # we use pairwise summation (https://github.com/JuliaLang/julia/pull/31020)
     result = if isfinite(minimum_d) && !(isfinite(maximum_d) && x >= div(minimum_d + maximum_d, 2))
-        c = sum(Base.Fix1(pdf, d), minimum_d:(max(x, minimum_d)))
+        c = sum(Broadcast.instantiate(Broadcast.broadcasted(pdf, (d,), minimum_d:(max(x, minimum_d)))))
         x < minimum_d ? zero(c) : c
     else
-        c = 1 - sum(Base.Fix1(pdf, d), (min(x + 1, maximum_d)):maximum_d)
+        c = 1 - sum(Broadcast.instantiate(Broadcast.broadcasted(pdf, (d,), min(x + 1, maximum_d):maximum_d)))
         x >= maximum_d ? one(c) : c
     end
 
@@ -565,11 +570,12 @@ function integerunitrange_ccdf(d::DiscreteUnivariateDistribution, x::Integer)
     minimum_d, maximum_d = extrema(d)
     isfinite(minimum_d) || isfinite(maximum_d) || error("support is unbounded")
 
+    # we use pairwise summation (https://github.com/JuliaLang/julia/pull/31020)
     result = if isfinite(minimum_d) && !(isfinite(maximum_d) && x >= div(minimum_d + maximum_d, 2))
-        c = 1 - sum(Base.Fix1(pdf, d), minimum_d:(max(x, minimum_d)))
+        c = 1 - sum(Broadcast.instantiate(Broadcast.broadcasted(pdf, (d,), minimum_d:(max(x, minimum_d)))))
         x < minimum_d ? one(c) : c
     else
-        c = sum(Base.Fix1(pdf, d), (min(x + 1, maximum_d)):maximum_d)
+        c = sum(Broadcast.instantiate(Broadcast.broadcasted(pdf, (d,), min(x + 1, maximum_d):maximum_d)))
         x >= maximum_d ? zero(c) : c
     end