diff --git a/docs/src/pma.md b/docs/src/pma.md
index 8a978b3f..bf2e2ed6 100644
--- a/docs/src/pma.md
+++ b/docs/src/pma.md
@@ -211,12 +211,6 @@ pairplot(
             PairPlots.MarginConfidenceLimits()
         ),
     labels=Dict(:mass=>"mass [Mⱼᵤₚ]", :a=>"sma. [au]"),
-    axis = (;
-        a = (;
-            scale=Makie.pseudolog10,
-            ticks=2 .^ (0:1:6)
-        )
-    )
 )
 ```
 
diff --git a/ext/OctofitterPigeonsExt.jl b/ext/OctofitterPigeonsExt.jl
index 481546b9..19889205 100644
--- a/ext/OctofitterPigeonsExt.jl
+++ b/ext/OctofitterPigeonsExt.jl
@@ -13,6 +13,15 @@ function Pigeons.initialization(model::Octofitter.LogDensityModel, rng::Abstract
     return initial_θ_t
 end
 
+# Valid for reference model only
+function Pigeons.sample_iid!(model_reference::Octofitter.LogDensityModel, replica, shared)
+    # This could in theory be done without any array allocations
+    θ = sample_priors(replica.rng, model_reference.system)
+    θ_t = model_reference.link(θ)
+    replica.state .= θ_t
+end
+
+
 
 """
 octofit_pigeons(model; nrounds, n_chains=[auto])