diff --git a/Project.toml b/Project.toml
index 2e59f4b6..692a5cb2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -19,7 +19,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [compat]
-AbstractGPs = "0.3, 0.4, 0.5"
+AbstractGPs = "0.6"
 ChainRulesCore = "1.7"
 Distributions = "0.25"
 FastGaussQuadrature = "0.4"
diff --git a/src/ApproximateGPs.jl b/src/ApproximateGPs.jl
index 83cb797e..30997653 100644
--- a/src/ApproximateGPs.jl
+++ b/src/ApproximateGPs.jl
@@ -19,6 +19,7 @@ export SparseVariationalApproximation
 export DefaultQuadrature, Analytic, GaussHermite, MonteCarlo
 
 include("utils.jl")
+include("latent_gp.jl")
 include("sparse_variational.jl")
 include("expected_loglik.jl")
 include("elbo.jl")
diff --git a/src/latent_gp.jl b/src/latent_gp.jl
new file mode 100644
index 00000000..27a6a776
--- /dev/null
+++ b/src/latent_gp.jl
@@ -0,0 +1,50 @@
+"""
+    LatentGP(f<:AbstractGP, lik, Σy)
+
+ - `f` is a `AbstractGP`.
+ - `lik` is the likelihood function which maps samples from `f` to the
+   corresponding conditional likelihood distributions (i.e., `lik` must return a
+   `Distribution` compatible with the observations).
+ - `Σy` is the noise under which the latent GP is "observed"; this represents
+   the jitter used to avoid numeric instability and should generally be small.
+"""
+struct LatentGP{Tf<:AbstractGP,Tlik,TΣy}
+    f::Tf
+    lik::Tlik
+    Σy::TΣy
+end
+
+"""
+    LatentFiniteGP(fx<:FiniteGP, lik)
+
+ - `fx` is a `FiniteGP`.
+ - `lik` is the likelihood function which maps samples from `f` to the
+   corresponding conditional likelihood distributions (i.e., `lik` must return a
+   `Distribution` compatible with the observations).
+"""
+struct LatentFiniteGP{Tfx<:FiniteGP,Tlik}
+    fx::Tfx
+    lik::Tlik
+end
+
+(lgp::LatentGP)(x) = LatentFiniteGP(lgp.f(x, lgp.Σy), lgp.lik)
+
+Base.length(lgpx::LatentFiniteGP) = length(lgpx.fx)
+
+function Distributions.rand(rng::AbstractRNG, lfgp::LatentFiniteGP)
+    f = rand(rng, lfgp.fx)
+    y = rand(rng, lfgp.lik(f))
+    return (f=f, y=y)
+end
+
+"""
+    logpdf(lfgp::LatentFiniteGP, y::NamedTuple{(:f, :y)})
+
+```math
+    log p(y, f; x)
+```
+The joint log density of the Gaussian process output `f` and observation `y`.
+"""
+function Distributions.logpdf(lfgp::LatentFiniteGP, y::NamedTuple{(:f, :y)})
+    return logpdf(lfgp.fx, y.f) + logpdf(lfgp.lik(y.f), y.y)
+end
diff --git a/test/latent_gp.jl b/test/latent_gp.jl
new file mode 100644
index 00000000..506dcf92
--- /dev/null
+++ b/test/latent_gp.jl
@@ -0,0 +1,19 @@
+@testset "latent_gp" begin
+    gp = GP(SqExponentialKernel())
+    x = rand(10)
+    y = rand(10)
+
+    lgp = LatentGP(gp, x -> MvNormal(x, 0.1), 1e-5)
+    @test lgp isa LatentGP
+    @test lgp.f isa AbstractGPs.AbstractGP
+    @test lgp.Σy isa Real
+
+    lfgp = lgp(x)
+    @test lfgp isa AbstractGPs.LatentFiniteGP
+    @test lfgp.fx isa AbstractGPs.FiniteGP
+    @test length(lfgp) == length(x)
+
+    f = rand(10)
+    @test logpdf(lfgp, (f=f, y=y)) isa Real
+    @test rand(lfgp) isa NamedTuple{(:f, :y)}
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 744be20f..2a46041c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -21,6 +21,10 @@ const PKGDIR = dirname(dirname(pathof(ApproximateGPs)))
 include("test_utils.jl")
 
 @testset "ApproximateGPs" begin
+    include("latent_gp.jl")
+    println(" ")
+    @info "Ran latent_gp tests"
+
     include("expected_loglik.jl")
     println(" ")
     @info "Ran expected_loglik tests"