diff --git a/Project.toml b/Project.toml
index 0ce5c9a..0a381fe 100644
--- a/Project.toml
+++ b/Project.toml
@@ -23,6 +23,9 @@ FastGaussQuadrature = "0.4, 0.5, 1"
 Functors = "0.1, 0.2, 0.3, 0.4"
 InverseFunctions = "0.1.2"
 IrrationalConstants = "0.1, 0.2"
+LinearAlgebra = "1"
+Random = "1"
 SpecialFunctions = "1, 2"
 StatsFuns = "0.9.13, 1"
+Test = "1"
 julia = "1.6"
diff --git a/README.md b/README.md
index d0e86a3..d97b9c9 100644
--- a/README.md
+++ b/README.md
@@ -6,6 +6,7 @@
 [![Codecov](https://codecov.io/gh/JuliaGaussianProcesses/GPLikelihoods.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/JuliaGaussianProcesses/GPLikelihoods.jl)
 [![Code Style: Blue](https://img.shields.io/badge/code%20style-blue-4495d1.svg)](https://github.com/invenia/BlueStyle)
 [![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor's%20Guide-blueviolet)](https://github.com/SciML/ColPrac)
+[![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
 
 GPLikelihoods.jl provides a collection of likelihoods to be used as building
 blocks for defining non-Gaussian problems.  It is intended to be mainly
diff --git a/src/expectations.jl b/src/expectations.jl
index 85db8c3..3354548 100644
--- a/src/expectations.jl
+++ b/src/expectations.jl
@@ -1,8 +1,10 @@
-using FastGaussQuadrature: gausshermite
+using FastGaussQuadrature: FastGaussQuadrature
 using SpecialFunctions: loggamma
 using ChainRulesCore: ChainRulesCore
 using IrrationalConstants: sqrt2, invsqrtπ
 
+gausshermite(n::Integer) = FastGaussQuadrature.gausshermite(n)
+
 struct DefaultExpectationMethod end
 
 struct AnalyticExpectation end
diff --git a/src/likelihoods/bernoulli.jl b/src/likelihoods/bernoulli.jl
index d14878e..e8bc08b 100644
--- a/src/likelihoods/bernoulli.jl
+++ b/src/likelihoods/bernoulli.jl
@@ -18,4 +18,4 @@ BernoulliLikelihood(l=logistic) = BernoulliLikelihood(link(l))
 
 (l::BernoulliLikelihood)(f::Real) = Bernoulli(l.invlink(f))
 
-(l::BernoulliLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs))
+(l::BernoulliLikelihood)(fs::AbstractVector{<:Real}) = product_distribution(map(l, fs))
diff --git a/src/likelihoods/categorical.jl b/src/likelihoods/categorical.jl
index 6caedbf..235310c 100644
--- a/src/likelihoods/categorical.jl
+++ b/src/likelihoods/categorical.jl
@@ -36,4 +36,4 @@ function (l::CategoricalLikelihood)(f::AbstractVector{<:Real})
     return Categorical(l.invlink(f))
 end
 
-(l::CategoricalLikelihood)(fs::AbstractVector) = Product(map(l, fs))
+(l::CategoricalLikelihood)(fs::AbstractVector) = product_distribution(map(l, fs))
diff --git a/src/likelihoods/exponential.jl b/src/likelihoods/exponential.jl
index 1b7dc60..9ab648e 100644
--- a/src/likelihoods/exponential.jl
+++ b/src/likelihoods/exponential.jl
@@ -15,7 +15,7 @@ ExponentialLikelihood(l=exp) = ExponentialLikelihood(link(l))
 
 (l::ExponentialLikelihood)(f::Real) = Exponential(l.invlink(f))
 
-(l::ExponentialLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs))
+(l::ExponentialLikelihood)(fs::AbstractVector{<:Real}) = product_distribution(map(l, fs))
 
 function expected_loglikelihood(
     ::AnalyticExpectation,
diff --git a/src/likelihoods/gamma.jl b/src/likelihoods/gamma.jl
index dcd0202..61bd4cb 100644
--- a/src/likelihoods/gamma.jl
+++ b/src/likelihoods/gamma.jl
@@ -20,7 +20,7 @@ GammaLikelihood(α::Real=1.0, l=exp) = GammaLikelihood(α, link(l))
 
 (l::GammaLikelihood)(f::Real) = Gamma(l.α, l.invlink(f))
 
-(l::GammaLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs))
+(l::GammaLikelihood)(fs::AbstractVector{<:Real}) = product_distribution(map(l, fs))
 
 function expected_loglikelihood(
     ::AnalyticExpectation,
diff --git a/src/likelihoods/negativebinomial.jl b/src/likelihoods/negativebinomial.jl
index d8f78a2..03c42d8 100644
--- a/src/likelihoods/negativebinomial.jl
+++ b/src/likelihoods/negativebinomial.jl
@@ -79,7 +79,9 @@ end
 
 @functor NegativeBinomialLikelihood
 
-(l::NegativeBinomialLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs))
+function (l::NegativeBinomialLikelihood)(fs::AbstractVector{<:Real})
+    return product_distribution(map(l, fs))
+end
 
 @doc raw"""
     NBParamSuccess(successes)
diff --git a/src/likelihoods/poisson.jl b/src/likelihoods/poisson.jl
index 9e31d6e..300f120 100644
--- a/src/likelihoods/poisson.jl
+++ b/src/likelihoods/poisson.jl
@@ -18,7 +18,7 @@ PoissonLikelihood(l=exp) = PoissonLikelihood(link(l))
 
 (l::PoissonLikelihood)(f::Real) = Poisson(l.invlink(f))
 
-(l::PoissonLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs))
+(l::PoissonLikelihood)(fs::AbstractVector{<:Real}) = product_distribution(map(l, fs))
 
 function expected_loglikelihood(
     ::AnalyticExpectation,
diff --git a/test/Project.toml b/test/Project.toml
index 6e2c775..feeb76e 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Functors = "d9f16b24-f501-4c13-a1f2-28368ffc5196"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -7,6 +8,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
+Aqua = "0.8"
 Distributions = "0.19, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25"
 Functors = "0.1, 0.2, 0.3, 0.4"
 StatsFuns = "0.9, 1"
diff --git a/test/runtests.jl b/test/runtests.jl
index 1552a99..05528af 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,11 +1,13 @@
 using GPLikelihoods
 using GPLikelihoods: GaussHermiteExpectation, MonteCarloExpectation
 using GPLikelihoods.TestInterface: test_interface
-using Test
-using Random
-using Functors
+
+using Aqua
 using Distributions
+using Functors
+using Random
 using StatsFuns
+using Test
 using Zygote
 
 @testset "GPLikelihoods.jl" begin
@@ -20,4 +22,9 @@ using Zygote
         include("likelihoods/negativebinomial.jl")
     end
     include("expectations.jl")
+    @testset "Code quality (Aqua.jl)" begin
+        Aqua.test_all(GPLikelihoods; ambiguities=false)
+        # Ref https://github.com/JuliaTesting/Aqua.jl/issues/77
+        Aqua.test_ambiguities(GPLikelihoods; recursive=false)
+    end
 end