values of logp in chain do not include logabsdetjac term

[penelopeysm]

Minimal working example

Note that the original bug report was with AdvancedMH.RWMH. However, the same is true of any sampler that uses linked VarInfo. For example, Turing.NUTS exhibits the same behaviour.

using Turing, AdvancedMH, Test

@model f() = x ~ LogNormal()
m = f()
c = sample(m, externalsampler(AdvancedMH.RWMH(2); unconstrained = true), 10);

# Fails
@test isapprox(c[:lp], logpdf.(LogNormal(), c[:x]))

# Passes
@test isapprox(c[:lp], logpdf.(Normal(), log.(c[:x])))

Description

Follow-up from #2583.

Julia version info

julia> versioninfo()
Julia Version 1.11.6
Commit 9615af0f269 (2025-07-09 12:58 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin24.0.0)
  CPU: 10 × Apple M1 Pro
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, apple-m1)
Threads: 1 default, 0 interactive, 1 GC (on 8 virtual cores)

Manifest

]st --manifest

(Turing) pkg> st --manifest
Project Turing v0.39.5
Status `~/ppl/lib/Manifest.toml`
  [47edcb42] ADTypes v1.15.0
  [621f4979] AbstractFFTs v1.5.0
  [80f14c24] AbstractMCMC v5.6.3
  [7a57a42e] AbstractPPL v0.12.0
  [1520ce14] AbstractTrees v0.4.5
  [7d9f7c33] Accessors v0.1.42
  [79e6a3ab] Adapt v4.3.0
  [0bf59076] AdvancedHMC v0.8.0
  [5b7e9947] AdvancedMH v0.8.8
⌅ [576499cb] AdvancedPS v0.6.2
  [b5ca4192] AdvancedVI v0.4.1
  [66dad0bd] AliasTables v1.1.3
  [dce04be8] ArgCheck v2.5.0
  [4fba245c] ArrayInterface v7.19.0
  [13072b0f] AxisAlgorithms v1.1.0
  [39de3d68] AxisArrays v0.4.7
  [198e06fe] BangBang v0.4.4
  [9718e550] Baselet v0.1.1
  [76274a88] Bijectors v0.15.8
  [082447d4] ChainRules v1.72.5
  [d360d2e6] ChainRulesCore v1.25.2
  [0ca39b1e] Chairmarks v1.3.1
  [9e997f8a] ChangesOfVariables v0.1.10
  [861a8166] Combinatorics v1.0.3
  [38540f10] CommonSolve v0.2.4
  [bbf7d656] CommonSubexpressions v0.3.1
  [34da2185] Compat v4.17.0
  [a33af91c] CompositionsBase v0.1.2
  [88cd18e8] ConsoleProgressMonitor v0.1.2
  [187b0558] ConstructionBase v1.6.0
  [a8cc5b0e] Crayons v4.1.1
  [9a962f9c] DataAPI v1.16.0
  [864edb3b] DataStructures v0.18.22
  [e2d170a0] DataValueInterfaces v1.0.0
  [244e2a9f] DefineSingletons v0.1.2
  [8bb1440f] DelimitedFiles v1.9.1
  [b429d917] DensityInterface v0.4.0
  [163ba53b] DiffResults v1.1.0
  [b552c78f] DiffRules v1.15.1
  [a0c0ee7d] DifferentiationInterface v0.7.2
  [31c24e10] Distributions v0.25.120
  [ced4e74d] DistributionsAD v0.6.58
  [ffbed154] DocStringExtensions v0.9.5
  [366bfd00] DynamicPPL v0.36.15
  [cad2338a] EllipticalSliceSampling v2.0.0
  [4e289a0a] EnumX v1.0.5
  [e2ba6199] ExprTools v0.1.10
  [55351af7] ExproniconLite v0.10.14
  [7a1cc6ca] FFTW v1.9.0
  [9aa1b823] FastClosures v0.3.2
  [1a297f60] FillArrays v1.13.0
  [6a86dc24] FiniteDiff v2.27.0
  [f6369f11] ForwardDiff v1.0.1
  [069b7b12] FunctionWrappers v1.1.3
  [77dc65aa] FunctionWrappersWrappers v0.1.3
  [d9f16b24] Functors v0.5.2
  [46192b85] GPUArraysCore v0.2.0
  [34004b35] HypergeometricFunctions v0.3.28
  [22cec73e] InitialValues v0.3.1
  [a98d9a8b] Interpolations v0.16.1
  [8197267c] IntervalSets v0.7.11
  [3587e190] InverseFunctions v0.1.17
  [41ab1584] InvertedIndices v1.3.1
  [92d709cd] IrrationalConstants v0.2.4
  [c8e1da08] IterTools v1.10.0
  [82899510] IteratorInterfaceExtensions v1.0.0
  [692b3bcd] JLLWrappers v1.7.0
  [682c06a0] JSON v0.21.4
  [ae98c720] Jieko v0.2.1
  [5ab0869b] KernelDensity v0.6.10
  [5be7bae1] LBFGSB v0.4.1
  [8ac3fa9e] LRUCache v1.6.2
  [b964fa9f] LaTeXStrings v1.4.0
  [1d6d02ad] LeftChildRightSiblingTrees v0.2.0
⌅ [6f1fad26] Libtask v0.8.8
  [d3d80556] LineSearches v7.4.0
  [6fdf6af0] LogDensityProblems v2.1.2
  [996a588d] LogDensityProblemsAD v1.13.1
  [2ab3a3ac] LogExpFunctions v0.3.29
  [e6f89c97] LoggingExtras v1.1.0
  [c7f686f2] MCMCChains v7.1.0
  [be115224] MCMCDiagnosticTools v0.3.14
  [e80e1ace] MLJModelInterface v1.11.1
  [1914dd2f] MacroTools v0.5.16
  [dbb5928d] MappedArrays v0.4.2
  [128add7d] MicroCollections v0.2.0
  [e1d29d7a] Missings v1.2.0
  [2e0e35c7] Moshi v0.3.7
  [d41bc354] NLSolversBase v7.10.0
  [77ba4419] NaNMath v1.1.3
  [86f7a689] NamedArrays v0.10.4
  [c020b1a1] NaturalSort v1.0.0
  [6fe1bfb0] OffsetArrays v1.17.0
  [429524aa] Optim v1.13.2
  [3bd65402] Optimisers v0.4.6
  [7f7a1694] Optimization v4.4.0
  [bca83a33] OptimizationBase v2.8.0
  [36348300] OptimizationOptimJL v0.4.3
  [bac558e1] OrderedCollections v1.8.1
  [90014a1f] PDMats v0.11.35
  [d96e819e] Parameters v0.12.3
  [69de0a69] Parsers v2.8.3
  [85a6dd25] PositiveFactorizations v0.2.4
⌅ [aea7be01] PrecompileTools v1.2.1
  [21216c6a] Preferences v1.4.3
  [08abe8d2] PrettyTables v2.4.0
  [33c8b6b6] ProgressLogging v0.1.5
  [92933f4c] ProgressMeter v1.10.4
  [43287f4e] PtrArrays v1.3.0
  [1fd47b50] QuadGK v2.11.2
  [74087812] Random123 v1.7.1
  [e6cf234a] RandomNumbers v1.6.0
  [b3c3ace0] RangeArrays v0.3.2
  [c84ed2f1] Ratios v0.4.5
  [c1ae055f] RealDot v0.1.0
  [3cdcf5f2] RecipesBase v1.3.4
  [731186ca] RecursiveArrayTools v3.34.1
  [189a3867] Reexport v1.2.2
  [ae029012] Requires v1.3.1
  [79098fc4] Rmath v0.8.0
  [f2b01f46] Roots v2.2.8
  [7e49a35a] RuntimeGeneratedFunctions v0.5.15
  [26aad666] SSMProblems v0.5.0
⌃ [0bca4576] SciMLBase v2.102.1
  [c0aeaf25] SciMLOperators v1.3.1
  [53ae85a6] SciMLStructures v1.7.0
  [30f210dd] ScientificTypesBase v3.0.0
  [efcf1570] Setfield v1.1.2
  [a2af1166] SortingAlgorithms v1.2.1
  [9f842d2f] SparseConnectivityTracer v0.6.21
  [dc90abb0] SparseInverseSubset v0.1.2
  [0a514795] SparseMatrixColorings v0.4.21
  [276daf66] SpecialFunctions v2.5.1
  [171d559e] SplittablesBase v0.1.15
  [90137ffa] StaticArrays v1.9.13
  [1e83bf80] StaticArraysCore v1.4.3
⌃ [64bff920] StatisticalTraits v3.4.0
  [10745b16] Statistics v1.11.1
  [82ae8749] StatsAPI v1.7.1
  [2913bbd2] StatsBase v0.34.5
  [4c63d2b9] StatsFuns v1.5.0
  [892a3eda] StringManipulation v0.4.1
  [09ab397b] StructArrays v0.7.1
  [2efcf032] SymbolicIndexingInterface v0.3.41
  [3783bdb8] TableTraits v1.0.1
  [bd369af6] Tables v1.12.1
  [5d786b92] TerminalLoggers v0.1.7
  [28d57a85] Transducers v0.4.84
  [3a884ed6] UnPack v1.0.2
  [efce3f68] WoodburyMatrices v1.0.0
  [700de1a5] ZygoteRules v0.2.7
  [f5851436] FFTW_jll v3.3.11+0
  [1d5cc7b8] IntelOpenMP_jll v2025.0.4+0
  [81d17ec3] L_BFGS_B_jll v3.0.1+0
  [856f044c] MKL_jll v2025.0.1+1
  [efe28fd5] OpenSpecFun_jll v0.5.6+0
  [f50d1b31] Rmath_jll v0.5.1+0
  [1317d2d5] oneTBB_jll v2022.0.0+0
  [0dad84c5] ArgTools v1.1.2
  [56f22d72] Artifacts v1.11.0
  [2a0f44e3] Base64 v1.11.0
  [ade2ca70] Dates v1.11.0
  [8ba89e20] Distributed v1.11.0
  [f43a241f] Downloads v1.6.0
  [7b1f6079] FileWatching v1.11.0
  [9fa8497b] Future v1.11.0
  [b77e0a4c] InteractiveUtils v1.11.0
  [4af54fe1] LazyArtifacts v1.11.0
  [b27032c2] LibCURL v0.6.4
  [76f85450] LibGit2 v1.11.0
  [8f399da3] Libdl v1.11.0
  [37e2e46d] LinearAlgebra v1.11.0
  [56ddb016] Logging v1.11.0
  [d6f4376e] Markdown v1.11.0
  [a63ad114] Mmap v1.11.0
  [ca575930] NetworkOptions v1.2.0
  [44cfe95a] Pkg v1.11.0
  [de0858da] Printf v1.11.0
  [3fa0cd96] REPL v1.11.0
  [9a3f8284] Random v1.11.0
  [ea8e919c] SHA v0.7.0
  [9e88b42a] Serialization v1.11.0
  [1a1011a3] SharedArrays v1.11.0
  [6462fe0b] Sockets v1.11.0
  [2f01184e] SparseArrays v1.11.0
  [f489334b] StyledStrings v1.11.0
  [4607b0f0] SuiteSparse
  [fa267f1f] TOML v1.0.3
  [a4e569a6] Tar v1.10.0
  [8dfed614] Test v1.11.0
  [cf7118a7] UUIDs v1.11.0
  [4ec0a83e] Unicode v1.11.0
  [e66e0078] CompilerSupportLibraries_jll v1.1.1+0
  [deac9b47] LibCURL_jll v8.6.0+0
  [e37daf67] LibGit2_jll v1.7.2+0
  [29816b5a] LibSSH2_jll v1.11.0+1
  [c8ffd9c3] MbedTLS_jll v2.28.6+0
  [14a3606d] MozillaCACerts_jll v2023.12.12
  [4536629a] OpenBLAS_jll v0.3.27+1
  [05823500] OpenLibm_jll v0.8.5+0
  [bea87d4a] SuiteSparse_jll v7.7.0+0
  [83775a58] Zlib_jll v1.2.13+1
  [8e850b90] libblastrampoline_jll v5.11.0+0
  [8e850ede] nghttp2_jll v1.59.0+0
  [3f19e933] p7zip_jll v17.4.0+2
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated -m`

[penelopeysm]

I know what the problem here is. From the perspective of the externalsampler, it's running in linked space and thus it returns linked values that have (as far as the externalsampler is concerned) the correct logp.

The problem is that we then set those values back in the VarInfo using unflatten — which is OK as long as the VarInfo has been linked! since values_as_in_model will then extract the correct, unlinked, values — but if we take logp from the externalsampler it is not correct.

There are a couple of possible solutions.

One is, since values_as_in_model will re-evaluate the model anyway, we can recalculate the correct logp during that evaluation.

I had hoped that we could manually calculate logabsdetjac and then subtract it. There is a function Bijectors.bijector(::Model) which is quite related, and it avoids re-evaluating the model if you provide a varinfo. The problem is that I think it doesn't deal with dynamic models correctly. For example, this sort of model will probably break it, because the stored distribution is based on an old value of b that isn't changed by unflatten:

@model function f()
    a ~ Normal()
    b ~ truncated(Normal(); lower = b)
end

[penelopeysm]

The amount of refactoring needed to get a simple logp value out of values_as_in_model in getparams is stunning. Grr.

[yebai]

We should probably store logp for constrained parameters (i.e. logp + logabsdetjac) in the returned chain.

For internal data bookkeeping, can we store

• logp_without_logabsdetjac
• logabsdetjac

together, then compute the logp on the fly if needed?

Also, does the new accumulator make extracting logp easier?

[penelopeysm]

Agreed! And indeed, I had wondered if we needed a separate accumulator for logjac. @mhauru maybe might have thoughts.

[mhauru]

I ended up implementing something very similar for MAP/MLE optimisation in https://github.com/TuringLang/Turing.jl/blob/7948bc55752b860be48e63cc3b66f6577ea7f72a/src/optimisation/Optimisation.jl, because optimisation needs to ignore logjac. So yeah, could write an accumulator just for logjac and use it to replace what I've put in for optimisation.

[penelopeysm]

SMC() is even more annoying, it just returns a constant lp of 0.0. The reason is because varinfo.logp[] isn't used at all (instead the logp is tracked via Libtask.produce and handled manually by AdvancedPS). I think there's some work here to either:

1. Unconditionally rerun the model once per MCMC step to obtain the correct logp
2. Implement some interface functions to indicate whether samplers return correct / wrong / no logp and only conditionally include lp in chain

[penelopeysm]

Current thinking:

• we need to unconditionally accumulate logjac i.e. there should be a default accumulator whose role is solely to accumulate logjac, which is 0 for unlinked variables and nonzero for linked variables. Note that tilde_assume!! always calculates logjac so there is no extra calculation needed for this.

note that the docstring needs to be incredibly clear about which direction of the logjac this is referring to. In DynamicPPL this is somewhat all over the place but it seems that accumulate_assume!! expects the inverse-transform logjac. I personally would prefer it to be the other way round but am not fussed as long as it is made VERY clear.

• The logprior accumulator doesn't need to be renamed, but it needs to be made clear e.g. in the docstring that logprior refers to the prior probability of potentially linked parameters, i.e., its value depends on whether the varinfo is linked. The same is true of the getlogprior and getlogjoint functions.

• loglikelihood is fine.

• add a new function getlogjac(vi) which returns the value of the logjac accumulator

• create additional functions getlogprior_unlinked and getlogjoint_unlinked which always refer to the "values as in model" / unlinked values, i.e., getlogprior_unlinked(vi) == getlogprior(vi) + getlogjac(vi) and likewise for the joint. This equation might have a sign error, need to be very careful with this (although it should be easy enough to determine the correct sign by testing). Names are open to bikeshedding.

Will present this at Monday's meeting.

(I have to say that the accumulator framework makes this so much easier!)

The main alternative would be to make getlogprior and getlogjoint always return the unlinked value. In its place there could be getlogprior_internal and getlogjoint_internal which return values depending on the linked status of the VarInfo.

In many ways I actually prefer this but it would mark a quite large shift in mentality from what came before this so I am not sure if it will trip other people up.

[mhauru]

I think I got confused somewhere there in the terminology. Would getlogprior return the value that includes the logjac if necessary? I would prefer that. E.g. if

vi = VarInfo(model)
vi = link!!(vi)
getlogprior(vi)

returns the statistically incorrect log prior, the one that doesn't understand the significance of linking, then that feels dangerous.

Also, while we are at it, I wonder if we should implement a simple abstract type to sit under AbstractAccumulator just to reduce code duplication for simple accumulators that hold one value (or more specifically one float). Could call it SimpleAccumulator or LogProbAccumulator. Just feels like there's a lot of trivial code being duplicated between LogPriorAcc, LogLikelihoodAcc, LogJacAcc, etc. To be decided when writing the code, see how much this would save and whether it's worth the abstraction, just leaving it here as an idea.

[penelopeysm]

Would getlogprior return the value that includes the logjac if necessary?

I would rather not. Reason is I consider linking to be an implementation detail. From a statistical point of view the logprior is the one in unlinked space. The user is never told about linking and never has to do it explicitly when sampling with NUTS.

But I'm not going to die on that hill.

I wonder if we should implement a simple abstract type to sit under AbstractAccumulator just to reduce code duplication for simple accumulators that hold one value (or more specifically one float).

Yup, sounds good to me.

[mhauru]

Reason is I consider linking to be an implementation detail. From a statistical point of view the logprior is the one in unlinked space. The user is never told about linking and never has to do it explicitly when sampling with NUTS.

Doesn't that make it even more important that if the user calls getlogprior(vi) they get the value that does not depend on whether some function that vi has passed through has, unbeknownst to the user, linked the vi or not?