diff --git a/.github/workflows/documenter.yml b/.github/workflows/documenter.yml
index 5e0e5321..3c24fb83 100644
--- a/.github/workflows/documenter.yml
+++ b/.github/workflows/documenter.yml
@@ -18,18 +18,17 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v2
+      - uses: quarto-dev/quarto-actions/setup@v2
+        with:
+          version: 1.3.353
       - uses: julia-actions/setup-julia@v1
       - uses: julia-actions/cache@v1
-      - name: Cache CondaPkg
-        id: cache-condaPkg
-        uses: actions/cache@v3
-        env:
-          cache-name: cache-condapkg
+      - name: Cache CmdStan
+        id: cache-cmdstan
+        uses: actions/cache@v2
         with:
-          path: docs/.CondaPkg
-          key: ${{ runner.os }}-${{ env.cache-name }}-${{ hashFiles('docs/CondaPkg.toml') }}
-          restore-keys: |
-            ${{ runner.os }}-${{ env.cache-name }}-
+          path: ${{ env.CMDSTAN_PATH }}
+          key: cmdstan-${{ env.CMDSTAN_VERSION }}-${{ runner.os }}
       - name: Download and build CmdStan
         if: steps.cache-cmdstan.outputs.cache-hit != 'true'
         run: |
diff --git a/docs/.gitignore b/docs/.gitignore
index ffba55c0..bbb3d394 100644
--- a/docs/.gitignore
+++ b/docs/.gitignore
@@ -1,9 +1,19 @@
 Manifest.toml
-.CondaPkg
-src/quickstart.md
-src/creating_custom_plots.md
+
+/.CondaPkg
+
+/src/.quarto/
+/src/_freeze/
+
+/src/quickstart.md
+/src/quickstart_files
+
+/src/creating_custom_plots.md
+/src/creating_custom_plots_files
+
 src/*.log
 src/assets/logo*.png
 src/assets/favicon.ico
 src/*.log
+
 build
diff --git a/docs/CondaPkg.toml b/docs/CondaPkg.toml
new file mode 100644
index 00000000..e5ff4db9
--- /dev/null
+++ b/docs/CondaPkg.toml
@@ -0,0 +1,2 @@
+[deps]
+jupyter = "1"
diff --git a/docs/Project.toml b/docs/Project.toml
index 5dda926a..b2ffe640 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -4,20 +4,20 @@ ArviZ = "131c737c-5715-5e2e-ad31-c244f01c1dc7"
 ArviZExampleData = "2f96bb34-afd9-46ae-bcd0-9b2d4372fe3c"
 ArviZPythonPlots = "4a6e88f0-2c8e-11ee-0601-e94153f0eada"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 EvoTrees = "f6006082-12f8-11e9-0c9c-0d5d367ab1e5"
+IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
 InferenceObjects = "b5cf5a8d-e756-4ee3-b014-01d49d192c00"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MCMCDiagnosticTools = "be115224-59cd-429b-ad48-344e309966f0"
 MLJBase = "a7f614a8-145f-11e9-1d2a-a57a1082229d"
 MLJIteration = "614be32b-d00c-4edb-bd02-1eb411ab5e55"
 PSIS = "ce719bf2-d5d0-4fb9-925d-10a81b42ad04"
-PlutoStaticHTML = "359b1769-a58e-495b-9770-312e911026ad"
-PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
 PosteriorStats = "7f36be82-ad55-44ba-a5c0-b8b5480d7aa5"
 StanSample = "c1514b29-d3a0-5178-b312-660c88baa699"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
@@ -30,16 +30,16 @@ ArviZ = ">=0.10"
 ArviZExampleData = "0.1.5"
 ArviZPythonPlots = "0.1"
 CairoMakie = "0.8.9, 0.9, 0.10"
+CondaPkg = "0.2"
 DataFrames = "1"
 DimensionalData = "0.23, 0.24"
 Distributions = "0.25"
 Documenter = "0.27"
 EvoTrees = "0.14.8, 0.15"
+IJulia = "1"
 MCMCDiagnosticTools = "0.3"
 MLJBase = "0.21.6"
 MLJIteration = "0.5.1"
-PlutoStaticHTML = "4.0.5, 5, 6"
-PlutoUI = "0.7"
 PosteriorStats = "0.1.1"
 StanSample = "7.0.1"
 StatsBase = "0.32, 0.33, 0.34"
diff --git a/docs/make.jl b/docs/make.jl
index 445d662b..18f50cbd 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,17 +1,47 @@
-using Documenter, Downloads, ArviZ
-using PlutoStaticHTML: PlutoStaticHTML
+using Pkg, CondaPkg, Documenter, Downloads, ArviZ
 
 const DOCS_SRC_PATH = joinpath(@__DIR__, "src")
 
-# generate markdown from Pluto notebooks
-output_format = PlutoStaticHTML.documenter_output
-build_opts = PlutoStaticHTML.BuildOptions(
-    DOCS_SRC_PATH;
-    previous_dir=DOCS_SRC_PATH,
-    output_format=output_format,
-    add_documenter_css=false,
-)
-PlutoStaticHTML.build_notebooks(build_opts)
+# generate markdown from Quarto files
+if Sys.which("quarto") !== nothing
+    CondaPkg.withenv() do
+        @info "Rendering Quarto files"
+        Pkg.build("IJulia")
+        run(`quarto render $(DOCS_SRC_PATH)`)
+    end
+else
+    @warn "Quarto not found, skipping rendering Quarto files"
+end
+
+function wrap_html_divs_in_raw_block(out_io, in_io)
+    level = 0
+    for line in eachline(in_io)
+        if contains(line, "<div")
+            level == 0 && write(out_io, "```@raw html\n")
+        end
+        write(out_io, line, '\n')
+        level += count("<div", line) - count("</div>", line)
+        if contains(line, "</div")
+            level == 0 && write(out_io, "```\n")
+        end
+    end
+    return nothing
+end
+
+# workaround for Julia markdown not supporting HTML
+for fn in readdir(DOCS_SRC_PATH; join=true)
+    root, ext = splitext(fn)
+    ext == ".qmd" || continue
+    md_file = root * ".md"
+    if isfile(md_file)
+        md_file_new, out_io = mktemp(; cleanup=false)
+        open(md_file, "r") do in_io
+            wrap_html_divs_in_raw_block(out_io, in_io)
+        end
+        close(out_io)
+        mv(md_file_new, md_file; force=true)
+    end
+end
 
 const ASSETS_DIR = joinpath(@__DIR__, "src", "assets")
 const ARVIZ_ASSETS_URL = "https://raw.githubusercontent.com/arviz-devs/arviz-project/main/arviz_logos"
@@ -98,8 +128,7 @@ makedocs(;
     doctestfilters,
     linkcheck=true,
     analytics="G-W1G68W77YV",
-    # allow linkcheck to fail so we can use pretty links to PlutoStaticHTML notebooks
-    strict=Documenter.except(:footnote, :linkcheck, :missing_docs),
+    strict=Documenter.except(:footnote, :missing_docs),
 )
 
 deploydocs(; repo="github.com/arviz-devs/ArviZ.jl.git", devbranch="main", push_preview=true)
diff --git a/docs/src/.gitignore b/docs/src/.gitignore
new file mode 100644
index 00000000..075b2542
--- /dev/null
+++ b/docs/src/.gitignore
@@ -0,0 +1 @@
+/.quarto/
diff --git a/docs/src/_quarto.yml b/docs/src/_quarto.yml
new file mode 100644
index 00000000..da7fe115
--- /dev/null
+++ b/docs/src/_quarto.yml
@@ -0,0 +1,16 @@
+project:
+  title: "ArviZ.jl"
+  output-dir: .
+  render:
+    - "*.qmd"
+
+execute:
+  freeze: auto
+  eval: true
+  echo: true
+  output: true
+
+format:
+  commonmark:
+    variant: -raw_html+tex_math_dollars
+    wrap: preserve
diff --git a/docs/src/creating_custom_plots.jl b/docs/src/creating_custom_plots.qmd
similarity index 52%
rename from docs/src/creating_custom_plots.jl
rename to docs/src/creating_custom_plots.qmd
index 68fc7974..a7724b69 100644
--- a/docs/src/creating_custom_plots.jl
+++ b/docs/src/creating_custom_plots.qmd
@@ -1,85 +1,57 @@
-### A Pluto.jl notebook ###
-# v0.19.16
-
-using Markdown
-using InteractiveUtils
-
-# ╔═╡ e5e15f39-027c-42aa-a765-588ec1cd0b63
-using Pkg, InteractiveUtils
-
-# ╔═╡ b1784f52-0cf2-11ed-32ba-69a8f63b48a9
-# Remove this cell to use release versions of dependencies
-# hideall
-let
-    docs_dir = dirname(@__DIR__)
-    pkg_dir = dirname(docs_dir)
-
-    Pkg.activate(docs_dir)
-    Pkg.develop(; path=pkg_dir)
-    Pkg.instantiate()
-end;
-
-# ╔═╡ 34ee7907-49e7-4e32-9bd5-1c5d490a13d8
-begin
-    using ArviZ, ArviZExampleData, DimensionalData, DataFrames, Statistics
-    using AlgebraOfGraphics, CairoMakie
-    using AlgebraOfGraphics: density
-    set_aog_theme!()
-end;
-
-# ╔═╡ 41b1ddca-20dc-4821-9cec-d7192ef34795
-using PlutoUI
-
-# ╔═╡ 70f66296-349c-48a5-9b83-cc5c6cdbd514
-md"""
-# Creating custom plots
-"""
-
-# ╔═╡ 47c7a1fd-0862-49d5-879d-d59d321ce014
-md"""
-While ArviZ includes many plotting functions for visualizing the data stored in `InferenceData` objects, you will often need to construct custom plots, or you may want to tweak some of our plots in your favorite plotting package.
+---
+title: "Creating custom plots"
+---
+
+While ArviZ includes many plotting functions for visualizing the data stored in [`InferenceData`](@ref) objects, you will often need to construct custom plots, or you may want to tweak some of our plots in your favorite plotting package.
 
 In this tutorial, we will show you a few useful techniques you can use to construct these plots using Julia's plotting packages.
-For demonstration purposes, we'll use Makie.jl and AlgebraOfGraphics.jl, which can consume `Dataset` objects since they implement the Tables interface.
+For demonstration purposes, we'll use Makie.jl and AlgebraOfGraphics.jl, which can consume [`Dataset`](@ref) objects since they implement the Tables interface.
 However, we could just as easily have used StatsPlots.jl.
-"""
 
-# ╔═╡ 7bdff7ba-408c-4893-9857-7cf328301a34
-md"""
+```{julia}
+#| echo: false
+#| code-fold: true
+#| output: false
+using Pkg;
+cd(@__DIR__)
+Pkg.activate(".."); # for reproducibility use the docs environment.
+```
+
+```{julia}
+using ArviZ, ArviZExampleData, DimensionalData, DataFrames, Statistics, AlgebraOfGraphics, CairoMakie
+using AlgebraOfGraphics: density
+set_aog_theme!()
+```
+
 We'll start by loading some draws from an implementation of the non-centered parameterization of the 8 schools model.
 In this parameterization, the model has some sampling issues.
-"""
 
-# ╔═╡ cbd81558-e355-4e7a-ba51-0ba9299cb558
+```{julia}
 idata = load_example_data("centered_eight")
+```
 
-# ╔═╡ 898d83d1-cd3a-47a4-84b1-f2cf0f7bf959
+```{julia}
 idata.posterior
+```
 
-# ╔═╡ da4e989e-a08c-45af-8e55-20ac3ed58745
-md"""
 The plotting functions we'll be using interact with a tabular view of a `Dataset`.
 Let's see what that view looks like for a `Dataset`:
-"""
 
-# ╔═╡ 52ef5ef5-2ea8-4e1a-8cd5-ba3aa5da6d9b
+```{julia}
 df = DataFrame(idata.posterior)
+```
 
-# ╔═╡ e77c2bc3-8c23-446d-9ec4-3a155aea23e9
-md"""
 The tabular view includes dimensions and variables as columns.
 
 When variables with different dimensions are flattened into a tabular form, there's always some duplication of values.
 As a simple case, note that `chain`, `draw`, and `school` all have repeated values in the above table.
 
 In this case, `theta` has the `school` dimension, but `tau` doesn't, so the values of `tau` will be repeated in the table for each value of `school`.
-"""
 
-# ╔═╡ 0d625d42-726d-4cb3-8456-c54e1001df6d
+```{julia}
 df[df.school .== Ref("Choate"), :].tau == df[df.school .== Ref("Deerfield"), :].tau
+```
 
-# ╔═╡ 9b1a7256-c875-4266-b2fb-47c29f2b13c9
-md"""
 In our first example, this will be important.
 
 Here, let's construct a trace plot.
@@ -89,17 +61,15 @@ Besides `idata`, all functions and types in the following cell are defined in Al
 - `visual` specifies how the data should be converted to a plot.
 - `Lines` is a plot type defined in Makie.
 - `draw` takes this combination and plots it.
-"""
 
-# ╔═╡ 0f4483b1-5820-45e2-8812-14550bff69e2
+```{julia}
 draw(
     data(idata.posterior.mu) *
     mapping(:draw, :mu; color=:chain => nonnumeric) *
     visual(Lines; alpha=0.8),
 )
+```
 
-# ╔═╡ ee4ab468-dbf6-46af-a7c5-44b82c031c2c
-md"""
 Note the line `idata.posterior.mu`.
 If we had just used `idata.posterior`, the plot would have looked more-or-less the same, but there would be artifacts due to `mu` being copied many times.
 By selecting `mu` directly, all other dimensions are discarded, so each value of `mu` appears in the plot exactly once.
@@ -107,67 +77,57 @@ By selecting `mu` directly, all other dimensions are discarded, so each value of
 When examining an MCMC trace plot, we want to see a "fuzzy caterpillar".
 Instead we see a few places where the Markov chains froze.
 We can do the same for `theta` as well, but it's more useful here to separate these draws by `school`.
-"""
 
-# ╔═╡ 9b590c11-b15a-4199-8d9f-da8e1735fed2
+```{julia}
 draw(
     data(idata.posterior) *
     mapping(:draw, :theta; layout=:school, color=:chain => nonnumeric) *
     visual(Lines; alpha=0.8),
 )
+```
 
-# ╔═╡ ca18a574-a5c6-4d83-a5d6-b5687f569522
-md"""
 Suppose we want to compare `tau` with `theta` for two different schools.
 To do so, we use `InferenceData`s indexing syntax to subset the data.
-"""
 
-# ╔═╡ 746785a2-c472-467f-973a-d2390ec3e0bb
+```{julia}
 draw(
     data(idata[:posterior, school=At(["Choate", "Deerfield"])]) *
     mapping(:theta, :tau; color=:school) *
     density() *
     visual(Contour; levels=10),
 )
+```
 
-# ╔═╡ e221c4b1-2256-4f97-bfd8-aeaf02fedc1b
-md"""
 We can also compare the density plots constructed from each chain for different schools.
-"""
 
-# ╔═╡ 663e5edb-a751-40e7-96ec-495685453515
+```{julia}
 draw(
     data(idata.posterior) *
     mapping(:theta; layout=:school, color=:chain => nonnumeric) *
     density(),
 )
+```
 
-# ╔═╡ 6e10a67d-045a-4fa9-a7ad-90fa79015ea8
-md"""
 If we want to compare many schools in a single plot, an ECDF plot is more convenient.
-"""
 
-# ╔═╡ 9d8a23e1-4961-44c6-b272-afb8883df1d7
+```{julia}
 draw(
     data(idata.posterior) * mapping(:theta; color=:school => nonnumeric) * visual(ECDFPlot);
     axis=(; ylabel="probability"),
 )
+```
 
-# ╔═╡ 7068ad2c-2260-4afc-b85d-229b2be6c207
-md"""
 So far we've just plotted data from one group, but we often want to combine data from multiple groups in one plot.
 The simplest way to do this is to create the plot out of multiple layers.
 Here we use this approach to plot the observations over the posterior predictive distribution.
-"""
 
-# ╔═╡ c860bb8d-c15c-4ac9-b2b7-143f6f12bf0b
+```{julia}
 draw(
     (data(idata.posterior_predictive) * mapping(:obs; layout=:school) * density()) +
     (data(idata.observed_data) * mapping(:obs, :obs => zero => ""; layout=:school)),
 )
+```
 
-# ╔═╡ 6017222f-769b-4db7-9918-af70c792a158
-md"""
 Another option is to combine the groups into a single dataset.
 
 Here we compare the prior and posterior.
@@ -177,9 +137,8 @@ This is not currently supported.
 We can then either plot the two distributions separately as we did before, or we can compare a single chain from each group.
 This is what we'll do here.
 To concatenate the two groups, we introduce a new named dimension using `DimensionalData.Dim`.
-"""
 
-# ╔═╡ b072dd93-3076-4630-bb2a-b9d490045a78
+```{julia}
 draw(
     data(
         cat(
@@ -191,15 +150,13 @@ draw(
     visual(; alpha=0.8);
     axis=(; ylabel="probability"),
 )
+```
 
-# ╔═╡ 82229ba7-101d-4fbc-96a9-ca56d899ebd5
-md"""
 From the trace plots, we suspected the geometry of this posterior was bad.
 Let's highlight divergent transitions.
 To do so, we merge `posterior` and `samplestats`, which can do with `merge` since they share no common variable names.
-"""
 
-# ╔═╡ 2a47f53a-a054-426d-b536-ccffcf62dd15
+```{julia}
 draw(
     data(merge(idata.posterior, idata.sample_stats)) * mapping(
         :theta,
@@ -209,77 +166,38 @@ draw(
         markersize=:diverging => (d -> d ? 5 : 2),
     ),
 )
+```
 
-# ╔═╡ 3c939c5a-90c6-4367-9a9b-2525796425ce
-md"""
 When we try building more complex plots, we may need to build new `Dataset`s from our existing ones.
 
 One example of this is the corner plot.
 To build this plot, we need to make a copy of `theta` with a copy of the `school` dimension.
-"""
-
-# ╔═╡ b3c8eed9-4083-457e-89dd-71d678b724ef
-let
-    theta = idata.posterior.theta[school=1:4]
-    theta2 = rebuild(set(theta; school=:school2); name=:theta2)
-    plot_data = Dataset(theta, theta2, idata.sample_stats.diverging)
-    draw(
-        data(plot_data) * mapping(
-            :theta,
-            :theta2 => "theta";
-            col=:school,
-            row=:school2,
-            color=:diverging,
-            markersize=:diverging => (d -> d ? 3 : 1),
-        );
-        figure=(; figsize=(5, 5)),
-        axis=(; aspect=1),
-    )
-end
-
-# ╔═╡ 7c12ac28-694c-4c12-af00-2e74bad04683
-md"""
+
+```{julia}
+theta = idata.posterior.theta[school=1:4]
+theta2 = rebuild(set(theta; school=:school2); name=:theta2)
+plot_data = Dataset(theta, theta2, idata.sample_stats.diverging)
+draw(
+    data(plot_data) * mapping(
+        :theta,
+        :theta2 => "theta";
+        col=:school,
+        row=:school2,
+        color=:diverging,
+        markersize=:diverging => (d -> d ? 3 : 1),
+    );
+    figure=(; figsize=(5, 5)),
+    axis=(; aspect=1),
+)
+```
+
 ## Environment
-"""
-
-# ╔═╡ ab01574e-0b29-4325-b42e-9a360267ba95
-with_terminal(Pkg.status; color=false)
-
-# ╔═╡ e3774f8d-9885-4d49-8503-f3fc94b8b113
-with_terminal(versioninfo)
-
-# ╔═╡ Cell order:
-# ╟─70f66296-349c-48a5-9b83-cc5c6cdbd514
-# ╠═b1784f52-0cf2-11ed-32ba-69a8f63b48a9
-# ╟─47c7a1fd-0862-49d5-879d-d59d321ce014
-# ╠═34ee7907-49e7-4e32-9bd5-1c5d490a13d8
-# ╟─7bdff7ba-408c-4893-9857-7cf328301a34
-# ╠═cbd81558-e355-4e7a-ba51-0ba9299cb558
-# ╠═898d83d1-cd3a-47a4-84b1-f2cf0f7bf959
-# ╟─da4e989e-a08c-45af-8e55-20ac3ed58745
-# ╠═52ef5ef5-2ea8-4e1a-8cd5-ba3aa5da6d9b
-# ╟─e77c2bc3-8c23-446d-9ec4-3a155aea23e9
-# ╠═0d625d42-726d-4cb3-8456-c54e1001df6d
-# ╟─9b1a7256-c875-4266-b2fb-47c29f2b13c9
-# ╠═0f4483b1-5820-45e2-8812-14550bff69e2
-# ╟─ee4ab468-dbf6-46af-a7c5-44b82c031c2c
-# ╠═9b590c11-b15a-4199-8d9f-da8e1735fed2
-# ╟─ca18a574-a5c6-4d83-a5d6-b5687f569522
-# ╠═746785a2-c472-467f-973a-d2390ec3e0bb
-# ╟─e221c4b1-2256-4f97-bfd8-aeaf02fedc1b
-# ╠═663e5edb-a751-40e7-96ec-495685453515
-# ╟─6e10a67d-045a-4fa9-a7ad-90fa79015ea8
-# ╠═9d8a23e1-4961-44c6-b272-afb8883df1d7
-# ╟─7068ad2c-2260-4afc-b85d-229b2be6c207
-# ╠═c860bb8d-c15c-4ac9-b2b7-143f6f12bf0b
-# ╟─6017222f-769b-4db7-9918-af70c792a158
-# ╠═b072dd93-3076-4630-bb2a-b9d490045a78
-# ╟─82229ba7-101d-4fbc-96a9-ca56d899ebd5
-# ╠═2a47f53a-a054-426d-b536-ccffcf62dd15
-# ╟─3c939c5a-90c6-4367-9a9b-2525796425ce
-# ╠═b3c8eed9-4083-457e-89dd-71d678b724ef
-# ╟─7c12ac28-694c-4c12-af00-2e74bad04683
-# ╠═e5e15f39-027c-42aa-a765-588ec1cd0b63
-# ╠═41b1ddca-20dc-4821-9cec-d7192ef34795
-# ╠═ab01574e-0b29-4325-b42e-9a360267ba95
-# ╠═e3774f8d-9885-4d49-8503-f3fc94b8b113
+
+```{julia}
+using Pkg
+Pkg.status()
+```
+
+```{julia}
+versioninfo()
+```
diff --git a/docs/src/index.md b/docs/src/index.md
index d2882954..b7269be3 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -25,7 +25,7 @@ pkg> add ArviZ
 
 ## [Usage](@id usage)
 
-See the [Quickstart](./quickstart) for example usage and the [API Overview](@ref api) for description of functions.
+See the [ArviZ Quickstart](@ref) for example usage and the [API Overview](@ref api) for description of functions.
 
 ## [Extending ArviZ.jl](@id extendingarviz)
 
diff --git a/docs/src/quickstart.jl b/docs/src/quickstart.jl
deleted file mode 100644
index 9ceab5f6..00000000
--- a/docs/src/quickstart.jl
+++ /dev/null
@@ -1,432 +0,0 @@
-### A Pluto.jl notebook ###
-# v0.19.16
-
-using Markdown
-using InteractiveUtils
-
-# ╔═╡ 2c0f3d74-8b1e-4648-b0fa-f8050035b0fd
-using Pkg, InteractiveUtils
-
-# ╔═╡ ac57d957-0bdd-457a-ac15-9a4f94f0c785
-# Remove this cell to use release versions of dependencies
-# hideall
-let
-    docs_dir = dirname(@__DIR__)
-    pkg_dir = dirname(docs_dir)
-
-    Pkg.activate(docs_dir)
-    Pkg.develop(; path=pkg_dir)
-    Pkg.instantiate()
-end;
-
-# ╔═╡ 467c2d13-6bfe-4feb-9626-fb14796168aa
-using ArviZ, ArviZPythonPlots, Distributions, LinearAlgebra, Random, StanSample, Turing
-
-# ╔═╡ 56a39a90-0594-48f4-ba04-f7b612019cd1
-using PlutoUI
-
-# ╔═╡ a23dfd65-50a8-4872-8c41-661e96585aca
-md"""
-# [ArviZ Quickstart](#quickstart)
-
-!!! note
-
-    This tutorial is adapted from [ArviZ's quickstart](https://python.arviz.org/en/latest/getting_started/Introduction.html).
-"""
-
-# ╔═╡ d2eedd48-48c6-4fcd-b179-6be7fe68d3d6
-md"""
-## [Setup](#setup)
-
-Here we add the necessary packages for this notebook and load a few we will use throughout.
-"""
-
-# ╔═╡ 06b00794-e97f-472b-b526-efe4815103f8
-# ArviZPythonPlots ships with style sheets!
-use_style("arviz-darkgrid")
-
-# ╔═╡ 5acbfd9a-cfea-4c3c-a3f0-1744eb7e4e27
-md"""
-## [Get started with plotting](#Get-started-with-plotting)
-
-To plot with ArviZ, we need to load the [ArviZPythonPlots](https://julia.arviz.org/ArviZPythonPlots) package.
-ArviZ is designed to be used with libraries like [Stan](https://github.com/StanJulia/Stan.jl), [Turing.jl](https://turinglang.org), and [Soss.jl](https://github.com/cscherrer/Soss.jl) but works fine with raw arrays.
-"""
-
-# ╔═╡ efb3f0af-9fac-48d8-bbb2-2dd6ebd5e4f6
-rng1 = Random.MersenneTwister(37772);
-
-# ╔═╡ 401e9b91-0bca-4369-8d36-3d9f0b3ad60b
-begin
-    plot_posterior(randn(rng1, 100_000))
-    gcf()
-end
-
-# ╔═╡ 2c718ea5-2800-4df6-b62c-e0a9e440a1c3
-md"""
-Plotting a dictionary of arrays, ArviZ will interpret each key as the name of a different random variable.
-Each row of an array is treated as an independent series of draws from the variable, called a _chain_.
-Below, we have 10 chains of 50 draws each for four different distributions.
-"""
-
-# ╔═╡ 49f19c17-ac1d-46b5-a655-4376b7713244
-let
-    s = (50, 10)
-    plot_forest((
-        normal=randn(rng1, s),
-        gumbel=rand(rng1, Gumbel(), s),
-        student_t=rand(rng1, TDist(6), s),
-        exponential=rand(rng1, Exponential(), s),
-    ),)
-    gcf()
-end
-
-# ╔═╡ a9789109-2b90-40f7-926c-e7c87025d15f
-md"""
-## [Plotting with MCMCChains.jl's `Chains` objects produced by Turing.jl](#Plotting-with-MCMCChains.jl's-Chains-objects-produced-by-Turing.jl)
-
-ArviZ is designed to work well with high dimensional, labelled data.
-Consider the [eight schools model](https://statmodeling.stat.columbia.edu/2014/01/21/everything-need-know-bayesian-statistics-learned-eight-schools/), which roughly tries to measure the effectiveness of SAT classes at eight different schools.
-To show off ArviZ's labelling, I give the schools the names of [a different eight schools](https://en.wikipedia.org/wiki/Eight_Schools_Association).
-
-This model is small enough to write down, is hierarchical, and uses labelling.
-Additionally, a centered parameterization causes [divergences](https://mc-stan.org/users/documentation/case-studies/divergences_and_bias.html) (which are interesting for illustration).
-
-First we create our data and set some sampling parameters.
-"""
-
-# ╔═╡ 69124a94-a6aa-43f6-8d4f-fa9a6b1cfef1
-begin
-    J = 8
-    y = [28.0, 8.0, -3.0, 7.0, -1.0, 1.0, 18.0, 12.0]
-    σ = [15.0, 10.0, 16.0, 11.0, 9.0, 11.0, 10.0, 18.0]
-    schools = [
-        "Choate",
-        "Deerfield",
-        "Phillips Andover",
-        "Phillips Exeter",
-        "Hotchkiss",
-        "Lawrenceville",
-        "St. Paul's",
-        "Mt. Hermon",
-    ]
-    ndraws = 1_000
-    ndraws_warmup = 1_000
-    nchains = 4
-end;
-
-# ╔═╡ 10986e0e-55f7-40ea-ab63-274fbab3126d
-md"""
-Now we write and run the model using Turing:
-"""
-
-# ╔═╡ f383d541-e22d-44b4-b8cb-28b3d67944a1
-Turing.@model function model_turing(y, σ, J=length(y))
-    μ ~ Normal(0, 5)
-    τ ~ truncated(Cauchy(0, 5), 0, Inf)
-    θ ~ filldist(Normal(μ, τ), J)
-    for i in 1:J
-        y[i] ~ Normal(θ[i], σ[i])
-    end
-end
-
-# ╔═╡ 86cb5e19-49e4-4e5e-8b89-e76936932055
-rng2 = Random.MersenneTwister(16653);
-
-# ╔═╡ 85bbcba7-c0f9-4c86-9cdf-a27055d3d448
-begin
-    param_mod_turing = model_turing(y, σ)
-    sampler = NUTS(ndraws_warmup, 0.8)
-
-    turing_chns = Turing.sample(
-        rng2, model_turing(y, σ), sampler, MCMCThreads(), ndraws, nchains
-    )
-end;
-
-# ╔═╡ bd4ab044-51ce-4af9-83b2-bd8fc827f810
-md"""
-Most ArviZ functions work fine with `Chains` objects from Turing:
-"""
-
-# ╔═╡ 500f4e0d-0a36-4b5c-8900-667560fbf1d4
-begin
-    plot_autocorr(turing_chns; var_names=(:μ, :τ))
-    gcf()
-end
-
-# ╔═╡ 1129ad94-f65a-4332-b354-21bcf7e53541
-md"""
-### Convert to `InferenceData`
-
-For much more powerful querying, analysis and plotting, we can use built-in ArviZ utilities to convert `Chains` objects to multidimensional data structures with named dimensions and indices.
-Note that for such dimensions, the information is not contained in `Chains`, so we need to provide it.
-
-ArviZ is built to work with [`InferenceData`](https://julia.arviz.org/InferenceObjects/stable/inference_data), and the more *groups* it has access to, the more powerful analyses it can perform.
-"""
-
-# ╔═╡ 803efdd8-656e-4e37-ba36-81195d064972
-idata_turing_post = from_mcmcchains(
-    turing_chns;
-    coords=(; school=schools),
-    dims=NamedTuple(k => (:school,) for k in (:y, :σ, :θ)),
-    library="Turing",
-)
-
-# ╔═╡ 79f342c8-0738-432b-bfd7-2da25e50fa91
-md"""
-Each group is an [`ArviZ.Dataset`](https://julia.arviz.org/InferenceObjects/stable/dataset), a `DimensionalData.AbstractDimStack` that can be used identically to a [`DimensionalData.Dimstack`](https://rafaqz.github.io/DimensionalData.jl/stable/reference/#DimensionalData.DimStack).
-We can view a summary of the dataset.
-"""
-
-# ╔═╡ 6209f947-5001-4507-b3e8-9747256f3328
-idata_turing_post.posterior
-
-# ╔═╡ 26692722-db0d-41de-b58c-339b639f6948
-md"""
-Here is a plot of the trace. Note the intelligent labels.
-"""
-
-# ╔═╡ 14046b83-9c0a-4d33-ae4e-36c7d6f1b2e6
-begin
-    plot_trace(idata_turing_post)
-    gcf()
-end
-
-# ╔═╡ 737f319c-1ddd-45f2-8d10-aaecdc1334be
-md"We can also generate summary stats..."
-
-# ╔═╡ 5be3131d-db8b-44c8-9fcc-571d68695148
-summarystats(idata_turing_post)
-
-# ╔═╡ 0787a842-63ca-407d-9d07-a3d949940f92
-md"...and examine the energy distribution of the Hamiltonian sampler."
-
-# ╔═╡ 6e8343c8-bee3-4e1d-82d6-1885bfd1dbec
-begin
-    plot_energy(idata_turing_post)
-    gcf()
-end
-
-# ╔═╡ cba6f6c9-82c4-4957-acc3-36e9f1c95d76
-md"""
-### Additional information in Turing.jl
-
-With a few more steps, we can use Turing to compute additional useful groups to add to the `InferenceData`.
-
-To sample from the prior, one simply calls `sample` but with the `Prior` sampler:
-"""
-
-# ╔═╡ 79374a8b-98c1-443c-acdd-0ee00bd42b38
-prior = Turing.sample(rng2, param_mod_turing, Prior(), ndraws);
-
-# ╔═╡ c23562ba-3b0b-49d9-bb41-2fedaa9e9500
-md"""
-To draw from the prior and posterior predictive distributions we can instantiate a "predictive model", i.e. a Turing model but with the observations set to `missing`, and then calling `predict` on the predictive model and the previously drawn samples:
-"""
-
-# ╔═╡ d0b447fe-26b2-48e0-bcec-d9f04697973a
-begin
-    # Instantiate the predictive model
-    param_mod_predict = model_turing(similar(y, Missing), σ)
-    # and then sample!
-    prior_predictive = Turing.predict(rng2, param_mod_predict, prior)
-    posterior_predictive = Turing.predict(rng2, param_mod_predict, turing_chns)
-end;
-
-# ╔═╡ 4d2fdcbe-c1d4-43b6-b382-f2e956b952a1
-md"""
-And to extract the pointwise log-likelihoods, which is useful if you want to compute metrics such as [`loo`](https://julia.arviz.org/ArviZ/stable/api/stats/#PosteriorStats.loo),
-"""
-
-# ╔═╡ 5a075722-232f-40fc-a499-8dc5b0c2424a
-log_likelihood = let
-    log_likelihood = Turing.pointwise_loglikelihoods(
-        param_mod_turing, MCMCChains.get_sections(turing_chns, :parameters)
-    )
-    # Ensure the ordering of the loglikelihoods matches the ordering of `posterior_predictive`
-    ynames = string.(keys(posterior_predictive))
-    log_likelihood_y = getindex.(Ref(log_likelihood), ynames)
-    (; y=cat(log_likelihood_y...; dims=3))
-end;
-
-# ╔═╡ 1b5af2c3-f2ce-4e9d-9ad7-ac287a9178e2
-md"This can then be included in the [`from_mcmcchains`](https://julia.arviz.org/ArviZ/stable/api/data/#ArviZ.from_mcmcchains) call from above:"
-
-# ╔═╡ b38c7a43-f00c-43c0-aa6b-9c581d6d0c73
-idata_turing = from_mcmcchains(
-    turing_chns;
-    posterior_predictive,
-    log_likelihood,
-    prior,
-    prior_predictive,
-    observed_data=(; y),
-    coords=(; school=schools),
-    dims=NamedTuple(k => (:school,) for k in (:y, :σ, :θ)),
-    library=Turing,
-)
-
-# ╔═╡ a3b71e4e-ac1f-4404-b8a2-35d03d326774
-md"""
-Then we can for example compute the expected *leave-one-out (LOO)* predictive density, which is an estimate of the out-of-distribution predictive fit of the model:
-"""
-
-# ╔═╡ f552b5b5-9744-41df-af90-46405367ea0b
-loo(idata_turing) # higher ELPD is better
-
-# ╔═╡ 9d3673f5-b57b-432e-944a-70b23643128a
-md"""
-If the model is well-calibrated, i.e. it replicates the true generative process well, the CDF of the pointwise LOO values should be similarly distributed to a uniform distribution.
-This can be inspected visually:
-"""
-
-# ╔═╡ 05c9be29-7758-4324-971c-5579f99aaf9d
-begin
-    plot_loo_pit(idata_turing; y=:y, ecdf=true)
-    gcf()
-end
-
-# ╔═╡ 98acc304-22e3-4e6b-a2f4-d22f6847145b
-md"""
-## [Plotting with Stan.jl outputs](#Plotting-with-Stan.jl-outputs)
-
-StanSample.jl comes with built-in support for producing `InferenceData` outputs.
-
-Here is the same centered eight schools model in Stan:
-"""
-
-# ╔═╡ b46af168-1ce3-4058-a014-b66c645a6e0d
-begin
-    schools_code = """
-    data {
-      int<lower=0> J;
-      real y[J];
-      real<lower=0> sigma[J];
-    }
-
-    parameters {
-      real mu;
-      real<lower=0> tau;
-      real theta[J];
-    }
-
-    model {
-      mu ~ normal(0, 5);
-      tau ~ cauchy(0, 5);
-      theta ~ normal(mu, tau);
-      y ~ normal(theta, sigma);
-    }
-
-    generated quantities {
-        vector[J] log_lik;
-        vector[J] y_hat;
-        for (j in 1:J) {
-            log_lik[j] = normal_lpdf(y[j] | theta[j], sigma[j]);
-            y_hat[j] = normal_rng(theta[j], sigma[j]);
-        }
-    }
-    """
-
-    schools_data = Dict("J" => J, "y" => y, "sigma" => σ)
-    idata_stan = mktempdir() do path
-        stan_model = SampleModel("schools", schools_code, path)
-        _ = stan_sample(
-            stan_model;
-            data=schools_data,
-            num_chains=nchains,
-            num_warmups=ndraws_warmup,
-            num_samples=ndraws,
-            seed=28983,
-            summary=false,
-        )
-        return StanSample.inferencedata(
-            stan_model;
-            posterior_predictive_var=:y_hat,
-            observed_data=(; y),
-            log_likelihood_var=:log_lik,
-            coords=(; school=schools),
-            dims=NamedTuple(
-                k => (:school,) for k in (:y, :sigma, :theta, :log_lik, :y_hat)
-            ),
-        )
-    end
-end
-
-# ╔═╡ ab145e41-b230-4cad-bef5-f31e0e0770d4
-begin
-    plot_density(idata_stan; var_names=(:mu, :tau))
-    gcf()
-end
-
-# ╔═╡ e44b260c-9d2f-43f8-a64b-04245a0a5658
-md"""Here is a plot showing where the Hamiltonian sampler had divergences:"""
-
-# ╔═╡ 5070bbbc-68d2-49b8-bd91-456dc0da4573
-begin
-    plot_pair(
-        idata_stan;
-        coords=Dict(:school => ["Choate", "Deerfield", "Phillips Andover"]),
-        divergences=true,
-    )
-    gcf()
-end
-
-# ╔═╡ ac2b4378-bd1c-4164-af05-d9a35b1bb08f
-md"## [Environment](#environment)"
-
-# ╔═╡ fd84237a-7a19-4bbb-a702-faa31075ecbc
-with_terminal(Pkg.status; color=false)
-
-# ╔═╡ ad29b7f3-6f5e-4b04-bf70-2308a7d110c0
-with_terminal(versioninfo)
-
-# ╔═╡ Cell order:
-# ╟─a23dfd65-50a8-4872-8c41-661e96585aca
-# ╟─d2eedd48-48c6-4fcd-b179-6be7fe68d3d6
-# ╠═ac57d957-0bdd-457a-ac15-9a4f94f0c785
-# ╠═467c2d13-6bfe-4feb-9626-fb14796168aa
-# ╠═06b00794-e97f-472b-b526-efe4815103f8
-# ╟─5acbfd9a-cfea-4c3c-a3f0-1744eb7e4e27
-# ╠═efb3f0af-9fac-48d8-bbb2-2dd6ebd5e4f6
-# ╠═401e9b91-0bca-4369-8d36-3d9f0b3ad60b
-# ╟─2c718ea5-2800-4df6-b62c-e0a9e440a1c3
-# ╠═49f19c17-ac1d-46b5-a655-4376b7713244
-# ╟─a9789109-2b90-40f7-926c-e7c87025d15f
-# ╠═69124a94-a6aa-43f6-8d4f-fa9a6b1cfef1
-# ╟─10986e0e-55f7-40ea-ab63-274fbab3126d
-# ╠═f383d541-e22d-44b4-b8cb-28b3d67944a1
-# ╠═86cb5e19-49e4-4e5e-8b89-e76936932055
-# ╠═85bbcba7-c0f9-4c86-9cdf-a27055d3d448
-# ╟─bd4ab044-51ce-4af9-83b2-bd8fc827f810
-# ╠═500f4e0d-0a36-4b5c-8900-667560fbf1d4
-# ╟─1129ad94-f65a-4332-b354-21bcf7e53541
-# ╠═803efdd8-656e-4e37-ba36-81195d064972
-# ╟─79f342c8-0738-432b-bfd7-2da25e50fa91
-# ╠═6209f947-5001-4507-b3e8-9747256f3328
-# ╟─26692722-db0d-41de-b58c-339b639f6948
-# ╠═14046b83-9c0a-4d33-ae4e-36c7d6f1b2e6
-# ╟─737f319c-1ddd-45f2-8d10-aaecdc1334be
-# ╠═5be3131d-db8b-44c8-9fcc-571d68695148
-# ╟─0787a842-63ca-407d-9d07-a3d949940f92
-# ╠═6e8343c8-bee3-4e1d-82d6-1885bfd1dbec
-# ╟─cba6f6c9-82c4-4957-acc3-36e9f1c95d76
-# ╠═79374a8b-98c1-443c-acdd-0ee00bd42b38
-# ╟─c23562ba-3b0b-49d9-bb41-2fedaa9e9500
-# ╠═d0b447fe-26b2-48e0-bcec-d9f04697973a
-# ╟─4d2fdcbe-c1d4-43b6-b382-f2e956b952a1
-# ╠═5a075722-232f-40fc-a499-8dc5b0c2424a
-# ╟─1b5af2c3-f2ce-4e9d-9ad7-ac287a9178e2
-# ╠═b38c7a43-f00c-43c0-aa6b-9c581d6d0c73
-# ╟─a3b71e4e-ac1f-4404-b8a2-35d03d326774
-# ╠═f552b5b5-9744-41df-af90-46405367ea0b
-# ╟─9d3673f5-b57b-432e-944a-70b23643128a
-# ╠═05c9be29-7758-4324-971c-5579f99aaf9d
-# ╟─98acc304-22e3-4e6b-a2f4-d22f6847145b
-# ╠═b46af168-1ce3-4058-a014-b66c645a6e0d
-# ╠═ab145e41-b230-4cad-bef5-f31e0e0770d4
-# ╟─e44b260c-9d2f-43f8-a64b-04245a0a5658
-# ╠═5070bbbc-68d2-49b8-bd91-456dc0da4573
-# ╟─ac2b4378-bd1c-4164-af05-d9a35b1bb08f
-# ╠═56a39a90-0594-48f4-ba04-f7b612019cd1
-# ╠═2c0f3d74-8b1e-4648-b0fa-f8050035b0fd
-# ╠═fd84237a-7a19-4bbb-a702-faa31075ecbc
-# ╠═ad29b7f3-6f5e-4b04-bf70-2308a7d110c0
diff --git a/docs/src/quickstart.qmd b/docs/src/quickstart.qmd
new file mode 100644
index 00000000..e9bfc9c5
--- /dev/null
+++ b/docs/src/quickstart.qmd
@@ -0,0 +1,305 @@
+---
+title: "ArviZ Quickstart"
+---
+
+## Set-up
+
+Here we add the necessary packages for this notebook and load a few we will use throughout.
+
+```{julia}
+#| echo: false
+#| code-fold: true
+#| output: false
+using Pkg;
+cd(@__DIR__)
+Pkg.activate(".."); # for reproducibility use the docs environment.
+```
+
+```{julia}
+using ArviZ, ArviZPythonPlots, Distributions, LinearAlgebra, Random, StanSample, Turing
+```
+
+```{julia}
+# ArviZPythonPlots ships with style sheets!
+use_style("arviz-darkgrid")
+```
+
+## Get started with plotting
+
+To plot with ArviZ, we need to load the [ArviZPythonPlots](https://julia.arviz.org/ArviZPythonPlots) package.
+ArviZ is designed to be used with libraries like [Stan](https://github.com/StanJulia/Stan.jl), [Turing.jl](https://turinglang.org), and [Soss.jl](https://github.com/cscherrer/Soss.jl) but works fine with raw arrays.
+
+```{julia}
+rng1 = Random.MersenneTwister(37772);
+```
+
+```{julia}
+plot_posterior(randn(rng1, 100_000));
+```
+
+Plotting a dictionary of arrays, ArviZ will interpret each key as the name of a different random variable.
+Each row of an array is treated as an independent series of draws from the variable, called a _chain_.
+Below, we have 10 chains of 50 draws each for four different distributions.
+
+```{julia}
+s = (50, 10)
+plot_forest((
+    normal=randn(rng1, s),
+    gumbel=rand(rng1, Gumbel(), s),
+    student_t=rand(rng1, TDist(6), s),
+    exponential=rand(rng1, Exponential(), s),
+));
+```
+
+## Plotting with MCMCChains.jl's `Chains` objects produced by Turing.jl
+
+ArviZ is designed to work well with high dimensional, labelled data.
+Consider the [eight schools model](https://statmodeling.stat.columbia.edu/2014/01/21/everything-need-know-bayesian-statistics-learned-eight-schools/), which roughly tries to measure the effectiveness of SAT classes at eight different schools.
+To show off ArviZ's labelling, I give the schools the names of [a different eight schools](https://en.wikipedia.org/wiki/Eight_Schools_Association).
+
+This model is small enough to write down, is hierarchical, and uses labelling.
+Additionally, a centered parameterization causes [divergences](https://mc-stan.org/users/documentation/case-studies/divergences_and_bias.html) (which are interesting for illustration).
+
+First we create our data and set some sampling parameters.
+
+```{julia}
+J = 8
+y = [28.0, 8.0, -3.0, 7.0, -1.0, 1.0, 18.0, 12.0]
+σ = [15.0, 10.0, 16.0, 11.0, 9.0, 11.0, 10.0, 18.0]
+schools = [
+    "Choate",
+    "Deerfield",
+    "Phillips Andover",
+    "Phillips Exeter",
+    "Hotchkiss",
+    "Lawrenceville",
+    "St. Paul's",
+    "Mt. Hermon",
+]
+ndraws = 1_000
+ndraws_warmup = 1_000
+nchains = 4;
+```
+
+Now we write and run the model using Turing:
+
+```{julia}
+Turing.@model function model_turing(y, σ, J=length(y))
+    μ ~ Normal(0, 5)
+    τ ~ truncated(Cauchy(0, 5), 0, Inf)
+    θ ~ filldist(Normal(μ, τ), J)
+    for i in 1:J
+        y[i] ~ Normal(θ[i], σ[i])
+    end
+end;
+```
+
+```{julia}
+rng2 = Random.MersenneTwister(16653);
+```
+
+```{julia}
+param_mod_turing = model_turing(y, σ)
+sampler = NUTS(ndraws_warmup, 0.8)
+
+turing_chns = Turing.sample(
+    rng2, model_turing(y, σ), sampler, MCMCThreads(), ndraws, nchains
+);
+```
+
+Most ArviZ functions work fine with `Chains` objects from Turing:
+
+```{julia}
+plot_autocorr(turing_chns; var_names=(:μ, :τ));
+```
+
+### Convert to `InferenceData`
+
+For much more powerful querying, analysis and plotting, we can use built-in ArviZ utilities to convert `Chains` objects to multidimensional data structures with named dimensions and indices.
+Note that for such dimensions, the information is not contained in `Chains`, so we need to provide it.
+
+ArviZ is built to work with [`InferenceData`](@ref), and the more *groups* it has access to, the more powerful analyses it can perform.
+
+```{julia}
+idata_turing_post = from_mcmcchains(
+    turing_chns;
+    coords=(; school=schools),
+    dims=NamedTuple(k => (:school,) for k in (:y, :σ, :θ)),
+    library="Turing",
+)
+```
+
+Each group is a [`Dataset`](@ref), a `DimensionalData.AbstractDimStack` that can be used identically to a [`DimensionalData.Dimstack`](https://rafaqz.github.io/DimensionalData.jl/stable/reference/#DimensionalData.DimStack).
+We can view a summary of the dataset.
+
+```{julia}
+idata_turing_post.posterior
+```
+
+Here is a plot of the trace. Note the intelligent labels.
+
+```{julia}
+plot_trace(idata_turing_post);
+```
+
+We can also generate summary stats...
+
+```{julia}
+summarystats(idata_turing_post)
+```
+
+...and examine the energy distribution of the Hamiltonian sampler.
+
+```{julia}
+plot_energy(idata_turing_post);
+```
+
+### Additional information in Turing.jl
+
+With a few more steps, we can use Turing to compute additional useful groups to add to the `InferenceData`.
+
+To sample from the prior, one simply calls `sample` but with the `Prior` sampler:
+
+```{julia}
+prior = Turing.sample(rng2, param_mod_turing, Prior(), ndraws);
+```
+
+To draw from the prior and posterior predictive distributions we can instantiate a "predictive model", i.e. a Turing model but with the observations set to `missing`, and then calling `predict` on the predictive model and the previously drawn samples:
+
+```{julia}
+# Instantiate the predictive model
+param_mod_predict = model_turing(similar(y, Missing), σ)
+# and then sample!
+prior_predictive = Turing.predict(rng2, param_mod_predict, prior)
+posterior_predictive = Turing.predict(rng2, param_mod_predict, turing_chns);
+```
+
+And to extract the pointwise log-likelihoods, which is useful if you want to compute metrics such as [`loo`](@ref),
+
+```{julia}
+log_likelihood = let
+    log_likelihood = Turing.pointwise_loglikelihoods(
+        param_mod_turing, MCMCChains.get_sections(turing_chns, :parameters)
+    )
+    # Ensure the ordering of the loglikelihoods matches the ordering of `posterior_predictive`
+    ynames = string.(keys(posterior_predictive))
+    log_likelihood_y = getindex.(Ref(log_likelihood), ynames)
+    (; y=cat(log_likelihood_y...; dims=3))
+end;
+```
+
+This can then be included in the [`from_mcmcchains`](@ref) call from above:
+
+```{julia}
+idata_turing = from_mcmcchains(
+    turing_chns;
+    posterior_predictive,
+    log_likelihood,
+    prior,
+    prior_predictive,
+    observed_data=(; y),
+    coords=(; school=schools),
+    dims=NamedTuple(k => (:school,) for k in (:y, :σ, :θ)),
+    library=Turing,
+)
+```
+
+Then we can for example compute the expected *leave-one-out (LOO)* predictive density, which is an estimate of the out-of-distribution predictive fit of the model:
+
+```{julia}
+loo(idata_turing) # higher ELPD is better
+```
+
+If the model is well-calibrated, i.e. it replicates the true generative process well, the CDF of the pointwise LOO values should be similarly distributed to a uniform distribution.
+This can be inspected visually:
+
+```{julia}
+plot_loo_pit(idata_turing; y=:y, ecdf=true);
+```
+
+## Plotting with Stan.jl outputs
+
+StanSample.jl comes with built-in support for producing `InferenceData` outputs.
+
+Here is the same centered eight schools model in Stan:
+
+```{julia}
+schools_code = """
+data {
+    int<lower=0> J;
+    real y[J];
+    real<lower=0> sigma[J];
+}
+
+parameters {
+    real mu;
+    real<lower=0> tau;
+    real theta[J];
+}
+
+model {
+    mu ~ normal(0, 5);
+    tau ~ cauchy(0, 5);
+    theta ~ normal(mu, tau);
+    y ~ normal(theta, sigma);
+}
+
+generated quantities {
+    vector[J] log_lik;
+    vector[J] y_hat;
+    for (j in 1:J) {
+        log_lik[j] = normal_lpdf(y[j] | theta[j], sigma[j]);
+        y_hat[j] = normal_rng(theta[j], sigma[j]);
+    }
+}
+"""
+
+schools_data = Dict("J" => J, "y" => y, "sigma" => σ)
+idata_stan = mktempdir() do path
+    stan_model = SampleModel("schools", schools_code, path)
+    _ = stan_sample(
+        stan_model;
+        data=schools_data,
+        num_chains=nchains,
+        num_warmups=ndraws_warmup,
+        num_samples=ndraws,
+        seed=28983,
+        summary=false,
+    )
+    return StanSample.inferencedata(
+        stan_model;
+        posterior_predictive_var=:y_hat,
+        observed_data=(; y),
+        log_likelihood_var=:log_lik,
+        coords=(; school=schools),
+        dims=NamedTuple(
+            k => (:school,) for k in (:y, :sigma, :theta, :log_lik, :y_hat)
+        ),
+    )
+end
+```
+
+```{julia}
+plot_density(idata_stan; var_names=(:mu, :tau));
+```
+
+Here is a plot showing where the Hamiltonian sampler had divergences:
+
+```{julia}
+plot_pair(
+    idata_stan;
+    coords=Dict(:school => ["Choate", "Deerfield", "Phillips Andover"]),
+    divergences=true,
+);
+```
+
+## Environment
+
+```{julia}
+using Pkg
+Pkg.status()
+```
+
+```{julia}
+versioninfo()
+```