diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
new file mode 100644
index 0000000..0d43ced
--- /dev/null
+++ b/.JuliaFormatter.toml
@@ -0,0 +1,7 @@
+
+style = "blue"
+align_assignment = true
+align_struct_field = true
+align_pair_arrow = true
+align_matrix = true
+align_conditional = true
diff --git a/.github/workflows/Format.yml b/.github/workflows/Format.yml
new file mode 100644
index 0000000..732212c
--- /dev/null
+++ b/.github/workflows/Format.yml
@@ -0,0 +1,26 @@
+name: Format suggestions
+
+on:
+  pull_request:
+  
+concurrency:
+  # Skip intermediate builds: always.
+  # Cancel intermediate builds: only if it is a pull request build.
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }}
+
+jobs:
+  format:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: 1
+      - run: |
+          julia  -e 'using Pkg; Pkg.add("JuliaFormatter")'
+          julia  -e 'using JuliaFormatter; format("."; verbose=true)'
+      - uses: reviewdog/action-suggester@v1
+        with:
+          tool_name: JuliaFormatter
+          fail_on_error: true
diff --git a/docs/make.jl b/docs/make.jl
index 3a27f0d..3fddeb7 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -20,11 +20,8 @@ makedocs(;
         "Univariate Slice Sampling"     => "univariate_slice.md",
         "Meta Multivariate Samplers"    => "meta_multivariate.md",
         "Latent Slice Sampling"         => "latent_slice.md",
-        "Gibbsian Polar Slice Sampling" => "gibbs_polar.md"
+        "Gibbsian Polar Slice Sampling" => "gibbs_polar.md",
     ],
 )
 
-deploydocs(;
-    repo="github.com/TuringLang/SliceSampling.jl",
-    push_preview=true
-)
+deploydocs(; repo="github.com/TuringLang/SliceSampling.jl", push_preview=true)
diff --git a/ext/SliceSamplingTuringExt.jl b/ext/SliceSamplingTuringExt.jl
index 213c281..fd3e6ab 100644
--- a/ext/SliceSamplingTuringExt.jl
+++ b/ext/SliceSamplingTuringExt.jl
@@ -6,7 +6,7 @@ if isdefined(Base, :get_extension)
     using Random
     using SliceSampling
     using Turing
-   # using Turing: Turing, Experimental
+    # using Turing: Turing, Experimental
 else
     using ..LogDensityProblemsAD
     using ..Random
@@ -17,46 +17,47 @@ end
 
 # Required for using the slice samplers as `externalsampler`s in Turing
 # begin
-Turing.Inference.getparams(
-          ::Turing.DynamicPPL.Model,
-    sample::SliceSampling.Transition
-) = sample.params
+function Turing.Inference.getparams(
+    ::Turing.DynamicPPL.Model, sample::SliceSampling.Transition
+)
+    return sample.params
+end
 # end
 
 # Required for using the slice samplers as `Experimental.Gibbs` samplers in Turing
 # begin
-Turing.Inference.getparams(
-         ::Turing.DynamicPPL.Model,
-    state::SliceSampling.UnivariateSliceState
-) = state.transition.params
+function Turing.Inference.getparams(
+    ::Turing.DynamicPPL.Model, state::SliceSampling.UnivariateSliceState
+)
+    return state.transition.params
+end
 
-Turing.Inference.getparams(
-         ::Turing.DynamicPPL.Model,
-    state::SliceSampling.GibbsState
-) = state.transition.params
+function Turing.Inference.getparams(
+    ::Turing.DynamicPPL.Model, state::SliceSampling.GibbsState
+)
+    return state.transition.params
+end
 
-Turing.Inference.getparams(
-         ::Turing.DynamicPPL.Model,
-    state::SliceSampling.HitAndRunState
-) = state.transition.params
+function Turing.Inference.getparams(
+    ::Turing.DynamicPPL.Model, state::SliceSampling.HitAndRunState
+)
+    return state.transition.params
+end
 
-Turing.Experimental.gibbs_requires_recompute_logprob(
+function Turing.Experimental.gibbs_requires_recompute_logprob(
     model_dst,
     ::Turing.DynamicPPL.Sampler{
-        <: Turing.Inference.ExternalSampler{
-            <: SliceSampling.AbstractSliceSampling, A, U
-        }
+        <:Turing.Inference.ExternalSampler{<:SliceSampling.AbstractSliceSampling,A,U}
     },
     sampler_src,
     state_dst,
-    state_src
-) where {A,U} = false
+    state_src,
+) where {A,U}
+    return false
+end
 # end
 
-function SliceSampling.initial_sample(
-    rng::Random.AbstractRNG,
-    ℓ  ::Turing.LogDensityFunction
-)
+function SliceSampling.initial_sample(rng::Random.AbstractRNG, ℓ::Turing.LogDensityFunction)
     model = ℓ.model
     spl   = Turing.SampleFromUniform()
     vi    = Turing.VarInfo(rng, model, spl)
@@ -67,14 +68,14 @@ function SliceSampling.initial_sample(
         if init_attempt_count == 10
             @warn "failed to find valid initial parameters in $(init_attempt_count) tries; consider providing explicit initial parameters using the `initial_params` keyword"
         end
-        
+
         # NOTE: This will sample in the unconstrained space.
         vi = last(DynamicPPL.evaluate!!(model, rng, vi, SampleFromUniform()))
         θ  = vi[spl]
 
         init_attempt_count += 1
     end
-    θ
+    return θ
 end
 
 end
diff --git a/src/SliceSampling.jl b/src/SliceSampling.jl
index cb1c82e..a795a16 100644
--- a/src/SliceSampling.jl
+++ b/src/SliceSampling.jl
@@ -32,7 +32,7 @@ Struct containing the results of the transition.
 - `lp::Real`: Log-target density of the samples.
 - `info::NamedTuple`: Named tuple containing information about the transition. 
 """
-struct Transition{P, L <: Real, I <: NamedTuple}
+struct Transition{P,L<:Real,I<:NamedTuple}
     "current state of the slice sampling chain"
     params::P
 
@@ -53,47 +53,44 @@ Return the initial sample for the `model` using the random number generator `rng
 - `model`: The target `LogDensityProblem`.
 """
 function initial_sample(::Random.AbstractRNG, ::Any)
-    error(
+    return error(
         "`initial_sample` is not implemented but an initialization wasn't provided. ",
-        "Consider supplying an initialization to `initial_params`."
+        "Consider supplying an initialization to `initial_params`.",
     )
 end
 
 # If target is from `LogDensityProblemsAD`, unwrap target before calling `initial_sample`.
 # This is necessary since Turing wraps `DynamicPPL.Model`s when passed to an `externalsampler`.
-initial_sample(
-    rng::Random.AbstractRNG,
-    wrap::LogDensityProblemsAD.ADGradientWrapper
-) = initial_sample(rng, parent(wrap))
+function initial_sample(
+    rng::Random.AbstractRNG, wrap::LogDensityProblemsAD.ADGradientWrapper
+)
+    return initial_sample(rng, parent(wrap))
+end
 
 function exceeded_max_prop(max_prop::Int)
-    error("Exceeded maximum number of proposal $(max_prop), ",
-          "which indicates an acceptance rate less than $(1/max_prop*100)%. ", 
-          "A quick fix is to increase `max_prop`, ",
-          "but an acceptance rate that is too low often indicates that there is a problem. ",
-          "Here are some possible causes:\n",
-          "- The model might be broken or degenerate (most likely cause).\n",
-          "- The tunable parameters of the sampler are suboptimal.\n",
-          "- The initialization is pathologic. (try supplying a (different) `initial_params`)\n",
-          "- There might be a bug in the sampler. (if this is suspected, file an issue to `SliceSampling`)\n"
-          )
+    return error(
+        "Exceeded maximum number of proposal $(max_prop), ",
+        "which indicates an acceptance rate less than $(1/max_prop*100)%. ",
+        "A quick fix is to increase `max_prop`, ",
+        "but an acceptance rate that is too low often indicates that there is a problem. ",
+        "Here are some possible causes:\n",
+        "- The model might be broken or degenerate (most likely cause).\n",
+        "- The tunable parameters of the sampler are suboptimal.\n",
+        "- The initialization is pathologic. (try supplying a (different) `initial_params`)\n",
+        "- There might be a bug in the sampler. (if this is suspected, file an issue to `SliceSampling`)\n",
+    )
 end
 
 ## Univariate Slice Sampling Algorithms
 export Slice, SliceSteppingOut, SliceDoublingOut
 
-abstract type AbstractUnivariateSliceSampling <: AbstractSliceSampling  end
+abstract type AbstractUnivariateSliceSampling <: AbstractSliceSampling end
 
-accept_slice_proposal(
-    ::AbstractSliceSampling,
-    ::Any,
-    ::Real,
-    ::Real,
-    ::Real,
-    ::Real,
-    ::Real,
-    ::Real,
-) = true
+function accept_slice_proposal(
+    ::AbstractSliceSampling, ::Any, ::Real, ::Real, ::Real, ::Real, ::Real, ::Real
+)
+    return true
+end
 
 function find_interval end
 
@@ -103,7 +100,7 @@ include("univariate/steppingout.jl")
 include("univariate/doublingout.jl")
 
 ## Multivariate slice sampling algorithms
-abstract type AbstractMultivariateSliceSampling <: AbstractSliceSampling  end
+abstract type AbstractMultivariateSliceSampling <: AbstractSliceSampling end
 
 # Meta Multivariate Samplers
 export RandPermGibbs, HitAndRun
diff --git a/src/multivariate/gibbspolar.jl b/src/multivariate/gibbspolar.jl
index 9600b03..f5be0d3 100644
--- a/src/multivariate/gibbspolar.jl
+++ b/src/multivariate/gibbspolar.jl
@@ -23,16 +23,18 @@ Gibbsian polar slice sampling algorithm by P. Schär, M. Habeck, and D. Rudolf [
 !!! info
     For Turing users: `Turing` might change `initial_params` to match the support of the posterior. This might lead to \$\$\\lVert x_0 \\rVert\$\$ being small, even though the vector you passed to`initial_params` has a sufficiently large norm. If this is suspected, simply try a different initialization value.
 """
-struct GibbsPolarSlice{W <: Real} <: AbstractMultivariateSliceSampling
+struct GibbsPolarSlice{W<:Real} <: AbstractMultivariateSliceSampling
     w::W
     max_proposals::Int
 end
 
-GibbsPolarSlice(w::Real; max_proposals::Int = DEFAULT_MAX_PROPOSALS) = GibbsPolarSlice(w, max_proposals)
+function GibbsPolarSlice(w::Real; max_proposals::Int=DEFAULT_MAX_PROPOSALS)
+    return GibbsPolarSlice(w, max_proposals)
+end
 
-struct GibbsPolarSliceState{T <: Transition, R <: Real, D <: AbstractVector} 
+struct GibbsPolarSliceState{T<:Transition,R<:Real,D<:AbstractVector}
     "Current [`Transition`](@ref)."
-    transition ::T
+    transition::T
 
     "direction (\$\\theta\$ in the original paper[^SHR2023])"
     direction::D
@@ -47,19 +49,21 @@ end
 
 function logdensity(target::GibbsPolarSliceTarget, x)
     d = length(x)
-    (d-1)*log(norm(x)) + LogDensityProblems.logdensity(target.model, x)
+    return (d - 1) * log(norm(x)) + LogDensityProblems.logdensity(target.model, x)
 end
 
-function AbstractMCMC.step(rng    ::Random.AbstractRNG,
-                           model  ::AbstractMCMC.LogDensityModel,
-                           sampler::GibbsPolarSlice;
-                           initial_params = nothing,
-                           kwargs...)
+function AbstractMCMC.step(
+    rng::Random.AbstractRNG,
+    model::AbstractMCMC.LogDensityModel,
+    sampler::GibbsPolarSlice;
+    initial_params=nothing,
+    kwargs...,
+)
     logdensitymodel = model.logdensity
-    x  = initial_params === nothing ? initial_sample(rng, logdensitymodel) : initial_params
-    d  = length(x)
+    x = initial_params === nothing ? initial_sample(rng, logdensitymodel) : initial_params
+    d = length(x)
     @assert d ≥ 2 "Gibbsian polar slice sampling works reliably only in dimension ≥2"
-    r  = norm(x)
+    r = norm(x)
     if r < 1e-5
         @warn "The norm of initial_params is smaller than 1e-5, which might be result in unstable behavior and the sampler might even get stuck indefinitely. If you are using Turing, this might be due to change of support through Bijectors."
     end
@@ -72,30 +76,30 @@ end
 function rand_subsphere(rng::Random.AbstractRNG, θ::AbstractVector)
     d  = length(θ)
     V1 = randn(rng, eltype(θ), d)
-    V2 = V1 - dot(θ, V1)*θ
-    V2 / max(norm(V2), eps(eltype(θ)))
+    V2 = V1 - dot(θ, V1) * θ
+    return V2 / max(norm(V2), eps(eltype(θ)))
 end
 
 function geodesic_shrinkage(
-    rng     ::Random.AbstractRNG,
-    ϱ1      ::GibbsPolarSliceTarget,
-    ℓT      ::F,
-    θ       ::AbstractVector{F},
-    r       ::F,
-    max_prop::Int
-) where {F <: Real}
-    y     = rand_subsphere(rng, θ)   
-    ω_max = convert(F, 2π)*rand(rng, F)
+    rng::Random.AbstractRNG,
+    ϱ1::GibbsPolarSliceTarget,
+    ℓT::F,
+    θ::AbstractVector{F},
+    r::F,
+    max_prop::Int,
+) where {F<:Real}
+    y     = rand_subsphere(rng, θ)
+    ω_max = convert(F, 2π) * rand(rng, F)
     ω_min = ω_max - convert(F, 2π)
 
     for n_props in 1:max_prop
         # `Uniform` had a type instability issue:
         # https://github.com/JuliaStats/Distributions.jl/pull/1860
         # ω = rand(rng, Uniform(ω_min, ω_max))
-        ω = ω_min + (ω_max - ω_min)*rand(rng, F)
-        θ′ = θ*cos(ω) + y*sin(ω)
+        ω = ω_min + (ω_max - ω_min) * rand(rng, F)
+        θ′ = θ * cos(ω) + y * sin(ω)
 
-        if logdensity(ϱ1, r*θ′) > ℓT
+        if logdensity(ϱ1, r * θ′) > ℓT
             return θ′, n_props
         end
 
@@ -105,27 +109,27 @@ function geodesic_shrinkage(
             ω_max = ω
         end
     end
-    exceeded_max_prop(max_prop)
+    return exceeded_max_prop(max_prop)
 end
 
 function radius_shrinkage(
-    rng     ::Random.AbstractRNG,
-    ϱ1      ::GibbsPolarSliceTarget,
-    ℓT      ::F,
-    θ       ::AbstractVector{F},
-    r       ::F,
-    w       ::Real,
-    max_prop::Int
-) where {F <: Real}
+    rng::Random.AbstractRNG,
+    ϱ1::GibbsPolarSliceTarget,
+    ℓT::F,
+    θ::AbstractVector{F},
+    r::F,
+    w::Real,
+    max_prop::Int,
+) where {F<:Real}
     u             = rand(rng, F)
     w             = convert(F, w)
-    r_min         = max(r - u*w, 0)
-    r_max         = r + (1 - u)*w
+    r_min         = max(r - u * w, 0)
+    r_max         = r + (1 - u) * w
     n_props_total = 0
 
     n_props = 0
-    while (r_min > 0) && logdensity(ϱ1, r_min*θ) > ℓT
-        r_min    = max(r_min - w, 0)
+    while (r_min > 0) && logdensity(ϱ1, r_min * θ) > ℓT
+        r_min = max(r_min - w, 0)
 
         n_props += 1
         if n_props > max_prop
@@ -135,7 +139,7 @@ function radius_shrinkage(
     n_props_total += n_props
 
     n_props = 0
-    while logdensity(ϱ1, r_max*θ) > ℓT
+    while logdensity(ϱ1, r_max * θ) > ℓT
         r_max = r_max + w
 
         n_props += 1
@@ -149,27 +153,27 @@ function radius_shrinkage(
         # `Uniform` had a type instability issue:
         # https://github.com/JuliaStats/Distributions.jl/pull/1860
         #r′ = rand(rng, Uniform{F}(r_min, r_max))
-        r′ = r_min + (r_max - r_min)*rand(rng, F)
+        r′ = r_min + (r_max - r_min) * rand(rng, F)
 
-        if logdensity(ϱ1, r′*θ) > ℓT
+        if logdensity(ϱ1, r′ * θ) > ℓT
             n_props_total += n_props
             return r′, n_props_total
         end
 
-        if r′ < r      
+        if r′ < r
             r_min = r′
         else
             r_max = r′
         end
     end
-    exceeded_max_prop(max_prop)
+    return exceeded_max_prop(max_prop)
 end
 
 function AbstractMCMC.step(
-    rng    ::Random.AbstractRNG,
-    model  ::AbstractMCMC.LogDensityModel, 
+    rng::Random.AbstractRNG,
+    model::AbstractMCMC.LogDensityModel,
     sampler::GibbsPolarSlice,
-    state  ::GibbsPolarSliceState;
+    state::GibbsPolarSliceState;
     kwargs...,
 )
     logdensitymodel = model.logdensity
@@ -183,16 +187,13 @@ function AbstractMCMC.step(
     ϱ1 = GibbsPolarSliceTarget(logdensitymodel)
 
     d  = length(x)
-    ℓT = ((d-1)*log(norm(x)) + ℓp) - Random.randexp(rng, eltype(ℓp))
+    ℓT = ((d - 1) * log(norm(x)) + ℓp) - Random.randexp(rng, eltype(ℓp))
 
-    θ, n_props_θ = geodesic_shrinkage(rng, ϱ1, ℓT, θ, r,    max_prop)
-    r, n_props_r = radius_shrinkage(  rng, ϱ1, ℓT, θ, r, w, max_prop)
-    x            = θ*r
+    θ, n_props_θ = geodesic_shrinkage(rng, ϱ1, ℓT, θ, r, max_prop)
+    r, n_props_r = radius_shrinkage(rng, ϱ1, ℓT, θ, r, w, max_prop)
+    x            = θ * r
 
     ℓp = LogDensityProblems.logdensity(logdensitymodel, x)
-    t  = Transition(x, ℓp, (
-        num_radius_proposals    = n_props_r,
-        num_direction_proposals = n_props_θ,
-    ))
-    t, GibbsPolarSliceState(t, θ, r)
+    t  = Transition(x, ℓp, (num_radius_proposals=n_props_r, num_direction_proposals=n_props_θ))
+    return t, GibbsPolarSliceState(t, θ, r)
 end
diff --git a/src/multivariate/hitandrun.jl b/src/multivariate/hitandrun.jl
index 4fa4b7d..f05ce72 100644
--- a/src/multivariate/hitandrun.jl
+++ b/src/multivariate/hitandrun.jl
@@ -8,36 +8,36 @@ This applies `unislice` along a random direction uniform sampled from the sphere
 # Arguments
 - `unislice::AbstractUnivariateSliceSampling`: Univariate slice sampling algorithm.
 """
-struct HitAndRun{
-    S <: AbstractUnivariateSliceSampling
-} <: AbstractMultivariateSliceSampling
+struct HitAndRun{S<:AbstractUnivariateSliceSampling} <: AbstractMultivariateSliceSampling
     unislice::S
 end
 
-struct HitAndRunState{T <: Transition}
+struct HitAndRunState{T<:Transition}
     "Current [`Transition`](@ref)."
     transition::T
 end
 
-struct HitAndRunTarget{Model, Vec <: AbstractVector}
-    model    ::Model
-    direction::Vec
-    reference::Vec
+struct HitAndRunTarget{Model,Vec<:AbstractVector}
+    model     :: Model
+    direction :: Vec
+    reference :: Vec
 end
 
 function LogDensityProblems.logdensity(target::HitAndRunTarget, λ)
     (; model, reference, direction) = target
-    LogDensityProblems.logdensity(model, reference + λ*direction)
+    return LogDensityProblems.logdensity(model, reference + λ * direction)
 end
 
-function AbstractMCMC.step(rng    ::Random.AbstractRNG,
-                           model  ::AbstractMCMC.LogDensityModel,
-                           sampler::HitAndRun;
-                           initial_params = nothing,
-                           kwargs...)
+function AbstractMCMC.step(
+    rng::Random.AbstractRNG,
+    model::AbstractMCMC.LogDensityModel,
+    sampler::HitAndRun;
+    initial_params=nothing,
+    kwargs...,
+)
     logdensitymodel = model.logdensity
-    θ  = isnothing(initial_params) ? initial_sample(rng, logdensitymodel) : initial_params
-    d  = length(θ)
+    θ = isnothing(initial_params) ? initial_sample(rng, logdensitymodel) : initial_params
+    d = length(θ)
     @assert d ≥ 2 "Hit-and-Run works reliably only in dimension ≥2"
     lp = LogDensityProblems.logdensity(logdensitymodel, θ)
     t  = Transition(θ, lp, NamedTuple())
@@ -46,14 +46,14 @@ end
 
 function rand_uniform_unit_sphere(rng::Random.AbstractRNG, type::Type, d::Int)
     x = randn(rng, type, d)
-    x / norm(x)
+    return x / norm(x)
 end
 
 function AbstractMCMC.step(
-    rng    ::Random.AbstractRNG,
-    model  ::AbstractMCMC.LogDensityModel, 
+    rng::Random.AbstractRNG,
+    model::AbstractMCMC.LogDensityModel,
     sampler::HitAndRun,
-    state  ::HitAndRunState;
+    state::HitAndRunState;
     kwargs...,
 )
     logdensitymodel = model.logdensity
@@ -65,10 +65,8 @@ function AbstractMCMC.step(
     direction    = rand_uniform_unit_sphere(rng, eltype(θ), d)
     hnrtarget    = HitAndRunTarget(logdensitymodel, direction, θ)
     λ            = zero(eltype(θ))
-    λ, ℓp, props = slice_sampling_univariate(
-        rng, unislice, hnrtarget, ℓp, λ
-    )
-    θ′ = θ + direction*λ
-    t = Transition(θ′, ℓp, (num_proposals=props,))
-    t, HitAndRunState(t)
+    λ, ℓp, props = slice_sampling_univariate(rng, unislice, hnrtarget, ℓp, λ)
+    θ′           = θ + direction * λ
+    t            = Transition(θ′, ℓp, (num_proposals=props,))
+    return t, HitAndRunState(t)
 end
diff --git a/src/multivariate/latent.jl b/src/multivariate/latent.jl
index d3f7051..9f0e288 100644
--- a/src/multivariate/latent.jl
+++ b/src/multivariate/latent.jl
@@ -10,48 +10,50 @@ Latent slice sampling algorithm by Li and Walker[^LW2023].
 # Keyword Arguments
 - `max_proposals::Int`: Maximum number of proposals allowed until throwing an error (default: `$(DEFAULT_MAX_PROPOSALS)`).
 """
-struct LatentSlice{B <: Real} <: AbstractMultivariateSliceSampling
-    beta         ::B
-    max_proposals::Int
+struct LatentSlice{B<:Real} <: AbstractMultivariateSliceSampling
+    beta          :: B
+    max_proposals :: Int
 end
 
-function LatentSlice(beta::Real; max_proposals::Int = DEFAULT_MAX_PROPOSALS)
+function LatentSlice(beta::Real; max_proposals::Int=DEFAULT_MAX_PROPOSALS)
     @assert beta > 0 "Beta must be strictly positive"
-    LatentSlice(beta, max_proposals)
+    return LatentSlice(beta, max_proposals)
 end
 
-struct LatentSliceState{T <: Transition, S <: AbstractVector}
+struct LatentSliceState{T<:Transition,S<:AbstractVector}
     "Current [`Transition`](@ref)."
-    transition ::T
+    transition::T
 
     "Auxiliary variables for adapting the slice window (\$s\$ in the original paper[^LW2023])"
     sliceparams::S
 end
 
-function AbstractMCMC.step(rng    ::Random.AbstractRNG,
-                           model  ::AbstractMCMC.LogDensityModel,
-                           sampler::LatentSlice;
-                           initial_params = nothing,
-                           kwargs...)
+function AbstractMCMC.step(
+    rng::Random.AbstractRNG,
+    model::AbstractMCMC.LogDensityModel,
+    sampler::LatentSlice;
+    initial_params=nothing,
+    kwargs...,
+)
     logdensitymodel = model.logdensity
-    y  = initial_params === nothing ? initial_sample(rng, logdensitymodel) : initial_params
-    β  = sampler.beta
-    d  = length(y)
+    y = initial_params === nothing ? initial_sample(rng, logdensitymodel) : initial_params
+    β = sampler.beta
+    d = length(y)
     lp = LogDensityProblems.logdensity(logdensitymodel, y)
-    s  = convert(Vector{eltype(y)}, rand(rng, Gamma(2, 1/β), d))
-    t  = Transition(y, lp, NamedTuple())
+    s = convert(Vector{eltype(y)}, rand(rng, Gamma(2, 1 / β), d))
+    t = Transition(y, lp, NamedTuple())
     return t, LatentSliceState(t, s)
 end
 
 function AbstractMCMC.step(
-    rng    ::Random.AbstractRNG,
-    model  ::AbstractMCMC.LogDensityModel, 
+    rng::Random.AbstractRNG,
+    model::AbstractMCMC.LogDensityModel,
     sampler::LatentSlice,
-    state  ::LatentSliceState;
+    state::LatentSliceState;
     kwargs...,
 )
     logdensitymodel = model.logdensity
-    max_proposals   = sampler.max_proposals 
+    max_proposals   = sampler.max_proposals
 
     β  = sampler.beta
     ℓp = state.transition.lp
@@ -61,16 +63,16 @@ function AbstractMCMC.step(
     ℓw = ℓp - Random.randexp(rng, eltype(y))
 
     u_l = rand(rng, eltype(y), d)
-    l   = (y - s/2) + u_l.*s
-    a   = l - s/2
-    b   = l + s/2
+    l   = (y - s / 2) + u_l .* s
+    a   = l - s / 2
+    b   = l + s / 2
 
     props = 0
     while true
         props += 1
 
         u_y    = rand(rng, eltype(y), d)
-        ystar  = a + u_y.*(b - a)
+        ystar  = a + u_y .* (b - a)
         ℓpstar = LogDensityProblems.logdensity(logdensitymodel, ystar)
 
         if ℓw < ℓpstar
@@ -83,7 +85,7 @@ function AbstractMCMC.step(
             exceeded_max_prop(max_proposals)
         end
 
-        @inbounds for i = 1:d
+        @inbounds for i in 1:d
             if ystar[i] < y[i]
                 a[i] = ystar[i]
             else
@@ -91,7 +93,7 @@ function AbstractMCMC.step(
             end
         end
     end
-    s = β*randexp(rng, eltype(y), d) + 2*abs.(l - y)
-    t = Transition(y, ℓp, (num_proposals = props,))
-    t, LatentSliceState(t, s)
+    s = β * randexp(rng, eltype(y), d) + 2 * abs.(l - y)
+    t = Transition(y, ℓp, (num_proposals=props,))
+    return t, LatentSliceState(t, s)
 end
diff --git a/src/multivariate/randpermgibbs.jl b/src/multivariate/randpermgibbs.jl
index c5b6372..366c320 100644
--- a/src/multivariate/randpermgibbs.jl
+++ b/src/multivariate/randpermgibbs.jl
@@ -12,38 +12,40 @@ When `unislice` is a vector of samplers, each slice sampler is applied to the co
 In that case, the `length(unislice)` must match the dimensionality of the posterior.
 """
 struct RandPermGibbs{
-    S <: Union{
-        <: AbstractUnivariateSliceSampling,
-        <: AbstractVector{<: AbstractUnivariateSliceSampling}
-    }
+    S<:Union{
+        <:AbstractUnivariateSliceSampling,
+        <:AbstractVector{<:AbstractUnivariateSliceSampling},
+    },
 } <: AbstractMultivariateSliceSampling
     unislice::S
 end
 
-struct GibbsState{T <: Transition}
+struct GibbsState{T<:Transition}
     "Current [`Transition`](@ref)."
     transition::T
 end
 
-struct GibbsTarget{Model, Idx <: Integer, Vec <: AbstractVector}
-    model::Model
-    idx  ::Idx
-    θ    ::Vec
+struct GibbsTarget{Model,Idx<:Integer,Vec<:AbstractVector}
+    model :: Model
+    idx   :: Idx
+    θ     :: Vec
 end
 
 function LogDensityProblems.logdensity(gibbs::GibbsTarget, θi)
     (; model, idx, θ) = gibbs
-    LogDensityProblems.logdensity(model, (@set θ[idx] = θi))
+    return LogDensityProblems.logdensity(model, (@set θ[idx] = θi))
 end
 
-function AbstractMCMC.step(rng    ::Random.AbstractRNG,
-                           model  ::AbstractMCMC.LogDensityModel,
-                           sampler::RandPermGibbs;
-                           initial_params = nothing,
-                           kwargs...)
+function AbstractMCMC.step(
+    rng::Random.AbstractRNG,
+    model::AbstractMCMC.LogDensityModel,
+    sampler::RandPermGibbs;
+    initial_params=nothing,
+    kwargs...,
+)
     logdensitymodel = model.logdensity
-    θ  = initial_params === nothing ? initial_sample(rng, logdensitymodel) : initial_params
-    d  = length(θ)
+    θ = initial_params === nothing ? initial_sample(rng, logdensitymodel) : initial_params
+    d = length(θ)
     if sampler.unislice isa AbstractVector
         @assert length(sampler.unislice) == d "Number of slice samplers does not match the dimensionality of the initial parameter."
     end
@@ -53,10 +55,10 @@ function AbstractMCMC.step(rng    ::Random.AbstractRNG,
 end
 
 function AbstractMCMC.step(
-    rng    ::Random.AbstractRNG,
-    model  ::AbstractMCMC.LogDensityModel, 
+    rng::Random.AbstractRNG,
+    model::AbstractMCMC.LogDensityModel,
     sampler::RandPermGibbs,
-    state  ::GibbsState;
+    state::GibbsState;
     kwargs...,
 )
     logdensitymodel = model.logdensity
@@ -72,13 +74,13 @@ function AbstractMCMC.step(
     props = zeros(Int, d)
     for i in shuffle(rng, 1:d)
         model_gibbs = GibbsTarget(logdensitymodel, i, θ)
-        unislice    = unislices[i]
+        unislice = unislices[i]
         θ′_coord, ℓp, props_coord = slice_sampling_univariate(
             rng, unislice, model_gibbs, ℓp, θ[i]
         )
         props[i] = props_coord
-        θ[i]     = θ′_coord
+        θ[i] = θ′_coord
     end
     t = Transition(θ, ℓp, (num_proposals=props,))
-    t, GibbsState(t)
+    return t, GibbsState(t)
 end
diff --git a/src/univariate/doublingout.jl b/src/univariate/doublingout.jl
index 3a9155e..17a77b0 100644
--- a/src/univariate/doublingout.jl
+++ b/src/univariate/doublingout.jl
@@ -11,40 +11,33 @@ Univariate slice sampling by automatically adapting the initial interval through
 - `max_doubling_out`: Maximum number of "doubling outs" (default: 8).
 - `max_proposals::Int`: Maximum number of proposals allowed until throwing an error (default: `$(DEFAULT_MAX_PROPOSALS)`).
 """
-struct SliceDoublingOut{W <: Real} <: AbstractUnivariateSliceSampling
-    window          ::W
-    max_doubling_out::Int
-    max_proposals   ::Int
+struct SliceDoublingOut{W<:Real} <: AbstractUnivariateSliceSampling
+    window           :: W
+    max_doubling_out :: Int
+    max_proposals    :: Int
 end
 
 function SliceDoublingOut(
-    window          ::Real;
-    max_doubling_out::Int = 8,
-    max_proposals   ::Int = DEFAULT_MAX_PROPOSALS,
+    window::Real; max_doubling_out::Int=8, max_proposals::Int=DEFAULT_MAX_PROPOSALS
 )
     @assert window > 0
-    SliceDoublingOut(window, max_doubling_out, max_proposals)
+    return SliceDoublingOut(window, max_doubling_out, max_proposals)
 end
 
 function find_interval(
-    rng  ::Random.AbstractRNG,
-    alg  ::SliceDoublingOut,
-    model,
-    w    ::Real,
-    ℓy   ::Real,
-    θ₀   ::F,
-) where {F <: Real}
+    rng::Random.AbstractRNG, alg::SliceDoublingOut, model, w::Real, ℓy::Real, θ₀::F
+) where {F<:Real}
     p = alg.max_doubling_out
 
     u = rand(rng, F)
-    L = θ₀ - w*u
+    L = θ₀ - w * u
     R = L + w
 
     ℓπ_L = LogDensityProblems.logdensity(model, L)
     ℓπ_R = LogDensityProblems.logdensity(model, R)
     K    = 2
 
-    for _ = 1:p
+    for _ in 1:p
         if ((ℓy ≥ ℓπ_L) && (ℓy ≥ ℓπ_R))
             break
         end
@@ -58,25 +51,18 @@ function find_interval(
         end
         K += 1
     end
-    L, R, K
+    return L, R, K
 end
 
 function accept_slice_proposal(
-         ::SliceDoublingOut,
-    model,
-    w    ::Real,
-    ℓy   ::Real,
-    θ₀   ::Real,
-    θ₁   ::Real,
-    L    ::Real,
-    R    ::Real,
-) 
+    ::SliceDoublingOut, model, w::Real, ℓy::Real, θ₀::Real, θ₁::Real, L::Real, R::Real
+)
     D    = false
     ℓπ_L = LogDensityProblems.logdensity(model, L)
     ℓπ_R = LogDensityProblems.logdensity(model, R)
 
-    while R - L > 1.1*w
-        M = (L + R)/2
+    while R - L > 1.1 * w
+        M = (L + R) / 2
         if (θ₀ < M && θ₁ ≥ M) || (θ₀ ≥ M && θ₁ < M)
             D = true
         end
@@ -93,5 +79,5 @@ function accept_slice_proposal(
             return false
         end
     end
-    true
+    return true
 end
diff --git a/src/univariate/fixedinterval.jl b/src/univariate/fixedinterval.jl
index b679ec3..df61f20 100644
--- a/src/univariate/fixedinterval.jl
+++ b/src/univariate/fixedinterval.jl
@@ -10,30 +10,21 @@ Univariate slice sampling with a fixed initial interval (Scheme 2 by Neal[^N2003
 # Keyword Arguments
 - `max_proposals::Int`: Maximum number of proposals allowed until throwing an error (default: `$(DEFAULT_MAX_PROPOSALS)`).
 """
-struct Slice{W <: Real} <: AbstractUnivariateSliceSampling
-    window       ::W
-    max_proposals::Int
+struct Slice{W<:Real} <: AbstractUnivariateSliceSampling
+    window        :: W
+    max_proposals :: Int
 end
 
-function Slice(
-    window       ::Real;
-    max_proposals::Int = DEFAULT_MAX_PROPOSALS, 
-)
+function Slice(window::Real; max_proposals::Int=DEFAULT_MAX_PROPOSALS)
     @assert window > 0
-    Slice(window, max_proposals)
+    return Slice(window, max_proposals)
 end
 
 function find_interval(
-    rng::Random.AbstractRNG,
-       ::Slice,
-       ::Any,
-    w  ::Real,
-       ::Real,
-    θ₀ ::F,
-) where {F <: Real}
+    rng::Random.AbstractRNG, ::Slice, ::Any, w::Real, ::Real, θ₀::F
+) where {F<:Real}
     u = rand(rng, F)
-    L = θ₀ - w*u
+    L = θ₀ - w * u
     R = L + w
-    L, R, 0
+    return L, R, 0
 end
-
diff --git a/src/univariate/steppingout.jl b/src/univariate/steppingout.jl
index 88d578e..8a81fcc 100644
--- a/src/univariate/steppingout.jl
+++ b/src/univariate/steppingout.jl
@@ -11,36 +11,29 @@ Univariate slice sampling by automatically adapting the initial interval through
 - `max_stepping_out::Int`: Maximum number of "stepping outs" (default: 32).
 - `max_proposals::Int`: Maximum number of proposals allowed until throwing an error (default: `$(DEFAULT_MAX_PROPOSALS)`).
 """
-struct SliceSteppingOut{W <: Real} <: AbstractUnivariateSliceSampling
-    window          ::W
-    max_stepping_out::Int
-    max_proposals   ::Int
+struct SliceSteppingOut{W<:Real} <: AbstractUnivariateSliceSampling
+    window           :: W
+    max_stepping_out :: Int
+    max_proposals    :: Int
 end
 
 function SliceSteppingOut(
-    window          ::Real;
-    max_stepping_out::Int = 32,
-    max_proposals   ::Int = DEFAULT_MAX_PROPOSALS,
+    window::Real; max_stepping_out::Int=32, max_proposals::Int=DEFAULT_MAX_PROPOSALS
 )
     @assert window > 0
-    SliceSteppingOut(window, max_stepping_out, max_proposals)
+    return SliceSteppingOut(window, max_stepping_out, max_proposals)
 end
 
 function find_interval(
-    rng  ::Random.AbstractRNG,
-    alg  ::SliceSteppingOut,
-    model,
-    w    ::Real,
-    ℓy   ::Real,
-    θ₀   ::F,
-) where {F <: Real}
+    rng::Random.AbstractRNG, alg::SliceSteppingOut, model, w::Real, ℓy::Real, θ₀::F
+) where {F<:Real}
     m      = alg.max_stepping_out
     u      = rand(rng, F)
-    L      = θ₀ - w*u
+    L      = θ₀ - w * u
     R      = L + w
     V      = rand(rng, F)
-    J      = floor(Int, m*V)
-    K      = (m - 1) - J 
+    J      = floor(Int, m * V)
+    K      = (m - 1) - J
     n_eval = 0
 
     while J > 0 && ℓy < LogDensityProblems.logdensity(model, L)
@@ -53,5 +46,5 @@ function find_interval(
         K = K - 1
         n_eval += 1
     end
-    L, R, n_eval
+    return L, R, n_eval
 end
diff --git a/src/univariate/univariate.jl b/src/univariate/univariate.jl
index d2719d8..0498b17 100644
--- a/src/univariate/univariate.jl
+++ b/src/univariate/univariate.jl
@@ -1,19 +1,15 @@
 
 function slice_sampling_univariate(
-    rng   ::Random.AbstractRNG,
-    alg   ::AbstractSliceSampling,
-    model,
-    ℓπ    ::Real,
-    θ     ::F,
-) where {F <: Real}
+    rng::Random.AbstractRNG, alg::AbstractSliceSampling, model, ℓπ::Real, θ::F
+) where {F<:Real}
     w, max_prop = alg.window, alg.max_proposals
     ℓy          = ℓπ - Random.randexp(rng, F)
     L, R, props = find_interval(rng, alg, model, w, ℓy, θ)
 
     for _ in 1:max_prop
         U     = rand(rng, F)
-        θ′     = L + U*(R - L)
-        ℓπ′    = LogDensityProblems.logdensity(model, θ′)
+        θ′    = L + U * (R - L)
+        ℓπ′   = LogDensityProblems.logdensity(model, θ′)
         props += 1
         if (ℓy < ℓπ′) && accept_slice_proposal(alg, model, w, ℓy, θ, θ′, L, R)
             return θ′, ℓπ′, props
@@ -25,21 +21,23 @@ function slice_sampling_univariate(
             R = θ′
         end
     end
-    exceeded_max_prop(max_prop)
+    return exceeded_max_prop(max_prop)
 end
 
-struct UnivariateSliceState{T <: Transition}
+struct UnivariateSliceState{T<:Transition}
     "Current [`Transition`](@ref)."
     transition::T
 end
 
-function AbstractMCMC.step(rng    ::Random.AbstractRNG,
-                           model  ::AbstractMCMC.LogDensityModel,
-                           sampler::AbstractUnivariateSliceSampling;
-                           initial_params = nothing,
-                           kwargs...)
+function AbstractMCMC.step(
+    rng::Random.AbstractRNG,
+    model::AbstractMCMC.LogDensityModel,
+    sampler::AbstractUnivariateSliceSampling;
+    initial_params=nothing,
+    kwargs...,
+)
     logdensitymodel = model.logdensity
-    θ  = isnothing(initial_params) ? initial_sample(rng, logdensitymodel) : initial_params
+    θ = isnothing(initial_params) ? initial_sample(rng, logdensitymodel) : initial_params
     @assert length(θ) == 1 "The dimensionality of the parameter should be 1."
     lp = LogDensityProblems.logdensity(logdensitymodel, θ)
     t  = Transition(θ, lp, NamedTuple())
@@ -51,14 +49,14 @@ struct UnivariateTarget{Model}
 end
 
 function LogDensityProblems.logdensity(uni::UnivariateTarget, θi)
-    LogDensityProblems.logdensity(uni.model, [θi])
+    return LogDensityProblems.logdensity(uni.model, [θi])
 end
 
 function AbstractMCMC.step(
-    rng    ::Random.AbstractRNG,
-    model  ::AbstractMCMC.LogDensityModel, 
+    rng::Random.AbstractRNG,
+    model::AbstractMCMC.LogDensityModel,
     sampler::AbstractUnivariateSliceSampling,
-    state  ::UnivariateSliceState;
+    state::UnivariateSliceState;
     kwargs...,
 )
     logdensitymodel = model.logdensity
@@ -69,5 +67,5 @@ function AbstractMCMC.step(
     )
 
     t = Transition([θ], ℓp, (num_proposals=props,))
-    t, UnivariateSliceState(t)
+    return t, UnivariateSliceState(t)
 end
diff --git a/test/maxprops.jl b/test/maxprops.jl
index d52911f..aa7711f 100644
--- a/test/maxprops.jl
+++ b/test/maxprops.jl
@@ -4,11 +4,11 @@ struct WrongModel end
 LogDensityProblems.logdensity(::WrongModel, θ) = -Inf
 
 function LogDensityProblems.capabilities(::Type{<:WrongModel})
-    LogDensityProblems.LogDensityOrder{0}()
+    return LogDensityProblems.LogDensityOrder{0}()
 end
 
 function LogDensityProblems.dimension(::WrongModel)
-    2
+    return 2
 end
 
 @testset "error handling" begin
@@ -18,7 +18,6 @@ end
         RandPermGibbs(Slice(1; max_proposals=32)),
         RandPermGibbs(SliceSteppingOut(1; max_proposals=32)),
         RandPermGibbs(SliceDoublingOut(1; max_proposals=32)),
-
         HitAndRun(Slice(1; max_proposals=32)),
         HitAndRun(SliceSteppingOut(1; max_proposals=32)),
         HitAndRun(SliceDoublingOut(1; max_proposals=32)),
@@ -31,7 +30,7 @@ end
     ]
         @testset "max proposal error" begin
             rng = Random.default_rng()
-            θ   = [1., 1.]
+            θ = [1.0, 1.0]
             _, init_state = AbstractMCMC.step(rng, model, sampler; initial_params=copy(θ))
 
             @test_throws ErrorException begin
diff --git a/test/multivariate.jl b/test/multivariate.jl
index 910fe63..5e081ae 100644
--- a/test/multivariate.jl
+++ b/test/multivariate.jl
@@ -1,64 +1,65 @@
 
-struct MultiModel{F <: Real, V <: AbstractVector}
+struct MultiModel{F<:Real,V<:AbstractVector}
     α::F
     β::F
     y::V
 end
 
 function MCMCTesting.sample_joint(
-    rng::AbstractRNG, model::MultiModel{F, V}
-) where {F <: Real, V <: AbstractVector}
+    rng::AbstractRNG, model::MultiModel{F,V}
+) where {F<:Real,V<:AbstractVector}
     α, β = model.α, model.β
 
     μ = rand(rng, Normal(zero(F), one(F)))
-    σ = rand(rng, InverseGamma(α, β)) |> F # InverseGamma is not type stable
+    σ = F(rand(rng, InverseGamma(α, β))) # InverseGamma is not type stable
     y = rand(rng, Normal(μ, σ), 10)
 
     θ = [μ, σ]
-    θ, y
+    return θ, y
 end
 
 function MCMCTesting.markovchain_transition(
-    rng    ::Random.AbstractRNG,
-    model  ::MultiModel,
+    rng::Random.AbstractRNG,
+    model::MultiModel,
     sampler::SliceSampling.AbstractSliceSampling,
-    θ, y
+    θ,
+    y,
 )
     model′ = AbstractMCMC.LogDensityModel(@set model.y = y)
     _, init_state = AbstractMCMC.step(rng, model′, sampler; initial_params=copy(θ))
     transition, _ = AbstractMCMC.step(rng, model′, sampler, init_state)
-    transition.params
+    return transition.params
 end
 
-function LogDensityProblems.logdensity(model::MultiModel{F, V}, θ) where {F <: Real, V}
+function LogDensityProblems.logdensity(model::MultiModel{F,V}, θ) where {F<:Real,V}
     α, β, y = model.α, model.β, model.y
 
     μ = θ[1]
     σ = θ[2]
 
     if σ ≤ 0
-       return typemin(F)
+        return typemin(F)
     end
 
-    logpdf(Normal(zero(F), one(F)), μ) +
-        logpdf(InverseGamma(α, β), σ) +
-        sum(Base.Fix1(logpdf, Normal(μ, σ)), y)
+    return logpdf(Normal(zero(F), one(F)), μ) +
+           logpdf(InverseGamma(α, β), σ) +
+           sum(Base.Fix1(logpdf, Normal(μ, σ)), y)
 end
 
 function SliceSampling.initial_sample(rng::Random.AbstractRNG, model::MultiModel)
-    randn(rng, LogDensityProblems.dimension(model))
+    return randn(rng, LogDensityProblems.dimension(model))
 end
 
 function LogDensityProblems.capabilities(::Type{<:MultiModel})
-    LogDensityProblems.LogDensityOrder{0}()
+    return LogDensityProblems.LogDensityOrder{0}()
 end
 
 function LogDensityProblems.dimension(model::MultiModel)
-    2
+    return 2
 end
 
 @testset "multivariate samplers" begin
-    model   = MultiModel(1., 1., [0.])
+    model = MultiModel(1.0, 1.0, [0.0])
     @testset for sampler in [
         # Vector-valued windows
         RandPermGibbs(Slice.(fill(1, LogDensityProblems.dimension(model)))),
@@ -69,7 +70,6 @@ end
         RandPermGibbs(Slice(1)),
         RandPermGibbs(SliceSteppingOut(1)),
         RandPermGibbs(SliceDoublingOut(1)),
-
         HitAndRun(Slice(1)),
         HitAndRun(SliceSteppingOut(1)),
         HitAndRun(SliceDoublingOut(1)),
@@ -82,17 +82,11 @@ end
     ]
         @testset "initial_params" begin
             model  = MultiModel(1.0, 1.0, [0.0])
-            θ, y  = MCMCTesting.sample_joint(Random.default_rng(), model)
+            θ, y   = MCMCTesting.sample_joint(Random.default_rng(), model)
             model′ = AbstractMCMC.LogDensityModel(@set model.y = y)
 
             θ0    = [1.0, 0.1]
-            chain = sample(
-                model′,
-                sampler,
-                10;
-                initial_params=θ0,
-                progress=false,
-            )
+            chain = sample(model′, sampler, 10; initial_params=θ0, progress=false)
             @test first(chain).params == θ0
         end
 
@@ -108,7 +102,7 @@ end
 
         @testset "determinism" begin
             model  = MultiModel(1.0, 1.0, [0.0])
-            θ, y  = MCMCTesting.sample_joint(Random.default_rng(), model)
+            θ, y   = MCMCTesting.sample_joint(Random.default_rng(), model)
             model′ = AbstractMCMC.LogDensityModel(@set model.y = y)
 
             rng           = StableRNG(1)
@@ -119,14 +113,14 @@ end
             rng           = StableRNG(1)
             _, init_state = AbstractMCMC.step(rng, model′, sampler; initial_params=copy(θ))
             transition, _ = AbstractMCMC.step(rng, model′, sampler, init_state)
-            θ′′            = transition.params
+            θ′′           = transition.params
             @test θ′ == θ′′
         end
 
         @testset "type stability $(type)" for type in [Float32, Float64]
-            rng   = Random.default_rng()
+            rng = Random.default_rng()
             model = MultiModel(one(type), one(type), [zero(type)])
-            θ, y  = MCMCTesting.sample_joint(Random.default_rng(), model)
+            θ, y = MCMCTesting.sample_joint(Random.default_rng(), model)
             model′ = AbstractMCMC.LogDensityModel(@set model.y = y)
 
             @test eltype(θ) == type
@@ -134,7 +128,7 @@ end
 
             _, init_state = AbstractMCMC.step(rng, model′, sampler; initial_params=copy(θ))
             transition, _ = AbstractMCMC.step(rng, model′, sampler, init_state)
-            θ′             = transition.params
+            θ′ = transition.params
 
             @test eltype(θ′) == type
         end
@@ -146,7 +140,7 @@ end
             n_mcmc_thin      = 10
             test             = ExactRankTest(n_samples, n_mcmc_steps, n_mcmc_thin)
 
-            model   = MultiModel(1., 1., [0.])
+            model   = MultiModel(1.0, 1.0, [0.0])
             subject = TestSubject(model, sampler)
             @test seqmcmctest(test, subject, 0.001, n_pvalue_samples; show_progress=false)
         end
diff --git a/test/runtests.jl b/test/runtests.jl
index ac126a8..23b4650 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,7 +4,7 @@ using Accessors
 using Distributions
 using LogDensityProblems
 using MCMCTesting
-using Random 
+using Random
 using Test
 using Turing
 using StableRNGs
diff --git a/test/turing.jl b/test/turing.jl
index 6df11e5..725828a 100644
--- a/test/turing.jl
+++ b/test/turing.jl
@@ -4,7 +4,7 @@
         s ~ InverseGamma(2, 3)
         m ~ Normal(0, sqrt(s))
         1.5 ~ Normal(m, sqrt(s))
-        2.0 ~ Normal(m, sqrt(s))
+        return 2.0 ~ Normal(m, sqrt(s))
     end
 
     n_samples = 1000
@@ -28,12 +28,7 @@
             progress=false,
         )
 
-        chain = sample(
-            model,
-            externalsampler(sampler),
-            n_samples;
-            progress=false,
-        )
+        chain = sample(model, externalsampler(sampler), n_samples; progress=false)
     end
 
     @testset "gibbs($sampler)" for sampler in [
@@ -46,13 +41,10 @@
     ]
         sample(
             model,
-            Turing.Experimental.Gibbs(
-                (
-                    s = externalsampler(sampler),
-                    m = externalsampler(sampler),
-                ),
-            ),
-            n_samples,
+            Turing.Experimental.Gibbs((
+                s=externalsampler(sampler), m=externalsampler(sampler)
+            ),),
+            n_samples;
             progress=false,
         )
     end
diff --git a/test/univariate.jl b/test/univariate.jl
index a18e783..55b9151 100644
--- a/test/univariate.jl
+++ b/test/univariate.jl
@@ -1,70 +1,60 @@
 
-
-struct UniModel{V <: AbstractVector}
+struct UniModel{V<:AbstractVector}
     y::V
 end
 
 function MCMCTesting.sample_joint(
     rng::AbstractRNG, ::UniModel{<:AbstractVector{F}}
-) where {F <: Real}
+) where {F<:Real}
     μ = rand(rng, Normal(zero(F), one(F)))
     y = rand(rng, Normal(μ, one(F)), 10)
-    [μ], y
+    return [μ], y
 end
 
 function MCMCTesting.markovchain_transition(
-    rng    ::Random.AbstractRNG,
-    model  ::UniModel,
+    rng::Random.AbstractRNG,
+    model::UniModel,
     sampler::SliceSampling.AbstractSliceSampling,
-    θ, y
+    θ,
+    y,
 )
     model′ = AbstractMCMC.LogDensityModel(@set model.y = y)
     _, init_state = AbstractMCMC.step(rng, model′, sampler; initial_params=copy(θ))
     transition, _ = AbstractMCMC.step(rng, model′, sampler, init_state)
-    transition.params
+    return transition.params
 end
 
 function LogDensityProblems.logdensity(
     model::UniModel{<:AbstractVector{F}}, θ
-) where {F <: Real}
+) where {F<:Real}
     y = model.y
     μ = only(θ)
 
-    logpdf(Normal(zero(F), one(F)), μ) + sum(Base.Fix1(logpdf, Normal(μ, one(F))), y)
+    return logpdf(Normal(zero(F), one(F)), μ) + sum(Base.Fix1(logpdf, Normal(μ, one(F))), y)
 end
 
 function SliceSampling.initial_sample(rng::Random.AbstractRNG, model::UniModel)
-    randn(rng, LogDensityProblems.dimension(model))
+    return randn(rng, LogDensityProblems.dimension(model))
 end
 
 function LogDensityProblems.capabilities(::Type{<:UniModel})
-    LogDensityProblems.LogDensityOrder{0}()
+    return LogDensityProblems.LogDensityOrder{0}()
 end
 
 function LogDensityProblems.dimension(model::UniModel)
-    1
+    return 1
 end
 
 @testset "multivariate samplers" begin
-    model   = UniModel([0.])
-    @testset for sampler in [
-        Slice(1),
-        SliceSteppingOut(1),
-        SliceDoublingOut(1),
-    ]
+    model = UniModel([0.0])
+    @testset for sampler in [Slice(1), SliceSteppingOut(1), SliceDoublingOut(1)]
         @testset "initialization" begin
             model  = UniModel([0.0])
-            θ, y  = MCMCTesting.sample_joint(Random.default_rng(), model)
+            θ, y   = MCMCTesting.sample_joint(Random.default_rng(), model)
             model′ = AbstractMCMC.LogDensityModel(@set model.y = y)
 
             θ0    = [1.0]
-            chain = sample(
-                model,
-                sampler,
-                10;
-                initial_params=θ0,
-                progress=false,
-            )
+            chain = sample(model, sampler, 10; initial_params=θ0, progress=false)
             @test first(chain).params == θ0
         end
 
@@ -80,7 +70,7 @@ end
 
         @testset "determinism" begin
             model  = UniModel([0.0])
-            θ, y  = MCMCTesting.sample_joint(Random.default_rng(), model)
+            θ, y   = MCMCTesting.sample_joint(Random.default_rng(), model)
             model′ = AbstractMCMC.LogDensityModel(@set model.y = y)
 
             rng           = StableRNG(1)
@@ -91,14 +81,14 @@ end
             rng           = StableRNG(1)
             _, init_state = AbstractMCMC.step(rng, model′, sampler; initial_params=copy(θ))
             transition, _ = AbstractMCMC.step(rng, model′, sampler, init_state)
-            θ′′            = transition.params
+            θ′′           = transition.params
             @test θ′ == θ′′
         end
 
         @testset "type stability $(type)" for type in [Float32, Float64]
-            rng   = Random.default_rng()
+            rng = Random.default_rng()
             model = UniModel([zero(type)])
-            θ, y  = MCMCTesting.sample_joint(Random.default_rng(), model)
+            θ, y = MCMCTesting.sample_joint(Random.default_rng(), model)
             model′ = AbstractMCMC.LogDensityModel(@set model.y = y)
 
             @test eltype(θ) == type
@@ -106,7 +96,7 @@ end
 
             _, init_state = AbstractMCMC.step(rng, model′, sampler; initial_params=copy(θ))
             transition, _ = AbstractMCMC.step(rng, model′, sampler, init_state)
-            θ′             = transition.params
+            θ′ = transition.params
 
             @test eltype(θ′) == type
         end
@@ -118,7 +108,7 @@ end
             n_mcmc_thin      = 10
             test             = ExactRankTest(n_samples, n_mcmc_steps, n_mcmc_thin)
 
-            model   = UniModel([0.])
+            model   = UniModel([0.0])
             subject = TestSubject(model, sampler)
             @test seqmcmctest(test, subject, 0.001, n_pvalue_samples; show_progress=false)
         end