diff --git a/ext/TuringDynamicHMCExt.jl b/ext/TuringDynamicHMCExt.jl
index 144f9e700..2c4bd0898 100644
--- a/ext/TuringDynamicHMCExt.jl
+++ b/ext/TuringDynamicHMCExt.jl
@@ -63,7 +63,7 @@ function DynamicPPL.initialstep(
 
     # Define log-density function.
     ℓ = DynamicPPL.LogDensityFunction(
-        model, DynamicPPL.getlogjoint, vi; adtype=spl.alg.adtype
+        model, DynamicPPL.getlogjoint_internal, vi; adtype=spl.alg.adtype
     )
 
     # Perform initial step.
@@ -73,14 +73,9 @@ function DynamicPPL.initialstep(
     steps = DynamicHMC.mcmc_steps(results.sampling_logdensity, results.final_warmup_state)
     Q, _ = DynamicHMC.mcmc_next_step(steps, results.final_warmup_state.Q)
 
-    # Update the variables.
-    vi = DynamicPPL.unflatten(vi, Q.q)
-    # TODO(DPPL0.37/penelopeysm): This is obviously incorrect. Fix this.
-    vi = DynamicPPL.setloglikelihood!!(vi, Q.ℓq)
-    vi = DynamicPPL.setlogprior!!(vi, 0.0)
-
     # Create first sample and state.
-    sample = Turing.Inference.Transition(model, vi)
+    vi = DynamicPPL.unflatten(vi, Q.q)
+    sample = Turing.Inference.Transition(model, vi, nothing)
     state = DynamicNUTSState(ℓ, vi, Q, steps.H.κ, steps.ϵ)
 
     return sample, state
@@ -99,12 +94,9 @@ function AbstractMCMC.step(
     steps = DynamicHMC.mcmc_steps(rng, spl.alg.sampler, state.metric, ℓ, state.stepsize)
     Q, _ = DynamicHMC.mcmc_next_step(steps, state.cache)
 
-    # Update the variables.
-    vi = DynamicPPL.unflatten(vi, Q.q)
-    vi = DynamicPPL.setlogp!!(vi, Q.ℓq)
-
     # Create next sample and state.
-    sample = Turing.Inference.Transition(model, vi)
+    vi = DynamicPPL.unflatten(vi, Q.q)
+    sample = Turing.Inference.Transition(model, vi, nothing)
     newstate = DynamicNUTSState(ℓ, vi, Q, state.metric, state.stepsize)
 
     return sample, newstate
diff --git a/ext/TuringOptimExt.jl b/ext/TuringOptimExt.jl
index 41a89d3b5..0f755988e 100644
--- a/ext/TuringOptimExt.jl
+++ b/ext/TuringOptimExt.jl
@@ -102,7 +102,7 @@ function Optim.optimize(
     options::Optim.Options=Optim.Options();
     kwargs...,
 )
-    f = Optimisation.OptimLogDensity(model, Optimisation.getlogjoint_without_jacobian)
+    f = Optimisation.OptimLogDensity(model, DynamicPPL.getlogjoint)
     init_vals = DynamicPPL.getparams(f.ldf)
     optimizer = Optim.LBFGS()
     return _map_optimize(model, init_vals, optimizer, options; kwargs...)
@@ -124,7 +124,7 @@ function Optim.optimize(
     options::Optim.Options=Optim.Options();
     kwargs...,
 )
-    f = Optimisation.OptimLogDensity(model, Optimisation.getlogjoint_without_jacobian)
+    f = Optimisation.OptimLogDensity(model, DynamicPPL.getlogjoint)
     init_vals = DynamicPPL.getparams(f.ldf)
     return _map_optimize(model, init_vals, optimizer, options; kwargs...)
 end
@@ -140,7 +140,7 @@ function Optim.optimize(
 end
 
 function _map_optimize(model::DynamicPPL.Model, args...; kwargs...)
-    f = Optimisation.OptimLogDensity(model, Optimisation.getlogjoint_without_jacobian)
+    f = Optimisation.OptimLogDensity(model, DynamicPPL.getlogjoint)
     return _optimize(f, args...; kwargs...)
 end
 
diff --git a/src/mcmc/Inference.jl b/src/mcmc/Inference.jl
index 299670a0d..327507252 100644
--- a/src/mcmc/Inference.jl
+++ b/src/mcmc/Inference.jl
@@ -17,8 +17,8 @@ using DynamicPPL:
     setindex!!,
     push!!,
     setlogp!!,
-    getlogp,
     getlogjoint,
+    getlogjoint_internal,
     VarName,
     getsym,
     getdist,
@@ -123,71 +123,94 @@ end
 ######################
 # Default Transition #
 ######################
-# Default
-getstats(t) = nothing
+getstats(::Any) = NamedTuple()
 
+# TODO(penelopeysm): Remove this abstract type by converting SGLDTransition,
+# SMCTransition, and PGTransition to Turing.Inference.Transition instead.
 abstract type AbstractTransition end
 
-struct Transition{T,F<:AbstractFloat,S<:Union{NamedTuple,Nothing}} <: AbstractTransition
+struct Transition{T,F<:AbstractFloat,N<:NamedTuple} <: AbstractTransition
     θ::T
-    lp::F # TODO: merge `lp` with `stat`
-    stat::S
-end
-
-Transition(θ, lp) = Transition(θ, lp, nothing)
-function Transition(model::DynamicPPL.Model, vi::AbstractVarInfo, t)
-    θ = getparams(model, vi)
-    lp = getlogjoint(vi)
-    return Transition(θ, lp, getstats(t))
-end
+    logprior::F
+    loglikelihood::F
+    stat::N
+
+    """
+        Transition(model::Model, vi::AbstractVarInfo, sampler_transition)
+
+    Construct a new `Turing.Inference.Transition` object using the outputs of a
+    sampler step.
+
+    Here, `vi` represents a VarInfo _for which the appropriate parameters have
+    already been set_. However, the accumulators (e.g. logp) may in general
+    have junk contents. The role of this method is to re-evaluate `model` and
+    thus set the accumulators to the correct values.
+
+    `sampler_transition` is the transition object returned by the sampler
+    itself and is only used to extract statistics of interest.
+    """
+    function Transition(model::DynamicPPL.Model, vi::AbstractVarInfo, sampler_transition)
+        vi = DynamicPPL.setaccs!!(
+            vi,
+            (
+                DynamicPPL.ValuesAsInModelAccumulator(true),
+                DynamicPPL.LogPriorAccumulator(),
+                DynamicPPL.LogLikelihoodAccumulator(),
+            ),
+        )
+        _, vi = DynamicPPL.evaluate!!(model, vi)
+
+        # Extract all the information we need
+        vals_as_in_model = DynamicPPL.getacc(vi, Val(:ValuesAsInModel)).values
+        logprior = DynamicPPL.getlogprior(vi)
+        loglikelihood = DynamicPPL.getloglikelihood(vi)
+
+        # Get additional statistics
+        stats = getstats(sampler_transition)
+        return new{typeof(vals_as_in_model),typeof(logprior),typeof(stats)}(
+            vals_as_in_model, logprior, loglikelihood, stats
+        )
+    end
 
-function metadata(t::Transition)
-    stat = t.stat
-    if stat === nothing
-        return (lp=t.lp,)
-    else
-        return merge((lp=t.lp,), stat)
+    function Transition(
+        model::DynamicPPL.Model,
+        untyped_vi::DynamicPPL.VarInfo{<:DynamicPPL.Metadata},
+        sampler_transition,
+    )
+        # Re-evaluating the model is unconscionably slow for untyped VarInfo. It's
+        # much faster to convert it to a typed varinfo first, hence this method.
+        # https://github.com/TuringLang/Turing.jl/issues/2604
+        return Transition(model, DynamicPPL.typed_varinfo(untyped_vi), sampler_transition)
     end
 end
 
-DynamicPPL.getlogjoint(t::Transition) = t.lp
-
-# Metadata of VarInfo object
-metadata(vi::AbstractVarInfo) = (lp=getlogjoint(vi),)
+function getstats_with_lp(t::Transition)
+    return merge(
+        t.stat,
+        (
+            lp=t.logprior + t.loglikelihood,
+            logprior=t.logprior,
+            loglikelihood=t.loglikelihood,
+        ),
+    )
+end
+function getstats_with_lp(vi::AbstractVarInfo)
+    return (
+        lp=DynamicPPL.getlogjoint(vi),
+        logprior=DynamicPPL.getlogprior(vi),
+        loglikelihood=DynamicPPL.getloglikelihood(vi),
+    )
+end
 
 ##########################
 # Chain making utilities #
 ##########################
 
-"""
-    getparams(model, t)
-
-Return a named tuple of parameters.
-"""
-getparams(model, t) = t.θ
-function getparams(model::DynamicPPL.Model, vi::DynamicPPL.VarInfo)
-    # NOTE: In the past, `invlink(vi, model)` + `values_as(vi, OrderedDict)` was used.
-    # Unfortunately, using `invlink` can cause issues in scenarios where the constraints
-    # of the parameters change depending on the realizations. Hence we have to use
-    # `values_as_in_model`, which re-runs the model and extracts the parameters
-    # as they are seen in the model, i.e. in the constrained space. Moreover,
-    # this means that the code below will work both of linked and invlinked `vi`.
-    # Ref: https://github.com/TuringLang/Turing.jl/issues/2195
-    # NOTE: We need to `deepcopy` here to avoid modifying the original `vi`.
-    return DynamicPPL.values_as_in_model(model, true, deepcopy(vi))
-end
-function getparams(
-    model::DynamicPPL.Model, untyped_vi::DynamicPPL.VarInfo{<:DynamicPPL.Metadata}
-)
-    # values_as_in_model is unconscionably slow for untyped VarInfo. It's
-    # much faster to convert it to a typed varinfo before calling getparams.
-    # https://github.com/TuringLang/Turing.jl/issues/2604
-    return getparams(model, DynamicPPL.typed_varinfo(untyped_vi))
+getparams(::DynamicPPL.Model, t::AbstractTransition) = t.θ
+function getparams(model::DynamicPPL.Model, vi::AbstractVarInfo)
+    t = Transition(model, vi, nothing)
+    return getparams(model, t)
 end
-function getparams(::DynamicPPL.Model, ::DynamicPPL.VarInfo{NamedTuple{(),Tuple{}}})
-    return Dict{VarName,Any}()
-end
-
 function _params_to_array(model::DynamicPPL.Model, ts::Vector)
     names_set = OrderedSet{VarName}()
     # Extract the parameter names and values from each transition.
@@ -203,7 +226,6 @@ function _params_to_array(model::DynamicPPL.Model, ts::Vector)
             iters = map(DynamicPPL.varname_and_value_leaves, keys(vals), values(vals))
             mapreduce(collect, vcat, iters)
         end
-
         nms = map(first, nms_and_vs)
         vs = map(last, nms_and_vs)
         for nm in nms
@@ -218,14 +240,9 @@ function _params_to_array(model::DynamicPPL.Model, ts::Vector)
     return names, vals
 end
 
-function get_transition_extras(ts::AbstractVector{<:VarInfo})
-    valmat = reshape([getlogjoint(t) for t in ts], :, 1)
-    return [:lp], valmat
-end
-
 function get_transition_extras(ts::AbstractVector)
-    # Extract all metadata.
-    extra_data = map(metadata, ts)
+    # Extract stats + log probabilities from each transition or VarInfo
+    extra_data = map(getstats_with_lp, ts)
     return names_values(extra_data)
 end
 
@@ -334,7 +351,7 @@ function AbstractMCMC.bundle_samples(
         vals = map(values(sym_to_vns)) do vns
             map(Base.Fix1(getindex, params), vns)
         end
-        return merge(NamedTuple(zip(keys(sym_to_vns), vals)), metadata(t))
+        return merge(NamedTuple(zip(keys(sym_to_vns), vals)), getstats_with_lp(t))
     end
 end
 
@@ -396,84 +413,4 @@ function DynamicPPL.get_matching_type(
     return Array{T,N}
 end
 
-##############
-# Utilities  #
-##############
-
-"""
-
-    transitions_from_chain(
-        [rng::AbstractRNG,]
-        model::Model,
-        chain::MCMCChains.Chains;
-        sampler = DynamicPPL.SampleFromPrior()
-    )
-
-Execute `model` conditioned on each sample in `chain`, and return resulting transitions.
-
-The returned transitions are represented in a `Vector{<:Turing.Inference.Transition}`.
-
-# Details
-
-In a bit more detail, the process is as follows:
-1. For every `sample` in `chain`
-   1. For every `variable` in `sample`
-      1. Set `variable` in `model` to its value in `sample`
-   2. Execute `model` with variables fixed as above, sampling variables NOT present
-      in `chain` using `SampleFromPrior`
-   3. Return sampled variables and log-joint
-
-# Example
-```julia-repl
-julia> using Turing
-
-julia> @model function demo()
-           m ~ Normal(0, 1)
-           x ~ Normal(m, 1)
-       end;
-
-julia> m = demo();
-
-julia> chain = Chains(randn(2, 1, 1), ["m"]); # 2 samples of `m`
-
-julia> transitions = Turing.Inference.transitions_from_chain(m, chain);
-
-julia> [Turing.Inference.getlogjoint(t) for t in transitions] # extract the logjoints
-2-element Array{Float64,1}:
- -3.6294991938628374
- -2.5697948166987845
-
-julia> [first(t.θ.x) for t in transitions] # extract samples for `x`
-2-element Array{Array{Float64,1},1}:
- [-2.0844148956440796]
- [-1.704630494695469]
-```
-"""
-function transitions_from_chain(
-    model::DynamicPPL.Model, chain::MCMCChains.Chains; kwargs...
-)
-    return transitions_from_chain(Random.default_rng(), model, chain; kwargs...)
-end
-
-function transitions_from_chain(
-    rng::Random.AbstractRNG,
-    model::DynamicPPL.Model,
-    chain::MCMCChains.Chains;
-    sampler=DynamicPPL.SampleFromPrior(),
-)
-    vi = Turing.VarInfo(model)
-
-    iters = Iterators.product(1:size(chain, 1), 1:size(chain, 3))
-    transitions = map(iters) do (sample_idx, chain_idx)
-        # Set variables present in `chain` and mark those NOT present in chain to be resampled.
-        DynamicPPL.setval_and_resample!(vi, chain, sample_idx, chain_idx)
-        model(rng, vi, sampler)
-
-        # Convert `VarInfo` into `NamedTuple` and save.
-        Transition(model, vi)
-    end
-
-    return transitions
-end
-
 end # module
diff --git a/src/mcmc/emcee.jl b/src/mcmc/emcee.jl
index cdbaf9fd0..076e61d7e 100644
--- a/src/mcmc/emcee.jl
+++ b/src/mcmc/emcee.jl
@@ -65,14 +65,14 @@ function AbstractMCMC.step(
     end
 
     # Compute initial transition and states.
-    transition = map(Base.Fix1(Transition, model), vis)
+    transition = [Transition(model, vi, nothing) for vi in vis]
 
     # TODO: Make compatible with immutable `AbstractVarInfo`.
     state = EmceeState(
         vis[1],
         map(vis) do vi
             vi = DynamicPPL.link!!(vi, model)
-            AMH.Transition(vi[:], DynamicPPL.getlogjoint(vi), false)
+            AMH.Transition(vi[:], DynamicPPL.getlogjoint_internal(vi), false)
         end,
     )
 
@@ -87,18 +87,17 @@ function AbstractMCMC.step(
     densitymodel = AMH.DensityModel(
         Base.Fix1(
             LogDensityProblems.logdensity,
-            DynamicPPL.LogDensityFunction(model, DynamicPPL.getlogjoint, vi),
+            DynamicPPL.LogDensityFunction(model, DynamicPPL.getlogjoint_internal, vi),
         ),
     )
 
     # Compute the next states.
-    states = last(AbstractMCMC.step(rng, densitymodel, spl.alg.ensemble, state.states))
+    t, states = AbstractMCMC.step(rng, densitymodel, spl.alg.ensemble, state.states)
 
     # Compute the next transition and state.
     transition = map(states) do _state
         vi = DynamicPPL.unflatten(vi, _state.params)
-        t = Transition(getparams(model, vi), _state.lp)
-        return t
+        return Transition(model, vi, t)
     end
     newstate = EmceeState(vi, states)
 
diff --git a/src/mcmc/ess.jl b/src/mcmc/ess.jl
index feb737a30..bbf900657 100644
--- a/src/mcmc/ess.jl
+++ b/src/mcmc/ess.jl
@@ -31,7 +31,7 @@ function DynamicPPL.initialstep(
         EllipticalSliceSampling.isgaussian(typeof(dist)) ||
             error("ESS only supports Gaussian prior distributions")
     end
-    return Transition(model, vi), vi
+    return Transition(model, vi, nothing), vi
 end
 
 function AbstractMCMC.step(
@@ -56,7 +56,7 @@ function AbstractMCMC.step(
     vi = DynamicPPL.unflatten(vi, sample)
     vi = DynamicPPL.setloglikelihood!!(vi, state.loglikelihood)
 
-    return Transition(model, vi), vi
+    return Transition(model, vi, nothing), vi
 end
 
 # Prior distribution of considered random variable
diff --git a/src/mcmc/external_sampler.jl b/src/mcmc/external_sampler.jl
index 992a2fb2d..200c59293 100644
--- a/src/mcmc/external_sampler.jl
+++ b/src/mcmc/external_sampler.jl
@@ -22,8 +22,6 @@ There are a few more optional functions which you can implement to improve the i
 - `Turing.Inference.isgibbscomponent(::MySampler)`: If you want your sampler to function as a component in Turing's Gibbs sampler, you should make this evaluate to `true`.
 
 - `Turing.Inference.requires_unconstrained_space(::MySampler)`: If your sampler requires unconstrained space, you should return `true`. This tells Turing to perform linking on the VarInfo before evaluation, and ensures that the parameter values passed to your sampler will always be in unconstrained (Euclidean) space.
-
-- `Turing.Inference.getlogp_external(external_transition, external_state)`: Tell Turing how to extract the log probability density associated with this transition (and state). If you do not specify these, Turing will simply re-evaluate the model with the parameters obtained from `getparams`, which can be inefficient. It is therefore recommended to store the log probability density in either the transition or the state (or both) and override this method.
 """
 struct ExternalSampler{S<:AbstractSampler,AD<:ADTypes.AbstractADType,Unconstrained} <:
        InferenceAlgorithm
@@ -85,27 +83,21 @@ function externalsampler(
     return ExternalSampler(sampler, adtype, Val(unconstrained))
 end
 
-"""
-    getlogp_external(external_transition, external_state)
-
-Get the log probability density associated with the external sampler's
-transition and state. Returns `missing` by default; in this case, an extra
-model evaluation will be needed to calculate the correct log density.
-"""
-getlogp_external(::Any, ::Any) = missing
-getlogp_external(mh::AdvancedMH.Transition, ::AdvancedMH.Transition) = mh.lp
-getlogp_external(hmc::AdvancedHMC.Transition, ::AdvancedHMC.HMCState) = hmc.stat.log_density
-
-struct TuringState{S,V1<:AbstractVarInfo,M,V}
+# TODO(penelopeysm): Can't we clean this up somehow?
+struct TuringState{S,V1,M,V}
     state::S
-    # Note that this varinfo has the correct parameters and logp obtained from
-    # the state, whereas `ldf.varinfo` will in general have junk inside it.
+    # Note that this varinfo must have the correct parameters set; but logp
+    # does not matter as it will be re-evaluated
     varinfo::V1
+    # Note that in general the VarInfo inside this LogDensityFunction will have
+    # junk parameters and logp. It only exists to provide structure
     ldf::DynamicPPL.LogDensityFunction{M,V}
 end
 
-varinfo(state::TuringState) = state.varinfo
-varinfo(state::AbstractVarInfo) = state
+# get_varinfo should return something from which the correct parameters can be
+# obtained, hence we use state.varinfo rather than state.ldf.varinfo
+get_varinfo(state::TuringState) = state.varinfo
+get_varinfo(state::AbstractVarInfo) = state
 
 getparams(::DynamicPPL.Model, transition::AdvancedHMC.Transition) = transition.z.θ
 function getparams(model::DynamicPPL.Model, state::AdvancedHMC.HMCState)
@@ -115,27 +107,6 @@ getstats(transition::AdvancedHMC.Transition) = transition.stat
 
 getparams(::DynamicPPL.Model, transition::AdvancedMH.Transition) = transition.params
 
-function make_updated_varinfo(
-    f::DynamicPPL.LogDensityFunction, external_transition, external_state
-)
-    # Set the parameters.
-    new_parameters = getparams(f.model, external_state)
-    new_varinfo = DynamicPPL.unflatten(f.varinfo, new_parameters)
-    # Set (or recalculate, if needed) the log density.
-    new_logp = getlogp_external(external_transition, external_state)
-    return if ismissing(new_logp)
-        last(DynamicPPL.evaluate!!(f.model, new_varinfo, f.context))
-    else
-        # TODO(DPPL0.37/penelopeysm) This is obviously wrong. Note that we
-        # have the same problem here as in HMC in that the sampler doesn't
-        # tell us about how logp is broken down into prior and likelihood.
-        # We should probably just re-evaluate unconditionally. A bit
-        # unfortunate.
-        DynamicPPL.setlogprior!!(new_varinfo, 0.0)
-        DynamicPPL.setloglikelihood!!(new_varinfo, new_logp)
-    end
-end
-
 # TODO: Do we also support `resume`, etc?
 function AbstractMCMC.step(
     rng::Random.AbstractRNG,
@@ -163,7 +134,7 @@ function AbstractMCMC.step(
 
     # Construct LogDensityFunction
     f = DynamicPPL.LogDensityFunction(
-        model, DynamicPPL.getlogjoint, varinfo; adtype=alg.adtype
+        model, DynamicPPL.getlogjoint_internal, varinfo; adtype=alg.adtype
     )
 
     # Then just call `AbstractMCMC.step` with the right arguments.
@@ -182,13 +153,10 @@ function AbstractMCMC.step(
         )
     end
 
-    # Get the parameters and log density, and set them in the varinfo.
-    new_varinfo = make_updated_varinfo(f, transition_inner, state_inner)
-
-    # Update the `state`
+    new_parameters = getparams(f.model, state_inner)
+    new_vi = DynamicPPL.unflatten(f.varinfo, new_parameters)
     return (
-        Transition(f.model, new_varinfo, transition_inner),
-        TuringState(state_inner, new_varinfo, f),
+        Transition(f.model, new_vi, transition_inner), TuringState(state_inner, new_vi, f)
     )
 end
 
@@ -207,12 +175,9 @@ function AbstractMCMC.step(
         rng, AbstractMCMC.LogDensityModel(f), sampler, state.state; kwargs...
     )
 
-    # Get the parameters and log density, and set them in the varinfo.
-    new_varinfo = make_updated_varinfo(f, transition_inner, state_inner)
-
-    # Update the `state`
+    new_parameters = getparams(f.model, state_inner)
+    new_vi = DynamicPPL.unflatten(f.varinfo, new_parameters)
     return (
-        Transition(f.model, new_varinfo, transition_inner),
-        TuringState(state_inner, new_varinfo, f),
+        Transition(f.model, new_vi, transition_inner), TuringState(state_inner, new_vi, f)
     )
 end
diff --git a/src/mcmc/gibbs.jl b/src/mcmc/gibbs.jl
index 58db29789..692748767 100644
--- a/src/mcmc/gibbs.jl
+++ b/src/mcmc/gibbs.jl
@@ -343,7 +343,7 @@ struct GibbsState{V<:DynamicPPL.AbstractVarInfo,S}
     states::S
 end
 
-varinfo(state::GibbsState) = state.vi
+get_varinfo(state::GibbsState) = state.vi
 
 """
 Initialise a VarInfo for the Gibbs sampler.
@@ -390,7 +390,7 @@ function AbstractMCMC.step(
         initial_params=initial_params,
         kwargs...,
     )
-    return Transition(model, vi), GibbsState(vi, states)
+    return Transition(model, vi, nothing), GibbsState(vi, states)
 end
 
 function AbstractMCMC.step_warmup(
@@ -415,7 +415,7 @@ function AbstractMCMC.step_warmup(
         initial_params=initial_params,
         kwargs...,
     )
-    return Transition(model, vi), GibbsState(vi, states)
+    return Transition(model, vi, nothing), GibbsState(vi, states)
 end
 
 """
@@ -465,7 +465,7 @@ function gibbs_initialstep_recursive(
         initial_params=initial_params_local,
         kwargs...,
     )
-    new_vi_local = varinfo(new_state)
+    new_vi_local = get_varinfo(new_state)
     # Merge in any new variables that were introduced during the step, but that
     # were not in the domain of the current sampler.
     vi = merge(vi, get_global_varinfo(context))
@@ -493,7 +493,7 @@ function AbstractMCMC.step(
     state::GibbsState;
     kwargs...,
 )
-    vi = varinfo(state)
+    vi = get_varinfo(state)
     alg = spl.alg
     varnames = alg.varnames
     samplers = alg.samplers
@@ -503,7 +503,7 @@ function AbstractMCMC.step(
     vi, states = gibbs_step_recursive(
         rng, model, AbstractMCMC.step, varnames, samplers, states, vi; kwargs...
     )
-    return Transition(model, vi), GibbsState(vi, states)
+    return Transition(model, vi, nothing), GibbsState(vi, states)
 end
 
 function AbstractMCMC.step_warmup(
@@ -513,7 +513,7 @@ function AbstractMCMC.step_warmup(
     state::GibbsState;
     kwargs...,
 )
-    vi = varinfo(state)
+    vi = get_varinfo(state)
     alg = spl.alg
     varnames = alg.varnames
     samplers = alg.samplers
@@ -523,7 +523,7 @@ function AbstractMCMC.step_warmup(
     vi, states = gibbs_step_recursive(
         rng, model, AbstractMCMC.step_warmup, varnames, samplers, states, vi; kwargs...
     )
-    return Transition(model, vi), GibbsState(vi, states)
+    return Transition(model, vi, nothing), GibbsState(vi, states)
 end
 
 """
@@ -541,14 +541,11 @@ function setparams_varinfo!!(model, ::Sampler, state, params::AbstractVarInfo)
 end
 
 function setparams_varinfo!!(
-    model::DynamicPPL.Model,
-    sampler::Sampler{<:MH},
-    state::AbstractVarInfo,
-    params::AbstractVarInfo,
+    model::DynamicPPL.Model, sampler::Sampler{<:MH}, state::MHState, params::AbstractVarInfo
 )
-    # The state is already a VarInfo, so we can just return `params`, but first we need to
-    # update its logprob.
-    return last(DynamicPPL.evaluate!!(model, params))
+    # Re-evaluate to update the logprob.
+    new_vi = last(DynamicPPL.evaluate!!(model, params))
+    return MHState(new_vi, DynamicPPL.getlogjoint_internal(new_vi))
 end
 
 function setparams_varinfo!!(
@@ -569,7 +566,7 @@ function setparams_varinfo!!(
     params::AbstractVarInfo,
 )
     logdensity = DynamicPPL.LogDensityFunction(
-        model, DynamicPPL.getlogjoint, state.ldf.varinfo; adtype=sampler.alg.adtype
+        model, DynamicPPL.getlogjoint_internal, state.ldf.varinfo; adtype=sampler.alg.adtype
     )
     new_inner_state = setparams_varinfo!!(
         AbstractMCMC.LogDensityModel(logdensity), sampler, state.state, params
@@ -608,7 +605,7 @@ state for this sampler. This is relevant when multilple samplers are sampling th
 variables, and one might need it to be linked while the other doesn't.
 """
 function match_linking!!(varinfo_local, prev_state_local, model)
-    prev_varinfo_local = varinfo(prev_state_local)
+    prev_varinfo_local = get_varinfo(prev_state_local)
     was_linked = DynamicPPL.istrans(prev_varinfo_local)
     is_linked = DynamicPPL.istrans(varinfo_local)
     if was_linked && !is_linked
@@ -690,10 +687,10 @@ function gibbs_step_recursive(
     # Take a step with the local sampler.
     new_state = last(step_function(rng, conditioned_model, sampler, state; kwargs...))
 
-    new_vi_local = varinfo(new_state)
+    new_vi_local = get_varinfo(new_state)
     # Merge the latest values for all the variables in the current sampler.
     new_global_vi = merge(get_global_varinfo(context), new_vi_local)
-    new_global_vi = setlogp!!(new_global_vi, getlogp(new_vi_local))
+    new_global_vi = DynamicPPL.setlogp!!(new_global_vi, DynamicPPL.getlogp(new_vi_local))
 
     new_states = (new_states..., new_state)
     return gibbs_step_recursive(
diff --git a/src/mcmc/hmc.jl b/src/mcmc/hmc.jl
index 18733f6a8..d80502f7e 100644
--- a/src/mcmc/hmc.jl
+++ b/src/mcmc/hmc.jl
@@ -25,7 +25,7 @@ end
 ### Hamiltonian Monte Carlo samplers.
 ###
 
-varinfo(state::HMCState) = state.vi
+get_varinfo(state::HMCState) = state.vi
 
 """
     HMC(ϵ::Float64, n_leapfrog::Int; adtype::ADTypes.AbstractADType = AutoForwardDiff())
@@ -193,7 +193,7 @@ function DynamicPPL.initialstep(
     metricT = getmetricT(spl.alg)
     metric = metricT(length(theta))
     ldf = DynamicPPL.LogDensityFunction(
-        model, DynamicPPL.getlogjoint, vi; adtype=spl.alg.adtype
+        model, DynamicPPL.getlogjoint_internal, vi; adtype=spl.alg.adtype
     )
     lp_func = Base.Fix1(LogDensityProblems.logdensity, ldf)
     lp_grad_func = Base.Fix1(LogDensityProblems.logdensity_and_gradient, ldf)
@@ -208,9 +208,6 @@ function DynamicPPL.initialstep(
     end
     theta = vi[:]
 
-    # Cache current log density. We will reuse this if the transition is rejected.
-    logp_old = DynamicPPL.getlogp(vi)
-
     # Find good eps if not provided one
     if iszero(spl.alg.ϵ)
         ϵ = AHMC.find_good_stepsize(rng, hamiltonian, theta)
@@ -234,22 +231,13 @@ function DynamicPPL.initialstep(
         )
     end
 
-    # Update VarInfo based on acceptance
-    if t.stat.is_accept
-        vi = DynamicPPL.unflatten(vi, t.z.θ)
-        # Re-evaluate to calculate log probability density.
-        # TODO(penelopeysm): This seems a little bit wasteful. Unfortunately,
-        # even though `t.stat.log_density` contains some kind of logp, this
-        # doesn't track prior and likelihood separately but rather a single
-        # log-joint (and in linked space), so which we have no way to decompose
-        # this back into prior and likelihood. I don't immediately see how to
-        # solve this without re-evaluating the model.
-        _, vi = DynamicPPL.evaluate!!(model, vi)
+    # Update VarInfo parameters based on acceptance
+    new_params = if t.stat.is_accept
+        t.z.θ
     else
-        # Reset VarInfo back to its original state.
-        vi = DynamicPPL.unflatten(vi, theta)
-        vi = DynamicPPL.setlogp!!(vi, logp_old)
+        theta
     end
+    vi = DynamicPPL.unflatten(vi, new_params)
 
     transition = Transition(model, vi, t)
     state = HMCState(vi, 1, kernel, hamiltonian, t.z, adaptor)
@@ -293,9 +281,6 @@ function AbstractMCMC.step(
     vi = state.vi
     if t.stat.is_accept
         vi = DynamicPPL.unflatten(vi, t.z.θ)
-        # Re-evaluate to calculate log probability density.
-        # TODO(penelopeysm): This seems a little bit wasteful. See note above.
-        _, vi = DynamicPPL.evaluate!!(model, vi)
     end
 
     # Compute next transition and state.
@@ -308,7 +293,7 @@ end
 function get_hamiltonian(model, spl, vi, state, n)
     metric = gen_metric(n, spl, state)
     ldf = DynamicPPL.LogDensityFunction(
-        model, DynamicPPL.getlogjoint, vi; adtype=spl.alg.adtype
+        model, DynamicPPL.getlogjoint_internal, vi; adtype=spl.alg.adtype
     )
     lp_func = Base.Fix1(LogDensityProblems.logdensity, ldf)
     lp_grad_func = Base.Fix1(LogDensityProblems.logdensity_and_gradient, ldf)
diff --git a/src/mcmc/is.jl b/src/mcmc/is.jl
index 5f2f1627f..319e424fc 100644
--- a/src/mcmc/is.jl
+++ b/src/mcmc/is.jl
@@ -31,25 +31,19 @@ DynamicPPL.initialsampler(sampler::Sampler{<:IS}) = sampler
 function DynamicPPL.initialstep(
     rng::AbstractRNG, model::Model, spl::Sampler{<:IS}, vi::AbstractVarInfo; kwargs...
 )
-    # Need to manually construct the Transition here because we only
-    # want to use the likelihood.
-    xs = Turing.Inference.getparams(model, vi)
-    lp = DynamicPPL.getloglikelihood(vi)
-    return Transition(xs, lp, nothing), nothing
+    return Transition(model, vi, nothing), nothing
 end
 
 function AbstractMCMC.step(
     rng::Random.AbstractRNG, model::Model, spl::Sampler{<:IS}, ::Nothing; kwargs...
 )
     vi = VarInfo(rng, model, spl)
-    xs = Turing.Inference.getparams(model, vi)
-    lp = DynamicPPL.getloglikelihood(vi)
-    return Transition(xs, lp, nothing), nothing
+    return Transition(model, vi, nothing), nothing
 end
 
 # Calculate evidence.
 function getlogevidence(samples::Vector{<:Transition}, ::Sampler{<:IS}, state)
-    return logsumexp(map(x -> x.lp, samples)) - log(length(samples))
+    return logsumexp(map(x -> x.loglikelihood, samples)) - log(length(samples))
 end
 
 function DynamicPPL.assume(rng, ::Sampler{<:IS}, dist::Distribution, vn::VarName, vi)
diff --git a/src/mcmc/mh.jl b/src/mcmc/mh.jl
index e0258fda4..eb5b3aa3e 100644
--- a/src/mcmc/mh.jl
+++ b/src/mcmc/mh.jl
@@ -153,6 +153,27 @@ function MH(model::Model; proposal_type=AMH.StaticProposal)
     return AMH.MetropolisHastings(priors)
 end
 
+"""
+    MHState(varinfo::AbstractVarInfo, logjoint_internal::Real)
+
+State for Metropolis-Hastings sampling.
+
+`varinfo` must have the correct parameters set inside it, but its other fields
+(e.g. accumulators, which track logp) can in general be missing or incorrect.
+
+`logjoint_internal` is the log joint probability of the model, evaluated using
+the parameters and linking status of `varinfo`. It should be equal to
+`DynamicPPL.getlogjoint_internal(varinfo)`. This information is returned by the
+MH sampler so we store this here to avoid re-evaluating the model
+unnecessarily.
+"""
+struct MHState{V<:AbstractVarInfo,L<:Real}
+    varinfo::V
+    logjoint_internal::L
+end
+
+get_varinfo(s::MHState) = s.varinfo
+
 #####################
 # Utility functions #
 #####################
@@ -297,14 +318,15 @@ end
 
 # Make a proposal if we don't have a covariance proposal matrix (the default).
 function propose!!(
-    rng::AbstractRNG, vi::AbstractVarInfo, model::Model, spl::Sampler{<:MH}, proposal
+    rng::AbstractRNG, prev_state::MHState, model::Model, spl::Sampler{<:MH}, proposal
 )
+    vi = prev_state.varinfo
     # Retrieve distribution and value NamedTuples.
     dt, vt = dist_val_tuple(spl, vi)
 
     # Create a sampler and the previous transition.
     mh_sampler = AMH.MetropolisHastings(dt)
-    prev_trans = AMH.Transition(vt, DynamicPPL.getlogjoint(vi), false)
+    prev_trans = AMH.Transition(vt, prev_state.logjoint_internal, false)
 
     # Make a new transition.
     spl_model = DynamicPPL.contextualize(
@@ -313,35 +335,35 @@ function propose!!(
     densitymodel = AMH.DensityModel(
         Base.Fix1(
             LogDensityProblems.logdensity,
-            DynamicPPL.LogDensityFunction(spl_model, DynamicPPL.getlogjoint, vi),
+            DynamicPPL.LogDensityFunction(spl_model, DynamicPPL.getlogjoint_internal, vi),
         ),
     )
     trans, _ = AbstractMCMC.step(rng, densitymodel, mh_sampler, prev_trans)
-
-    # TODO: Make this compatible with immutable `VarInfo`.
-    # Update the values in the VarInfo.
-    # TODO(DPPL0.37/penelopeysm): This is obviously incorrect. We need to
-    # re-evaluate the model.
+    # trans.params isa NamedTuple
     set_namedtuple!(vi, trans.params)
-    vi = DynamicPPL.setloglikelihood!!(vi, trans.lp)
-    return DynamicPPL.setlogprior!!(vi, 0.0)
+    # Here, `trans.lp` is equal to `getlogjoint_internal(vi)`. We don't know
+    # how to set this back inside vi (without re-evaluating). However, the next
+    # MH step will require this information to calculate the acceptance
+    # probability, so we return it together with vi.
+    return MHState(vi, trans.lp)
 end
 
 # Make a proposal if we DO have a covariance proposal matrix.
 function propose!!(
     rng::AbstractRNG,
-    vi::AbstractVarInfo,
+    prev_state::MHState,
     model::Model,
     spl::Sampler{<:MH},
     proposal::AdvancedMH.RandomWalkProposal,
 )
+    vi = prev_state.varinfo
     # If this is the case, we can just draw directly from the proposal
     # matrix.
     vals = vi[:]
 
     # Create a sampler and the previous transition.
     mh_sampler = AMH.MetropolisHastings(spl.alg.proposals)
-    prev_trans = AMH.Transition(vals, DynamicPPL.getlogjoint(vi), false)
+    prev_trans = AMH.Transition(vals, prev_state.logjoint_internal, false)
 
     # Make a new transition.
     spl_model = DynamicPPL.contextualize(
@@ -350,16 +372,17 @@ function propose!!(
     densitymodel = AMH.DensityModel(
         Base.Fix1(
             LogDensityProblems.logdensity,
-            DynamicPPL.LogDensityFunction(spl_model, DynamicPPL.getlogjoint, vi),
+            DynamicPPL.LogDensityFunction(spl_model, DynamicPPL.getlogjoint_internal, vi),
         ),
     )
     trans, _ = AbstractMCMC.step(rng, densitymodel, mh_sampler, prev_trans)
-
-    # TODO(DPPL0.37/penelopeysm): This is obviously incorrect. We need to
-    # re-evaluate the model.
+    # trans.params isa AbstractVector
     vi = DynamicPPL.unflatten(vi, trans.params)
-    vi = DynamicPPL.setloglikelihood!!(vi, trans.lp)
-    return DynamicPPL.setlogprior!!(vi, 0.0)
+    # Here, `trans.lp` is equal to `getlogjoint_internal(vi)`. We don't know
+    # how to set this back inside vi (without re-evaluating). However, the next
+    # MH step will require this information to calculate the acceptance
+    # probability, so we return it together with vi.
+    return MHState(vi, trans.lp)
 end
 
 function DynamicPPL.initialstep(
@@ -373,18 +396,18 @@ function DynamicPPL.initialstep(
     # just link everything before sampling.
     vi = maybe_link!!(vi, spl, spl.alg.proposals, model)
 
-    return Transition(model, vi), vi
+    return Transition(model, vi, nothing), MHState(vi, DynamicPPL.getlogjoint_internal(vi))
 end
 
 function AbstractMCMC.step(
-    rng::AbstractRNG, model::Model, spl::Sampler{<:MH}, vi::AbstractVarInfo; kwargs...
+    rng::AbstractRNG, model::Model, spl::Sampler{<:MH}, state::MHState; kwargs...
 )
     # Cases:
     # 1. A covariance proposal matrix
     # 2. A bunch of NamedTuples that specify the proposal space
-    vi = propose!!(rng, vi, model, spl, spl.alg.proposals)
+    new_state = propose!!(rng, state, model, spl, spl.alg.proposals)
 
-    return Transition(model, vi), vi
+    return Transition(model, new_state.varinfo, nothing), new_state
 end
 
 ####
diff --git a/src/mcmc/particle_mcmc.jl b/src/mcmc/particle_mcmc.jl
index f7e40f432..6959e22cc 100644
--- a/src/mcmc/particle_mcmc.jl
+++ b/src/mcmc/particle_mcmc.jl
@@ -146,13 +146,11 @@ end
 
 function SMCTransition(model::DynamicPPL.Model, vi::AbstractVarInfo, weight)
     theta = getparams(model, vi)
-    lp = DynamicPPL.getlogjoint(vi)
+    lp = DynamicPPL.getlogjoint_internal(vi)
     return SMCTransition(theta, lp, weight)
 end
 
-metadata(t::SMCTransition) = (lp=t.lp, weight=t.weight)
-
-DynamicPPL.getlogp(t::SMCTransition) = t.lp
+getstats_with_lp(t::SMCTransition) = (lp=t.lp, weight=t.weight)
 
 struct SMCState{P,F<:AbstractFloat}
     particles::P
@@ -317,17 +315,15 @@ struct PGState
     rng::Random.AbstractRNG
 end
 
-varinfo(state::PGState) = state.vi
+get_varinfo(state::PGState) = state.vi
 
 function PGTransition(model::DynamicPPL.Model, vi::AbstractVarInfo, logevidence)
     theta = getparams(model, vi)
-    lp = DynamicPPL.getlogjoint(vi)
+    lp = DynamicPPL.getlogjoint_internal(vi)
     return PGTransition(theta, lp, logevidence)
 end
 
-metadata(t::PGTransition) = (lp=t.lp, logevidence=t.logevidence)
-
-DynamicPPL.getlogp(t::PGTransition) = t.lp
+getstats_with_lp(t::PGTransition) = (lp=t.lp, logevidence=t.logevidence)
 
 function getlogevidence(samples, sampler::Sampler{<:PG}, state::PGState)
     return mean(x.logevidence for x in samples)
diff --git a/src/mcmc/prior.jl b/src/mcmc/prior.jl
index eadeaceb3..6d7463c2f 100644
--- a/src/mcmc/prior.jl
+++ b/src/mcmc/prior.jl
@@ -17,11 +17,7 @@ function AbstractMCMC.step(
         model, DynamicPPL.SamplingContext(rng, DynamicPPL.SampleFromPrior(), model.context)
     )
     _, vi = DynamicPPL.evaluate!!(sampling_model, VarInfo())
-    # Need to manually construct the Transition here because we only
-    # want to use the prior probability.
-    xs = Turing.Inference.getparams(model, vi)
-    lp = DynamicPPL.getlogprior(vi)
-    return Transition(xs, lp, nothing), nothing
+    return Transition(model, vi, nothing), nothing
 end
 
 DynamicPPL.default_chain_type(sampler::Prior) = MCMCChains.Chains
diff --git a/src/mcmc/sghmc.jl b/src/mcmc/sghmc.jl
index 6eeddfefa..5ca351643 100644
--- a/src/mcmc/sghmc.jl
+++ b/src/mcmc/sghmc.jl
@@ -58,19 +58,15 @@ function DynamicPPL.initialstep(
     vi::AbstractVarInfo;
     kwargs...,
 )
-    # Transform the samples to unconstrained space and compute the joint log probability.
+    # Transform the samples to unconstrained space.
     if !DynamicPPL.islinked(vi)
         vi = DynamicPPL.link!!(vi, model)
-        vi = last(DynamicPPL.evaluate!!(model, vi))
     end
 
     # Compute initial sample and state.
-    sample = Transition(model, vi)
+    sample = Transition(model, vi, nothing)
     ℓ = DynamicPPL.LogDensityFunction(
-        model,
-        vi,
-        DynamicPPL.SamplingContext(spl, DynamicPPL.DefaultContext());
-        adtype=spl.alg.adtype,
+        model, DynamicPPL.getlogjoint_internal, vi; adtype=spl.alg.adtype
     )
     state = SGHMCState(ℓ, vi, zero(vi[:]))
 
@@ -98,12 +94,11 @@ function AbstractMCMC.step(
     α = spl.alg.momentum_decay
     newv = (1 - α) .* v .+ η .* grad .+ sqrt(2 * η * α) .* randn(rng, eltype(v), length(v))
 
-    # Save new variables and recompute log density.
+    # Save new variables.
     vi = DynamicPPL.unflatten(vi, θ)
-    vi = last(DynamicPPL.evaluate!!(model, vi))
 
     # Compute next sample and state.
-    sample = Transition(model, vi)
+    sample = Transition(model, vi, nothing)
     newstate = SGHMCState(ℓ, vi, newv)
 
     return sample, newstate
@@ -200,13 +195,11 @@ end
 
 function SGLDTransition(model::DynamicPPL.Model, vi::AbstractVarInfo, stepsize)
     theta = getparams(model, vi)
-    lp = DynamicPPL.getlogjoint(vi)
+    lp = DynamicPPL.getlogjoint_internal(vi)
     return SGLDTransition(theta, lp, stepsize)
 end
 
-metadata(t::SGLDTransition) = (lp=t.lp, SGLD_stepsize=t.stepsize)
-
-DynamicPPL.getlogp(t::SGLDTransition) = t.lp
+getstats_with_lp(t::SGLDTransition) = (lp=t.lp, SGLD_stepsize=t.stepsize)
 
 struct SGLDState{L,V<:AbstractVarInfo}
     logdensity::L
@@ -221,19 +214,15 @@ function DynamicPPL.initialstep(
     vi::AbstractVarInfo;
     kwargs...,
 )
-    # Transform the samples to unconstrained space and compute the joint log probability.
+    # Transform the samples to unconstrained space.
     if !DynamicPPL.islinked(vi)
         vi = DynamicPPL.link!!(vi, model)
-        vi = last(DynamicPPL.evaluate!!(model, vi))
     end
 
     # Create first sample and state.
     sample = SGLDTransition(model, vi, zero(spl.alg.stepsize(0)))
     ℓ = DynamicPPL.LogDensityFunction(
-        model,
-        vi,
-        DynamicPPL.SamplingContext(spl, DynamicPPL.DefaultContext());
-        adtype=spl.alg.adtype,
+        model, DynamicPPL.getlogjoint_internal, vi; adtype=spl.alg.adtype
     )
     state = SGLDState(ℓ, vi, 1)
 
@@ -252,9 +241,8 @@ function AbstractMCMC.step(
     stepsize = spl.alg.stepsize(step)
     θ .+= (stepsize / 2) .* grad .+ sqrt(stepsize) .* randn(rng, eltype(θ), length(θ))
 
-    # Save new variables and recompute log density.
+    # Save new variables.
     vi = DynamicPPL.unflatten(vi, θ)
-    vi = last(DynamicPPL.evaluate!!(model, vi))
 
     # Compute next sample and state.
     sample = SGLDTransition(model, vi, stepsize)
diff --git a/src/optimisation/Optimisation.jl b/src/optimisation/Optimisation.jl
index 29ea06726..19c52c381 100644
--- a/src/optimisation/Optimisation.jl
+++ b/src/optimisation/Optimisation.jl
@@ -43,113 +43,6 @@ Concrete type for maximum a posteriori estimation. Only used for the Optim.jl in
 """
 struct MAP <: ModeEstimator end
 
-# Most of these functions for LogPriorWithoutJacobianAccumulator are copied from
-# LogPriorAccumulator. The only one that is different is the accumulate_assume!! one.
-"""
-    LogPriorWithoutJacobianAccumulator{T} <: DynamicPPL.AbstractAccumulator
-
-Exactly like DynamicPPL.LogPriorAccumulator, but does not include the log determinant of the
-Jacobian of any variable transformations.
-
-Used for MAP optimisation.
-"""
-struct LogPriorWithoutJacobianAccumulator{T} <: DynamicPPL.AbstractAccumulator
-    logp::T
-end
-
-"""
-    LogPriorWithoutJacobianAccumulator{T}()
-
-Create a new `LogPriorWithoutJacobianAccumulator` accumulator with the log prior initialized to zero.
-"""
-LogPriorWithoutJacobianAccumulator{T}() where {T<:Real} =
-    LogPriorWithoutJacobianAccumulator(zero(T))
-function LogPriorWithoutJacobianAccumulator()
-    return LogPriorWithoutJacobianAccumulator{DynamicPPL.LogProbType}()
-end
-
-function Base.show(io::IO, acc::LogPriorWithoutJacobianAccumulator)
-    return print(io, "LogPriorWithoutJacobianAccumulator($(repr(acc.logp)))")
-end
-
-function DynamicPPL.accumulator_name(::Type{<:LogPriorWithoutJacobianAccumulator})
-    return :LogPriorWithoutJacobian
-end
-
-Base.copy(acc::LogPriorWithoutJacobianAccumulator) = acc
-
-function DynamicPPL.split(::LogPriorWithoutJacobianAccumulator{T}) where {T}
-    return LogPriorWithoutJacobianAccumulator(zero(T))
-end
-
-function DynamicPPL.combine(
-    acc::LogPriorWithoutJacobianAccumulator, acc2::LogPriorWithoutJacobianAccumulator
-)
-    return LogPriorWithoutJacobianAccumulator(acc.logp + acc2.logp)
-end
-
-function Base.:+(
-    acc1::LogPriorWithoutJacobianAccumulator, acc2::LogPriorWithoutJacobianAccumulator
-)
-    return LogPriorWithoutJacobianAccumulator(acc1.logp + acc2.logp)
-end
-
-function Base.zero(acc::LogPriorWithoutJacobianAccumulator)
-    return LogPriorWithoutJacobianAccumulator(zero(acc.logp))
-end
-
-function DynamicPPL.accumulate_assume!!(
-    acc::LogPriorWithoutJacobianAccumulator, val, logjac, vn, right
-)
-    return acc + LogPriorWithoutJacobianAccumulator(Distributions.logpdf(right, val))
-end
-function DynamicPPL.accumulate_observe!!(
-    acc::LogPriorWithoutJacobianAccumulator, right, left, vn
-)
-    return acc
-end
-
-function Base.convert(
-    ::Type{LogPriorWithoutJacobianAccumulator{T}}, acc::LogPriorWithoutJacobianAccumulator
-) where {T}
-    return LogPriorWithoutJacobianAccumulator(convert(T, acc.logp))
-end
-
-function DynamicPPL.convert_eltype(
-    ::Type{T}, acc::LogPriorWithoutJacobianAccumulator
-) where {T}
-    return LogPriorWithoutJacobianAccumulator(convert(T, acc.logp))
-end
-
-function getlogprior_without_jacobian(vi::DynamicPPL.AbstractVarInfo)
-    acc = DynamicPPL.getacc(vi, Val(:LogPriorWithoutJacobian))
-    return acc.logp
-end
-
-function getlogjoint_without_jacobian(vi::DynamicPPL.AbstractVarInfo)
-    return getlogprior_without_jacobian(vi) + DynamicPPL.getloglikelihood(vi)
-end
-
-# This is called when constructing a LogDensityFunction, and ensures the VarInfo has the
-# right accumulators.
-function DynamicPPL.ldf_default_varinfo(
-    model::DynamicPPL.Model, ::typeof(getlogprior_without_jacobian)
-)
-    vi = DynamicPPL.VarInfo(model)
-    vi = DynamicPPL.setaccs!!(vi, (LogPriorWithoutJacobianAccumulator(),))
-    return vi
-end
-
-function DynamicPPL.ldf_default_varinfo(
-    model::DynamicPPL.Model, ::typeof(getlogjoint_without_jacobian)
-)
-    vi = DynamicPPL.VarInfo(model)
-    vi = DynamicPPL.setaccs!!(
-        vi, (LogPriorWithoutJacobianAccumulator(), DynamicPPL.LogLikelihoodAccumulator())
-    )
-    return vi
-end
-
 """
     OptimLogDensity{
         M<:DynamicPPL.Model,
@@ -628,8 +521,10 @@ function estimate_mode(
 
     # Create an OptimLogDensity object that can be used to evaluate the objective function,
     # i.e. the negative log density.
-    getlogdensity =
-        estimator isa MAP ? getlogjoint_without_jacobian : DynamicPPL.getloglikelihood
+    # Note that we use `getlogjoint` rather than `getlogjoint_internal`: this
+    # is intentional, because even though the VarInfo may be linked, the
+    # optimisation target should not take the Jacobian term into account.
+    getlogdensity = estimator isa MAP ? DynamicPPL.getlogjoint : DynamicPPL.getloglikelihood
 
     # Set its VarInfo to the initial parameters.
     # TODO(penelopeysm): Unclear if this is really needed? Any time that logp is calculated
diff --git a/test/essential/container.jl b/test/essential/container.jl
index 8ebce270d..124637aab 100644
--- a/test/essential/container.jl
+++ b/test/essential/container.jl
@@ -28,14 +28,11 @@ using Turing
         @test trace.model.ctask.taped_globals.other === trace
 
         res = AdvancedPS.advance!(trace, false)
-        @test DynamicPPL.get_num_produce(trace.model.f.varinfo) == 1
         @test res ≈ -log(2)
 
         # Catch broken copy, espetially for RNG / VarInfo
         newtrace = AdvancedPS.fork(trace)
         res2 = AdvancedPS.advance!(trace)
-        @test DynamicPPL.get_num_produce(trace.model.f.varinfo) == 2
-        @test DynamicPPL.get_num_produce(newtrace.model.f.varinfo) == 1
     end
 
     @testset "fork" begin
diff --git a/test/mcmc/Inference.jl b/test/mcmc/Inference.jl
index 5d3d265c7..8c26e2a22 100644
--- a/test/mcmc/Inference.jl
+++ b/test/mcmc/Inference.jl
@@ -5,7 +5,7 @@ using ..NumericalTests: check_gdemo, check_numerical
 using Distributions: Bernoulli, Beta, InverseGamma, Normal
 using Distributions: sample
 import DynamicPPL
-using DynamicPPL: Sampler, getlogp
+using DynamicPPL: Sampler
 import ForwardDiff
 using LinearAlgebra: I
 import MCMCChains
@@ -116,11 +116,9 @@ using Turing
     @testset "Prior" begin
         N = 10_000
 
-        # Note that all chains contain 3 values per sample: 2 variables + log probability
         @testset "Single-threaded vanilla" begin
             chains = sample(StableRNG(seed), gdemo_d(), Prior(), N)
             @test chains isa MCMCChains.Chains
-            @test size(chains) == (N, 3, 1)
             @test mean(chains, :s) ≈ 3 atol = 0.11
             @test mean(chains, :m) ≈ 0 atol = 0.1
         end
@@ -128,7 +126,6 @@ using Turing
         @testset "Multi-threaded" begin
             chains = sample(StableRNG(seed), gdemo_d(), Prior(), MCMCThreads(), N, 4)
             @test chains isa MCMCChains.Chains
-            @test size(chains) == (N, 3, 4)
             @test mean(chains, :s) ≈ 3 atol = 0.11
             @test mean(chains, :m) ≈ 0 atol = 0.1
         end
@@ -139,8 +136,9 @@ using Turing
             )
             @test chains isa Vector{<:NamedTuple}
             @test length(chains) == N
-            @test all(length(x) == 3 for x in chains)
             @test all(haskey(x, :lp) for x in chains)
+            @test all(haskey(x, :logprior) for x in chains)
+            @test all(haskey(x, :loglikelihood) for x in chains)
             @test mean(x[:s][1] for x in chains) ≈ 3 atol = 0.11
             @test mean(x[:m][1] for x in chains) ≈ 0 atol = 0.1
         end
diff --git a/test/mcmc/ess.jl b/test/mcmc/ess.jl
index e918b3a51..1e1be9b45 100644
--- a/test/mcmc/ess.jl
+++ b/test/mcmc/ess.jl
@@ -54,7 +54,7 @@ using Turing
 
         @testset "gdemo with CSMC + ESS" begin
             alg = Gibbs(:s => CSMC(15), :m => ESS())
-            chain = sample(StableRNG(seed), gdemo(1.5, 2.0), alg, 2000)
+            chain = sample(StableRNG(seed), gdemo(1.5, 2.0), alg, 3_000)
             check_numerical(chain, [:s, :m], [49 / 24, 7 / 6]; atol=0.1)
         end
 
diff --git a/test/mcmc/external_sampler.jl b/test/mcmc/external_sampler.jl
index 6a6aebddb..5127c628e 100644
--- a/test/mcmc/external_sampler.jl
+++ b/test/mcmc/external_sampler.jl
@@ -21,7 +21,7 @@ function initialize_nuts(model::DynamicPPL.Model)
 
     # Create a LogDensityFunction
     f = DynamicPPL.LogDensityFunction(
-        model, DynamicPPL.getlogjoint, linked_vi; adtype=Turing.DEFAULT_ADTYPE
+        model, DynamicPPL.getlogjoint_internal, linked_vi; adtype=Turing.DEFAULT_ADTYPE
     )
 
     # Choose parameter dimensionality and initial parameter value
diff --git a/test/mcmc/gibbs.jl b/test/mcmc/gibbs.jl
index 092cc71f7..0fd76be3a 100644
--- a/test/mcmc/gibbs.jl
+++ b/test/mcmc/gibbs.jl
@@ -295,7 +295,7 @@ end
         vi::T
     end
 
-    Turing.Inference.varinfo(state::VarInfoState) = state.vi
+    Turing.Inference.get_varinfo(state::VarInfoState) = state.vi
     function Turing.Inference.setparams_varinfo!!(
         ::DynamicPPL.Model,
         ::DynamicPPL.Sampler,
@@ -312,8 +312,8 @@ end
         kwargs...,
     )
         spl.alg.non_warmup_init_count += 1
-        return Turing.Inference.Transition(nothing, 0.0),
-        VarInfoState(DynamicPPL.VarInfo(model))
+        vi = DynamicPPL.VarInfo(model)
+        return (Turing.Inference.Transition(model, vi, nothing), VarInfoState(vi))
     end
 
     function AbstractMCMC.step_warmup(
@@ -323,30 +323,30 @@ end
         kwargs...,
     )
         spl.alg.warmup_init_count += 1
-        return Turing.Inference.Transition(nothing, 0.0),
-        VarInfoState(DynamicPPL.VarInfo(model))
+        vi = DynamicPPL.VarInfo(model)
+        return (Turing.Inference.Transition(model, vi, nothing), VarInfoState(vi))
     end
 
     function AbstractMCMC.step(
         ::Random.AbstractRNG,
-        ::DynamicPPL.Model,
+        model::DynamicPPL.Model,
         spl::DynamicPPL.Sampler{<:WarmupCounter},
         s::VarInfoState;
         kwargs...,
     )
         spl.alg.non_warmup_count += 1
-        return Turing.Inference.Transition(nothing, 0.0), s
+        return Turing.Inference.Transition(model, s.vi, nothing), s
     end
 
     function AbstractMCMC.step_warmup(
         ::Random.AbstractRNG,
-        ::DynamicPPL.Model,
+        model::DynamicPPL.Model,
         spl::DynamicPPL.Sampler{<:WarmupCounter},
         s::VarInfoState;
         kwargs...,
     )
         spl.alg.warmup_count += 1
-        return Turing.Inference.Transition(nothing, 0.0), s
+        return Turing.Inference.Transition(model, s.vi, nothing), s
     end
 
     @model f() = x ~ Normal()
@@ -886,7 +886,9 @@ end
         function check_logp_correct(sampler)
             @testset "logp is set correctly" begin
                 @model logp_check() = x ~ Normal()
-                chn = sample(logp_check(), Gibbs(@varname(x) => sampler), 100)
+                chn = sample(
+                    logp_check(), Gibbs(@varname(x) => sampler), 100; progress=false
+                )
                 @test isapprox(logpdf.(Normal(), chn[:x]), chn[:lp])
             end
         end
diff --git a/test/mcmc/is.jl b/test/mcmc/is.jl
index 44fbe9201..2811e9c86 100644
--- a/test/mcmc/is.jl
+++ b/test/mcmc/is.jl
@@ -47,11 +47,11 @@ using Turing
 
         Random.seed!(seed)
         chain = sample(model, alg, n; check_model=false)
-        sampled = get(chain, [:a, :b, :lp])
+        sampled = get(chain, [:a, :b, :loglikelihood])
 
         @test vec(sampled.a) == ref.as
         @test vec(sampled.b) == ref.bs
-        @test vec(sampled.lp) == ref.logps
+        @test vec(sampled.loglikelihood) == ref.logps
         @test chain.logevidence == ref.logevidence
     end
 
diff --git a/test/optimisation/Optimisation.jl b/test/optimisation/Optimisation.jl
index 9909ee149..6b93e7629 100644
--- a/test/optimisation/Optimisation.jl
+++ b/test/optimisation/Optimisation.jl
@@ -41,12 +41,10 @@ using Turing
         @testset "With ConditionContext" begin
             m1 = model1(x)
             m2 = model2() | (x=x,)
-            # Doesn't matter if we use getlogjoint or getlogjoint_without_jacobian since the
+            # Doesn't matter if we use getlogjoint or getlogjoint_internal since the
             # VarInfo isn't linked.
-            ld1 = Turing.Optimisation.OptimLogDensity(
-                m1, Turing.Optimisation.getlogjoint_without_jacobian
-            )
-            ld2 = Turing.Optimisation.OptimLogDensity(m2, DynamicPPL.getlogjoint)
+            ld1 = Turing.Optimisation.OptimLogDensity(m1, DynamicPPL.getlogjoint)
+            ld2 = Turing.Optimisation.OptimLogDensity(m2, DynamicPPL.getlogjoint_internal)
             @test ld1(w) == ld2(w)
         end
 
@@ -54,22 +52,16 @@ using Turing
             vn = @varname(inner)
             m1 = prefix(model1(x), vn)
             m2 = prefix((model2() | (x=x,)), vn)
-            ld1 = Turing.Optimisation.OptimLogDensity(
-                m1, Turing.Optimisation.getlogjoint_without_jacobian
-            )
-            ld2 = Turing.Optimisation.OptimLogDensity(m2, DynamicPPL.getlogjoint)
+            ld1 = Turing.Optimisation.OptimLogDensity(m1, DynamicPPL.getlogjoint)
+            ld2 = Turing.Optimisation.OptimLogDensity(m2, DynamicPPL.getlogjoint_internal)
             @test ld1(w) == ld2(w)
         end
 
         @testset "Joint, prior, and likelihood" begin
             m1 = model1(x)
             a = [0.3]
-            ld_joint = Turing.Optimisation.OptimLogDensity(
-                m1, Turing.Optimisation.getlogjoint_without_jacobian
-            )
-            ld_prior = Turing.Optimisation.OptimLogDensity(
-                m1, Turing.Optimisation.getlogprior_without_jacobian
-            )
+            ld_joint = Turing.Optimisation.OptimLogDensity(m1, DynamicPPL.getlogjoint)
+            ld_prior = Turing.Optimisation.OptimLogDensity(m1, DynamicPPL.getlogprior)
             ld_likelihood = Turing.Optimisation.OptimLogDensity(
                 m1, DynamicPPL.getloglikelihood
             )