diff --git a/src/abstract_gp.jl b/src/abstract_gp.jl
index cba746f1..648f09a9 100644
--- a/src/abstract_gp.jl
+++ b/src/abstract_gp.jl
@@ -200,6 +200,8 @@ function logpdf(f::FiniteGP, Y::AbstractMatrix{<:Real})
     return -((size(Y, 1) * T(log(2π)) + logdet(C)) .+ diag_Xt_invA_X(C, Y .- μ)) ./ 2
 end
 
+Distributions.loglikelihood(f::FiniteGP, y::AbstractVector{<:Real}) = logpdf(f, y)
+Distributions.loglikelihood(f::FiniteGP, y::AbstractMatrix{<:Real}) = sum(logpdf(f, y))
 """
    elbo(f::FiniteGP, y::AbstractVector{<:Real}, u::FiniteGP)
 
diff --git a/test/Project.toml b/test/Project.toml
index 0494cf24..40027763 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -10,6 +10,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
@@ -20,5 +21,6 @@ Documenter = "0.25"
 FiniteDifferences = "0.10"
 Flux = "0.9, 0.10, 0.11"
 TimerOutputs = "0.5"
+Turing = "0.14.10"
 Zygote = "0.5"
 julia = "1.3"
diff --git a/test/integration/turing.jl b/test/integration/turing.jl
new file mode 100644
index 00000000..be8042b0
--- /dev/null
+++ b/test/integration/turing.jl
@@ -0,0 +1,24 @@
+# Ensure that using Stheno inside of Turing doesn't error.
+@testset "turing" begin
+
+    # make some fake data
+    l = 0.4
+    σ² = 1.3
+    σ²_n = 0.05 # noise
+    x_ = collect(range(-4.0, 4.0; length=10))
+    y_ = rand(GP(σ² * stretch(Matern52(), 1 / l), GPC())(x_, σ²_n))
+
+    # prior model
+    @model gp0(y,x) = begin
+        σ² ~ LogNormal(0, 1)
+        l ~ LogNormal(0, 1)
+        σ²_n ~ LogNormal(0, 1)  
+        k =  σ² * stretch(Matern52(), 1 ./ l)
+        f = GP(k, GPC())
+        y ~ f(x, σ²_n + 1e-3)
+    end
+
+    # sample from posterior
+    m = gp0(y_, x_)
+    chain = sample(m, HMC(0.01, 100), 2)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 7069371d..35398fed 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,6 +9,7 @@ using Statistics
 using Stheno
 using Test
 using TimerOutputs
+using Turing
 using Zygote
 
 using Stheno: ew, pw, mean_vector, cov, cov_diag
@@ -67,7 +68,7 @@ include("test_util.jl")
         include(joinpath("flux", "neural_kernel_network.jl"))
     end
 
-    println("doctests")
+    println("doctests:")
     @timedtestset "doctests" begin
         DocMeta.setdocmeta!(
             Stheno,
@@ -77,6 +78,11 @@ include("test_util.jl")
         )
         doctest(Stheno)
     end
+
+    println("turing:")
+    @timedtestset "turing" begin
+        include(joinpath("integration", "turing.jl"))
+    end
 end
 
 display(to)