diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index f719c3b..022d807 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -48,7 +48,7 @@ jobs:
         with:
           file: lcov.info
 
-  # slack:
+  slack:
     name: Notify Slack Failure
     needs: test
     runs-on: ubuntu-latest
@@ -79,5 +79,6 @@ jobs:
             Pkg.instantiate()
             include("docs/make.jl")'
         env:
+          JULIA_PKG_SERVER: ""
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
           DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index b81699f..52ed95d 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -1,12 +1,87 @@
 # This file is machine-generated - editing it directly is not advised
 
+[[ANSIColoredPrinters]]
+git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
+uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9"
+version = "0.0.1"
+
+[[AbstractGPs]]
+deps = ["ChainRulesCore", "Distributions", "FillArrays", "KernelFunctions", "LinearAlgebra", "Random", "RecipesBase", "Reexport", "Statistics", "StatsBase", "Test"]
+git-tree-sha1 = "d1b311e6f0358fa73c82473daf23de6ba1148900"
+uuid = "99985d1d-32ba-4be9-9821-2ec096f28918"
+version = "0.3.10"
+
+[[ArgTools]]
+uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
+
+[[Artifacts]]
+uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+
 [[Base64]]
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 
+[[BlockDiagonals]]
+deps = ["ChainRulesCore", "FillArrays", "FiniteDifferences", "LinearAlgebra"]
+git-tree-sha1 = "0bcce621536e48a8a766909d5c0c1946bf956790"
+uuid = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
+version = "0.1.21"
+
+[[ChainRulesCore]]
+deps = ["Compat", "LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "bdc0937269321858ab2a4f288486cb258b9a0af7"
+uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+version = "1.3.0"
+
+[[Compat]]
+deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
+git-tree-sha1 = "79b9563ef3f2cc5fc6d3046a5ee1a57c9de52495"
+uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
+version = "3.33.0"
+
+[[CompilerSupportLibraries_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
+
+[[CompositionsBase]]
+git-tree-sha1 = "455419f7e328a1a2493cabc6428d79e951349769"
+uuid = "a33af91c-f02d-484b-be07-31d278c5ca2b"
+version = "0.1.1"
+
+[[DataAPI]]
+git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
+uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
+version = "1.7.0"
+
+[[DataStructures]]
+deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
+git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
+uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+version = "0.18.10"
+
 [[Dates]]
 deps = ["Printf"]
 uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
 
+[[DelimitedFiles]]
+deps = ["Mmap"]
+uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+
+[[Distances]]
+deps = ["LinearAlgebra", "Statistics", "StatsAPI"]
+git-tree-sha1 = "abe4ad222b26af3337262b8afb28fab8d215e9f8"
+uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+version = "0.10.3"
+
+[[Distributed]]
+deps = ["Random", "Serialization", "Sockets"]
+uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+
+[[Distributions]]
+deps = ["FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"]
+git-tree-sha1 = "3889f646423ce91dd1055a76317e9a1d3a23fff1"
+uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
+version = "0.25.11"
+
 [[DocStringExtensions]]
 deps = ["LibGit2"]
 git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
@@ -14,10 +89,32 @@ uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 version = "0.8.5"
 
 [[Documenter]]
-deps = ["Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
-git-tree-sha1 = "47f13b6305ab195edb73c86815962d84e31b0f48"
+deps = ["ANSIColoredPrinters", "Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
+git-tree-sha1 = "350dced36c11f794c6c4da5dc6493ec894e50c16"
 uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "0.27.3"
+version = "0.27.5"
+
+[[Downloads]]
+deps = ["ArgTools", "LibCURL", "NetworkOptions"]
+uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+
+[[FillArrays]]
+deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
+git-tree-sha1 = "8c8eac2af06ce35973c3eadb4ab3243076a408e7"
+uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
+version = "0.12.1"
+
+[[FiniteDifferences]]
+deps = ["ChainRulesCore", "LinearAlgebra", "Printf", "Random", "Richardson", "StaticArrays"]
+git-tree-sha1 = "9a586f04a21e6945f4cbee0d0fb6aebd7b86aa8f"
+uuid = "26cc04aa-876d-5657-8c51-4c34ba976000"
+version = "0.12.18"
+
+[[Functors]]
+deps = ["MacroTools"]
+git-tree-sha1 = "4cd9e70bf8fce05114598b663ad79dfe9ae432b3"
+uuid = "d9f16b24-f501-4c13-a1f2-28368ffc5196"
+version = "0.2.3"
 
 [[IOCapture]]
 deps = ["Logging", "Random"]
@@ -29,44 +126,141 @@ version = "0.2.2"
 deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 
+[[IrrationalConstants]]
+git-tree-sha1 = "f76424439413893a832026ca355fe273e93bce94"
+uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
+version = "0.1.0"
+
+[[JLLWrappers]]
+deps = ["Preferences"]
+git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
+uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
+version = "1.3.0"
+
 [[JSON]]
 deps = ["Dates", "Mmap", "Parsers", "Unicode"]
-git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4"
+git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
 uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
-version = "0.21.1"
+version = "0.21.2"
+
+[[KernelFunctions]]
+deps = ["ChainRulesCore", "Compat", "CompositionsBase", "Distances", "FillArrays", "Functors", "IrrationalConstants", "LinearAlgebra", "LogExpFunctions", "Random", "Requires", "SpecialFunctions", "StatsBase", "TensorCore", "Test", "ZygoteRules"]
+git-tree-sha1 = "6f46fb7fa868699dfbb6ae7973ba2825d3558ade"
+uuid = "ec8451be-7e33-11e9-00cf-bbf324bd1392"
+version = "0.10.16"
+
+[[LibCURL]]
+deps = ["LibCURL_jll", "MozillaCACerts_jll"]
+uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
+
+[[LibCURL_jll]]
+deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
+uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
 
 [[LibGit2]]
 deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
 uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
 
+[[LibSSH2_jll]]
+deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
+uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
+
+[[Libdl]]
+uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
+
+[[LinearAlgebra]]
+deps = ["Libdl"]
+uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+
 [[LinearMixingModels]]
-path = ".."
+deps = ["AbstractGPs", "BlockDiagonals", "ChainRulesCore", "Distributions", "FillArrays", "KernelFunctions", "LinearAlgebra", "Random", "Statistics"]
+git-tree-sha1 = "820f493cba1fe57a3848bbed039949c87bcab718"
+repo-rev = "master"
+repo-url = "https://github.com/invenia/LinearMixingModels.jl.git"
 uuid = "b8ce4b42-e81b-4a39-a84a-67f74a9a16dd"
 version = "0.1.0"
 
+[[LogExpFunctions]]
+deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
+git-tree-sha1 = "3d682c07e6dd250ed082f883dc88aee7996bf2cc"
+uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+version = "0.3.0"
+
 [[Logging]]
 uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
 
+[[MacroTools]]
+deps = ["Markdown", "Random"]
+git-tree-sha1 = "0fb723cd8c45858c22169b2e42269e53271a6df7"
+uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
+version = "0.5.7"
+
 [[Markdown]]
 deps = ["Base64"]
 uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 
+[[MbedTLS_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
+
+[[Missings]]
+deps = ["DataAPI"]
+git-tree-sha1 = "4ea90bd5d3985ae1f9a908bd4500ae88921c5ce7"
+uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
+version = "1.0.0"
+
 [[Mmap]]
 uuid = "a63ad114-7e13-5084-954f-fe012c677804"
 
+[[MozillaCACerts_jll]]
+uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
+
 [[NetworkOptions]]
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
 
+[[OpenSpecFun_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
+uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
+version = "0.5.5+0"
+
+[[OrderedCollections]]
+git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
+uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
+version = "1.4.1"
+
+[[PDMats]]
+deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
+git-tree-sha1 = "4dd403333bcf0909341cfe57ec115152f937d7d8"
+uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
+version = "0.11.1"
+
 [[Parsers]]
 deps = ["Dates"]
-git-tree-sha1 = "c8abc88faa3f7a3950832ac5d6e690881590d6dc"
+git-tree-sha1 = "477bf42b4d1496b454c10cce46645bb5b8a0cf2c"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "1.1.0"
+version = "2.0.2"
+
+[[Pkg]]
+deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
+uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+
+[[Preferences]]
+deps = ["TOML"]
+git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
+uuid = "21216c6a-2e73-6563-6e65-726566657250"
+version = "1.2.2"
 
 [[Printf]]
 deps = ["Unicode"]
 uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
+[[QuadGK]]
+deps = ["DataStructures", "LinearAlgebra"]
+git-tree-sha1 = "12fbe86da16df6679be7521dfb39fbc861e1dc7b"
+uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+version = "2.4.1"
+
 [[REPL]]
 deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
 uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
@@ -75,18 +269,139 @@ uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 deps = ["Serialization"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
+[[RecipesBase]]
+git-tree-sha1 = "b3fb709f3c97bfc6e948be68beeecb55a0b340ae"
+uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+version = "1.1.1"
+
+[[Reexport]]
+git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
+uuid = "189a3867-3050-52da-a836-e630ba90ab69"
+version = "1.1.0"
+
+[[Requires]]
+deps = ["UUIDs"]
+git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
+uuid = "ae029012-a4dd-5104-9daa-d747884805df"
+version = "1.1.3"
+
+[[Richardson]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "e03ca566bec93f8a3aeb059c8ef102f268a38949"
+uuid = "708f8203-808e-40c0-ba2d-98a6953ed40d"
+version = "1.4.0"
+
+[[Rmath]]
+deps = ["Random", "Rmath_jll"]
+git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
+uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
+version = "0.7.0"
+
+[[Rmath_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
+uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
+version = "0.3.0+0"
+
 [[SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 
 [[Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
 
+[[SharedArrays]]
+deps = ["Distributed", "Mmap", "Random", "Serialization"]
+uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+
 [[Sockets]]
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 
+[[SortingAlgorithms]]
+deps = ["DataStructures"]
+git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
+uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
+version = "1.0.1"
+
+[[SparseArrays]]
+deps = ["LinearAlgebra", "Random"]
+uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+
+[[SpecialFunctions]]
+deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"]
+git-tree-sha1 = "a322a9493e49c5f3a10b50df3aedaf1cdb3244b7"
+uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
+version = "1.6.1"
+
+[[StaticArrays]]
+deps = ["LinearAlgebra", "Random", "Statistics"]
+git-tree-sha1 = "3240808c6d463ac46f1c1cd7638375cd22abbccb"
+uuid = "90137ffa-7385-5640-81b9-e52037218182"
+version = "1.2.12"
+
+[[Statistics]]
+deps = ["LinearAlgebra", "SparseArrays"]
+uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[[StatsAPI]]
+git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
+uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
+version = "1.0.0"
+
+[[StatsBase]]
+deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
+git-tree-sha1 = "fed1ec1e65749c4d96fc20dd13bea72b55457e62"
+uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+version = "0.33.9"
+
+[[StatsFuns]]
+deps = ["IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
+git-tree-sha1 = "20d1bb720b9b27636280f751746ba4abb465f19d"
+uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
+version = "0.9.9"
+
+[[SuiteSparse]]
+deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
+uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
+
+[[TOML]]
+deps = ["Dates"]
+uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+
+[[Tar]]
+deps = ["ArgTools", "SHA"]
+uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
+
+[[TensorCore]]
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6"
+uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50"
+version = "0.1.1"
+
 [[Test]]
 deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
+[[UUIDs]]
+deps = ["Random", "SHA"]
+uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
+
 [[Unicode]]
 uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
+
+[[Zlib_jll]]
+deps = ["Libdl"]
+uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
+
+[[ZygoteRules]]
+deps = ["MacroTools"]
+git-tree-sha1 = "9e7a1e8ca60b742e508a315c17eef5211e7fbfd7"
+uuid = "700de1a5-db45-46bc-99cf-38207098b444"
+version = "0.2.1"
+
+[[nghttp2_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
+
+[[p7zip_jll]]
+deps = ["Artifacts", "Libdl"]
+uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
diff --git a/docs/Project.toml b/docs/Project.toml
index 9d246ee..0e1c594 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,3 +1,7 @@
 [deps]
+AbstractGPs = "99985d1d-32ba-4be9-9821-2ec096f28918"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+KernelFunctions = "ec8451be-7e33-11e9-00cf-bbf324bd1392"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMixingModels = "b8ce4b42-e81b-4a39-a84a-67f74a9a16dd"
diff --git a/docs/make.jl b/docs/make.jl
index 2ca78f0..845da44 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,7 +1,12 @@
 using LinearMixingModels
 using Documenter
 
-DocMeta.setdocmeta!(LinearMixingModels, :DocTestSetup, :(using LinearMixingModels); recursive=true)
+DocMeta.setdocmeta!(
+    LinearMixingModels,
+    :DocTestSetup,
+    :(using LinearMixingModels, KernelFunctions, FillArrays, LinearAlgebra, AbstractGPs);
+    recursive=true,
+)
 
 makedocs(;
     modules=[LinearMixingModels],
diff --git a/src/LinearMixingModels.jl b/src/LinearMixingModels.jl
index 0dc601a..c77dbfe 100644
--- a/src/LinearMixingModels.jl
+++ b/src/LinearMixingModels.jl
@@ -21,7 +21,6 @@ include("oilmm.jl")
 export ILMM
 export IndependentMOGP, independent_mogp
 export Orthogonal, OILMM
-
 export get_latent_gp
 
 end
diff --git a/src/ilmm.jl b/src/ilmm.jl
index 45ef6c8..ea49a86 100644
--- a/src/ilmm.jl
+++ b/src/ilmm.jl
@@ -25,68 +25,24 @@ Returns the underlying latent space AbstractGP belonging to `f`.
 
 
 ```jldoctest
-julia> f = ILMM(IndependentMOGP([GP(SEKernel())]), rand(2,2));
+julia> f = ILMM(independent_mogp([GP(SEKernel())]), rand(2,2));
 
 julia> latent_f = get_latent_gp(f);
 
-julia> latent_f == IndependentMOGP([GP(SEKernel())])
+julia> latent_f isa IndependentMOGP
 true
-```
-"""
-get_latent_gp(f::ILMM) = f.f
-
-"""
-    noise_var(Σ)
-
-Return the diagonal element of Σ.
-
-```jldoctest
-julia> Σ = Diagonal(Fill(2, 3));
 
-julia> noise_var(Σ) == 2
+julia> latent_f.fs == [GP(SEKernel())]
 true
 ```
 """
-noise_var(Σ::Diagonal{<:Real, <:Fill}) = FillArrays.getindex_value(Σ.diag)
-
-
-"""
-    reshape_y(y, N)
-
-Reshape `y` in to an adjoint Matrix of dimension (length(y)/N, N)`
-
-```jldoctest
-julia> y = rand(16);
+get_latent_gp(f::ILMM) = f.f
 
-julia> size(reshape_y(y,8)) == (2, 8)
-true
+noise_var(Σ::Diagonal{<:Real,<:Fill}) = FillArrays.getindex_value(Σ.diag)
 
-julia> size(reshape_y(y,2)) == (8, 2)
-true
-```
-"""
 reshape_y(y::AbstractVector{<:Real}, N::Int) = reshape(y, N, :)'
 
-"""
-    unpack(fx)
-
-Collect the relevant underlying fields of the Finite ILMM. This includes
-the latent space GP, the mixing matrix, the noise and the observations.
-
-```jldoctest
-julia> fs = independent_mogp([GP(Matern32Kernel())]);
-
-julia> H = rand(2,1);
-
-julia> x = MOInputIsotopicByOutputs(rand(2,2), 2);
-
-julia> ilmmx = ILMM(fs, H)(x, 0.1);
-
-julia> (fs, H, 0.1, x.x) == unpack(ilmmx)
-true
-```
-"""
-function unpack(fx::FiniteGP{<:ILMM, <:MOInputIsotopicByOutputs, <:Diagonal{<:Real, <:Fill}})
+function unpack(fx::FiniteGP{<:ILMM,<:MOInputIsotopicByOutputs,<:Diagonal{<:Real,<:Fill}})
     f = fx.f.f
     H = fx.f.H
     σ² = noise_var(fx.Σy)
diff --git a/src/independent_mogp.jl b/src/independent_mogp.jl
index 9ddfd2a..3b89138 100644
--- a/src/independent_mogp.jl
+++ b/src/independent_mogp.jl
@@ -17,7 +17,14 @@ end
 Returns an IndependentMOGP given a list of single output GPs `fs`.
 
 ```jldoctest
-julia> independent_mogp([GP(SEKernel())]) == IndependentMOGP([GP(SEKernel())])
+julia> ind_mogp1 = independent_mogp([GP(KernelFunctions.SEKernel())]);
+
+julia> ind_mogp2 = IndependentMOGP([GP(KernelFunctions.SEKernel())]);
+
+julia> typeof(ind_mogp1) == typeof(ind_mogp2)
+true
+
+julia> ind_mogp1.fs == ind_mogp2.fs
 true
 ```
 """
@@ -110,11 +117,10 @@ Posterior implementation for isotopic inputs, given diagonal Σy (OILMM).
 See AbstractGPs.jl API docs.
 """
 function AbstractGPs.posterior(
-    ft::IsotropicFiniteIndependentMOGP,
-    y::AbstractVector{<:Real},
+    ft::IsotropicFiniteIndependentMOGP, y::AbstractVector{<:Real}
 )
     finiteGPs = finite_gps(ft, ft.Σy[1])
-    ys = collect(eachcol(reshape(y, (length(ft.x.x),:))))
+    ys = collect(eachcol(reshape(y, (length(ft.x.x), :))))
     ind_posts = [AbstractGPs.posterior(fx, y_i) for (fx, y_i) in zip(finiteGPs, ys)]
-    return IndependentMOGP(ind_posts)
+    return independent_mogp(ind_posts)
 end
diff --git a/src/oilmm.jl b/src/oilmm.jl
index cff968e..64c0d52 100644
--- a/src/oilmm.jl
+++ b/src/oilmm.jl
@@ -124,6 +124,11 @@ function AbstractGPs.posterior(fx::FiniteGP{<:OILMM}, y::AbstractVector{<:Real})
     # Condition each latent process on the projected observations.
     y_rows = collect(eachrow(Ty))
     ΣT_rows = collect(eachrow(repeat(ΣT, 1, size(Ty, 2))))
-    fs_posterior = map((f, s, y) -> AbstractGPs.posterior(f(x, collect(s)), collect(y)), fs.fs, ΣT_rows, y_rows)
-    return ILMM(IndependentMOGP(fs_posterior), H)
+    fs_posterior = map(
+        (f, s, y) -> AbstractGPs.posterior(f(x, collect(s)), collect(y)),
+        fs.fs,
+        ΣT_rows,
+        y_rows,
+    )
+    return ILMM(independent_mogp(fs_posterior), H)
 end
diff --git a/test/Project.toml b/test/Project.toml
index 7656bb0..40653de 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -2,14 +2,11 @@
 AbstractGPs = "99985d1d-32ba-4be9-9821-2ec096f28918"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 KernelFunctions = "ec8451be-7e33-11e9-00cf-bbf324bd1392"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMixingModels = "b8ce4b42-e81b-4a39-a84a-67f74a9a16dd"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-Stheno = "8188c328-b5d6-583d-959b-9690869a5511"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
-
-[compat]
-Stheno = "0.7"
diff --git a/test/ilmm.jl b/test/ilmm.jl
index 1c07f19..b610fb6 100644
--- a/test/ilmm.jl
+++ b/test/ilmm.jl
@@ -2,17 +2,17 @@
     rng = Random.seed!(04161999)
     atol = 1e-2
     @testset "Full Rank, Dense H" begin
-        H = rand(3,3)
+        H = rand(3, 3)
         k1, k2, k3 = SEKernel(), Matern32Kernel(), Matern32Kernel()
 
-        ilmm = ILMM(IndependentMOGP([GP(k1), GP(k2), GP(k3)]), H)
+        ilmm = ILMM(independent_mogp([GP(k1), GP(k2), GP(k3)]), H)
         n_ilmm = GP(LinearMixingModelKernel([k1, k2, k3], H'))
 
-        x = range(0,10;length=5);
+        x = range(0, 10; length=5)
         ys = rand(rng, GP(SEKernel())(x, 1e-6), 3)
-        y1 = ys[:,1]
-        y2 = ys[:,2]
-        y3 = ys[:,3]
+        y1 = ys[:, 1]
+        y2 = ys[:, 2]
+        y3 = ys[:, 3]
         indices = randcycle(rng, 5)
         x_train = zeros(3)
         y_1_train = zeros(3)
@@ -38,25 +38,25 @@
         x_train = MOInputIsotopicByOutputs(x_train, 3)
         x_test = MOInputIsotopicByOutputs(x_test, 3)
         y_train = vcat(y_1_train, y_2_train, y_3_train)
-        y_test = vcat(y_1_test, y_2_test, y_3_test);
+        y_test = vcat(y_1_test, y_2_test, y_3_test)
 
         ilmmx = ilmm(x_train, 1e-6)
         n_ilmmx = n_ilmm(x_train, 1e-6)
 
-        @test isapprox(mean(ilmmx), mean(n_ilmmx), atol=atol)
-        @test isapprox(var(ilmmx), var(n_ilmmx), atol=atol)
-        @test isapprox(logpdf(ilmmx, y_train), logpdf(n_ilmmx, y_train), atol=atol)
+        @test isapprox(mean(ilmmx), mean(n_ilmmx); atol=atol)
+        @test isapprox(var(ilmmx), var(n_ilmmx); atol=atol)
+        @test isapprox(logpdf(ilmmx, y_train), logpdf(n_ilmmx, y_train); atol=atol)
         @test _is_approx(marginals(ilmmx), marginals(n_ilmmx))
 
-        p_ilmmx = posterior(ilmmx, y_train);
-        p_n_ilmmx = posterior(n_ilmmx, y_train);
+        p_ilmmx = posterior(ilmmx, y_train)
+        p_n_ilmmx = posterior(n_ilmmx, y_train)
 
-        pi = p_ilmmx(x_test, 1e-6);
-        pni = p_n_ilmmx(x_test, 1e-6);
+        pi = p_ilmmx(x_test, 1e-6)
+        pni = p_n_ilmmx(x_test, 1e-6)
 
-        @test isapprox(mean(pi), mean(pni), atol=atol)
-        @test isapprox(var(pi), var(pni), atol=atol)
-        @test isapprox(logpdf(pi, y_test), logpdf(pni, y_test), atol=atol)
+        @test isapprox(mean(pi), mean(pni); atol=atol)
+        @test isapprox(var(pi), var(pni); atol=atol)
+        @test isapprox(logpdf(pi, y_test), logpdf(pni, y_test); atol=atol)
         @test _is_approx(marginals(pi), marginals(pni))
 
         @testset "primary_public_interface" begin
@@ -70,14 +70,14 @@
 
         k1, k2 = SEKernel(), Matern32Kernel()
 
-        ilmm = ILMM(IndependentMOGP([GP(k1), GP(k2)]), H)
+        ilmm = ILMM(independent_mogp([GP(k1), GP(k2)]), H)
         n_ilmm = GP(LinearMixingModelKernel([k1, k2], H'))
 
-        x = range(0,10;length=5);
+        x = range(0, 10; length=5)
         ys = rand(rng, GP(SEKernel())(x, 1e-6), 3)
-        y1 = ys[:,1]
-        y2 = ys[:,2]
-        y3 = ys[:,3]
+        y1 = ys[:, 1]
+        y2 = ys[:, 2]
+        y3 = ys[:, 3]
         indices = randcycle(rng, 5)
         x_train = zeros(3)
         y_1_train = zeros(3)
@@ -103,25 +103,25 @@
         x_train = MOInputIsotopicByOutputs(x_train, 3)
         x_test = MOInputIsotopicByOutputs(x_test, 3)
         y_train = vcat(y_1_train, y_2_train, y_3_train)
-        y_test = vcat(y_1_test, y_2_test, y_3_test);
+        y_test = vcat(y_1_test, y_2_test, y_3_test)
 
         ilmmx = ilmm(x_train, 1e-6)
         n_ilmmx = n_ilmm(x_train, 1e-6)
 
-        @test isapprox(mean(ilmmx), mean(n_ilmmx), atol=atol)
-        @test isapprox(var(ilmmx), var(n_ilmmx), atol=atol)
-        @test isapprox(logpdf(ilmmx, y_train), logpdf(n_ilmmx, y_train), atol=atol)
+        @test isapprox(mean(ilmmx), mean(n_ilmmx); atol=atol)
+        @test isapprox(var(ilmmx), var(n_ilmmx); atol=atol)
+        @test isapprox(logpdf(ilmmx, y_train), logpdf(n_ilmmx, y_train); atol=atol)
         @test _is_approx(marginals(ilmmx), marginals(n_ilmmx))
 
-        p_ilmmx = posterior(ilmmx, y_train);
-        p_n_ilmmx = posterior(n_ilmmx, y_train);
+        p_ilmmx = posterior(ilmmx, y_train)
+        p_n_ilmmx = posterior(n_ilmmx, y_train)
 
-        pi = p_ilmmx(x_test, 1e-6);
-        pni = p_n_ilmmx(x_test, 1e-6);
+        pi = p_ilmmx(x_test, 1e-6)
+        pni = p_n_ilmmx(x_test, 1e-6)
 
-        @test isapprox(mean(pi), mean(pni), atol=atol)
-        @test isapprox(var(pi), var(pni), atol=atol)
-        @test isapprox(logpdf(pi, y_test), logpdf(pni, y_test), atol=atol)
+        @test isapprox(mean(pi), mean(pni); atol=atol)
+        @test isapprox(var(pi), var(pni); atol=atol)
+        @test isapprox(logpdf(pi, y_test), logpdf(pni, y_test); atol=atol)
         @test _is_approx(marginals(pi), marginals(pni))
 
         @testset "primary_public_interface" begin
@@ -135,14 +135,14 @@
 
         k = SEKernel()
 
-        ilmm = ILMM(IndependentMOGP([GP(k)]), H)
+        ilmm = ILMM(independent_mogp([GP(k)]), H)
         n_ilmm = GP(LinearMixingModelKernel([k], H'))
 
-        x = range(0,10;length=5);
+        x = range(0, 10; length=5)
         ys = rand(rng, GP(SEKernel())(x, 1e-6), 3)
-        y1 = ys[:,1]
-        y2 = ys[:,2]
-        y3 = ys[:,3]
+        y1 = ys[:, 1]
+        y2 = ys[:, 2]
+        y3 = ys[:, 3]
         indices = randcycle(rng, 5)
         x_train = zeros(3)
         y_1_train = zeros(3)
@@ -168,27 +168,45 @@
         x_train = MOInputIsotopicByOutputs(x_train, 3)
         x_test = MOInputIsotopicByOutputs(x_test, 3)
         y_train = vcat(y_1_train, y_2_train, y_3_train)
-        y_test = vcat(y_1_test, y_2_test, y_3_test);
+        y_test = vcat(y_1_test, y_2_test, y_3_test)
 
         ilmmx = ilmm(x_train, 1e-6)
         n_ilmmx = n_ilmm(x_train, 1e-6)
 
-        @test isapprox(mean(ilmmx), mean(n_ilmmx), atol=atol)
-        @test isapprox(var(ilmmx), var(n_ilmmx), atol=atol)
-        @test isapprox(logpdf(ilmmx, y_train), logpdf(n_ilmmx, y_train), atol=atol)
+        @test isapprox(mean(ilmmx), mean(n_ilmmx); atol=atol)
+        @test isapprox(var(ilmmx), var(n_ilmmx); atol=atol)
+        @test isapprox(logpdf(ilmmx, y_train), logpdf(n_ilmmx, y_train); atol=atol)
         @test marginals(ilmmx) == marginals(n_ilmmx)
 
-        p_ilmmx = posterior(ilmmx, y_train);
-        p_n_ilmmx = posterior(n_ilmmx, y_train);
+        p_ilmmx = posterior(ilmmx, y_train)
+        p_n_ilmmx = posterior(n_ilmmx, y_train)
 
-        pi = p_ilmmx(x_test, 1e-6);
-        pni = p_n_ilmmx(x_test, 1e-6);
+        pi = p_ilmmx(x_test, 1e-6)
+        pni = p_n_ilmmx(x_test, 1e-6)
 
-        @test isapprox(mean(pi), mean(pni), atol=atol)
-        @test isapprox(var(pi), var(pni), atol=atol)
-        @test isapprox(logpdf(pi, y_test), logpdf(pni, y_test), atol=atol)
+        @test isapprox(mean(pi), mean(pni); atol=atol)
+        @test isapprox(var(pi), var(pni); atol=atol)
+        @test isapprox(logpdf(pi, y_test), logpdf(pni, y_test); atol=atol)
         @test _is_approx(marginals(pi), marginals(pni))
 
+        @testset "util" begin
+            Σ = Diagonal(Fill(2, 3))
+            @test noise_var(Σ) == 2
+
+            y = rand(16)
+            @test size(reshape_y(y, 8)) == (2, 8)
+            @test size(reshape_y(y, 2)) == (8, 2)
+
+            fs = independent_mogp([GP(Matern32Kernel())])
+            H = rand(2, 1)
+            x = MOInputIsotopicByOutputs(ColVecs(rand(2, 2)), 2)
+            ilmm = ILMM(fs, H)
+            ilmmx = ilmm(x, 0.1)
+            @test (fs, H, 0.1, x.x) == unpack(ilmmx)
+
+            @test get_latent_gp(ilmm) == fs
+        end
+
         @testset "primary_public_interface" begin
             test_finitegp_primary_public_interface(rng, ilmmx)
             test_finitegp_primary_public_interface(rng, pi)
diff --git a/test/oilmm.jl b/test/oilmm.jl
index 1d242dd..6571822 100644
--- a/test/oilmm.jl
+++ b/test/oilmm.jl
@@ -3,17 +3,17 @@
     atol = 1e-2
     @testset "Full Rank, Dense H" begin
         U, S, _ = svd(rand(rng, 3, 3))
-        H = Orthogonal(U, Diagonal(S));
-        fs = IndependentMOGP([GP(SEKernel()), GP(Matern32Kernel()), GP(Matern32Kernel())])
+        H = Orthogonal(U, Diagonal(S))
+        fs = independent_mogp([GP(SEKernel()), GP(Matern32Kernel()), GP(Matern32Kernel())])
 
         ilmm = ILMM(fs, collect(H))
         oilmm = ILMM(fs, H)
 
-        x = range(0,10;length=5);
+        x = range(0, 10; length=5)
         ys = rand(rng, GP(SEKernel())(x, 1e-6), 3)
-        y1 = ys[:,1]
-        y2 = ys[:,2]
-        y3 = ys[:,3]
+        y1 = ys[:, 1]
+        y2 = ys[:, 2]
+        y3 = ys[:, 3]
         indices = randcycle(rng, 5)
         x_train = zeros(3)
         y_1_train = zeros(3)
@@ -39,25 +39,25 @@
         x_train = MOInputIsotopicByOutputs(x_train, 3)
         x_test = MOInputIsotopicByOutputs(x_test, 3)
         y_train = vcat(y_1_train, y_2_train, y_3_train)
-        y_test = vcat(y_1_test, y_2_test, y_3_test);
+        y_test = vcat(y_1_test, y_2_test, y_3_test)
 
         ilmmx = ilmm(x_train, 0.1)
         oilmmx = oilmm(x_train, 0.1)
 
-        @test isapprox(mean(ilmmx), mean(oilmmx), atol=atol)
-        @test isapprox(var(ilmmx), var(oilmmx), atol=atol)
-        @test isapprox(logpdf(ilmmx, y_train), logpdf(oilmmx, y_train), atol=atol)
+        @test isapprox(mean(ilmmx), mean(oilmmx); atol=atol)
+        @test isapprox(var(ilmmx), var(oilmmx); atol=atol)
+        @test isapprox(logpdf(ilmmx, y_train), logpdf(oilmmx, y_train); atol=atol)
         @test _is_approx(marginals(ilmmx), marginals(oilmmx))
 
-        p_ilmmx = posterior(ilmmx, y_train);
-        p_oilmmx = posterior(oilmmx, y_train);
+        p_ilmmx = posterior(ilmmx, y_train)
+        p_oilmmx = posterior(oilmmx, y_train)
 
-        pi = p_ilmmx(x_test, 0.1);
-        po = p_oilmmx(x_test, 0.1);
+        pi = p_ilmmx(x_test, 0.1)
+        po = p_oilmmx(x_test, 0.1)
 
-        @test isapprox(mean(pi), mean(po), atol=atol)
-        @test isapprox(var(pi), var(po), atol=atol)
-        @test isapprox(logpdf(pi, y_test), logpdf(po, y_test), atol=atol)
+        @test isapprox(mean(pi), mean(po); atol=atol)
+        @test isapprox(var(pi), var(po); atol=atol)
+        @test isapprox(logpdf(pi, y_test), logpdf(po, y_test); atol=atol)
         @test _is_approx(marginals(pi), marginals(po))
 
         @testset "primary_public_interface" begin
@@ -68,16 +68,16 @@
 
     @testset "M Latent Processes" begin
         U, S, _ = svd(rand(rng, 3, 2))
-        H = Orthogonal(U, Diagonal(S));
-        fs = IndependentMOGP([GP(SEKernel()), GP(Matern32Kernel())])
+        H = Orthogonal(U, Diagonal(S))
+        fs = independent_mogp([GP(SEKernel()), GP(Matern32Kernel())])
         ilmm = ILMM(fs, collect(H))
         oilmm = ILMM(fs, H)
 
-        x = range(0,10;length=5);
+        x = range(0, 10; length=5)
         ys = rand(rng, GP(SEKernel())(x, 1e-6), 3)
-        y1 = ys[:,1]
-        y2 = ys[:,2]
-        y3 = ys[:,3]
+        y1 = ys[:, 1]
+        y2 = ys[:, 2]
+        y3 = ys[:, 3]
         indices = randcycle(rng, 5)
         x_train = zeros(3)
         y_1_train = zeros(3)
@@ -103,25 +103,25 @@
         x_train = MOInputIsotopicByOutputs(x_train, 3)
         x_test = MOInputIsotopicByOutputs(x_test, 3)
         y_train = vcat(y_1_train, y_2_train, y_3_train)
-        y_test = vcat(y_1_test, y_2_test, y_3_test);
+        y_test = vcat(y_1_test, y_2_test, y_3_test)
 
         ilmmx = ilmm(x_train, 0.1)
         oilmmx = oilmm(x_train, 0.1)
 
-        @test isapprox(mean(ilmmx), mean(oilmmx), atol=atol)
-        @test isapprox(var(ilmmx), var(oilmmx), atol=atol)
-        @test isapprox(logpdf(ilmmx, y_train), logpdf(oilmmx, y_train), atol=atol)
+        @test isapprox(mean(ilmmx), mean(oilmmx); atol=atol)
+        @test isapprox(var(ilmmx), var(oilmmx); atol=atol)
+        @test isapprox(logpdf(ilmmx, y_train), logpdf(oilmmx, y_train); atol=atol)
         @test _is_approx(marginals(ilmmx), marginals(oilmmx))
 
-        p_ilmmx = posterior(ilmmx, y_train);
-        p_oilmmx = posterior(oilmmx, y_train);
+        p_ilmmx = posterior(ilmmx, y_train)
+        p_oilmmx = posterior(oilmmx, y_train)
 
-        pi = p_ilmmx(x_test, 0.1);
-        po = p_oilmmx(x_test, 0.1);
+        pi = p_ilmmx(x_test, 0.1)
+        po = p_oilmmx(x_test, 0.1)
 
-        @test isapprox(mean(pi), mean(po), atol=atol)
-        @test isapprox(var(pi), var(po), atol=atol)
-        @test isapprox(logpdf(pi, y_test), logpdf(po, y_test), atol=atol)
+        @test isapprox(mean(pi), mean(po); atol=atol)
+        @test isapprox(var(pi), var(po); atol=atol)
+        @test isapprox(logpdf(pi, y_test), logpdf(po, y_test); atol=atol)
         @test _is_approx(marginals(pi), marginals(po))
 
         @testset "primary_public_interface" begin
@@ -132,17 +132,17 @@
 
     @testset "1 Latent Processes" begin
         U, S, _ = svd(rand(rng, 3, 1))
-        H = Orthogonal(U, Diagonal(S));
-        fs = IndependentMOGP([GP(SEKernel())])
+        H = Orthogonal(U, Diagonal(S))
+        fs = independent_mogp([GP(SEKernel())])
 
         ilmm = ILMM(fs, collect(H))
         oilmm = ILMM(fs, H)
 
-        x = range(0,10;length=5);
+        x = range(0, 10; length=5)
         ys = rand(rng, GP(SEKernel())(x, 1e-6), 3)
-        y1 = ys[:,1]
-        y2 = ys[:,2]
-        y3 = ys[:,3]
+        y1 = ys[:, 1]
+        y2 = ys[:, 2]
+        y3 = ys[:, 3]
         indices = randcycle(rng, 5)
         x_train = zeros(3)
         y_1_train = zeros(3)
@@ -168,25 +168,25 @@
         x_train = MOInputIsotopicByOutputs(x_train, 3)
         x_test = MOInputIsotopicByOutputs(x_test, 3)
         y_train = vcat(y_1_train, y_2_train, y_3_train)
-        y_test = vcat(y_1_test, y_2_test, y_3_test);
+        y_test = vcat(y_1_test, y_2_test, y_3_test)
 
         ilmmx = ilmm(x_train, 0.1)
         oilmmx = oilmm(x_train, 0.1)
 
-        @test isapprox(mean(ilmmx), mean(oilmmx), atol=atol)
-        @test isapprox(var(ilmmx), var(oilmmx), atol=atol)
-        @test isapprox(logpdf(ilmmx, y_train), logpdf(oilmmx, y_train), atol=atol)
+        @test isapprox(mean(ilmmx), mean(oilmmx); atol=atol)
+        @test isapprox(var(ilmmx), var(oilmmx); atol=atol)
+        @test isapprox(logpdf(ilmmx, y_train), logpdf(oilmmx, y_train); atol=atol)
         @test _is_approx(marginals(ilmmx), marginals(oilmmx))
 
-        p_ilmmx = posterior(ilmmx, y_train);
-        p_oilmmx = posterior(oilmmx, y_train);
+        p_ilmmx = posterior(ilmmx, y_train)
+        p_oilmmx = posterior(oilmmx, y_train)
 
-        pi = p_ilmmx(x_test, 0.1);
-        po = p_oilmmx(x_test, 0.1);
+        pi = p_ilmmx(x_test, 0.1)
+        po = p_oilmmx(x_test, 0.1)
 
-        @test isapprox(mean(pi), mean(po), atol=atol)
-        @test isapprox(var(pi), var(po), atol=atol)
-        @test isapprox(logpdf(pi, y_test), logpdf(po, y_test), atol=atol)
+        @test isapprox(mean(pi), mean(po); atol=atol)
+        @test isapprox(var(pi), var(po); atol=atol)
+        @test isapprox(logpdf(pi, y_test), logpdf(po, y_test); atol=atol)
         @test _is_approx(marginals(pi), marginals(po))
 
         @testset "primary_public_interface" begin
diff --git a/test/runtests.jl b/test/runtests.jl
index 8d60746..ca03ba1 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,6 +1,7 @@
 using AbstractGPs
 using Distributions
 using Documenter
+using FillArrays
 using KernelFunctions
 using LinearAlgebra
 using LinearMixingModels