From 3f06635dd58932b9bc9af7219290dd83cc1d59bf Mon Sep 17 00:00:00 2001 From: Hong Ge Date: Wed, 19 Aug 2026 18:50:51 +0100 Subject: [PATCH 1/5] particle_mcmc: Make `delete_retained!` mutate the trace `PG`/`CSMC` posteriors were biased. `AdvancedPS.fork` marks a fork of the reference as no longer retained by calling `delete_retained!` for its side effect and discarding the return value, but our hook returned a fresh `TracedModel` instead of mutating, so `resample` stayed `false`. Every descendant of the reference therefore kept replaying the retained values in `tilde_assume!!` -- it was a copy of the reference, not a branch off it -- and the sweep lost the diversity that makes particle Gibbs valid. On a two-state HMM with ten observations, `PG(16)` state marginals were up to seven Monte Carlo standard errors from the exact forward-backward values, and are now within one. Present since v0.41.0. Co-Authored-By: Claude Code --- HISTORY.md | 7 +++++++ Project.toml | 2 +- src/mcmc/particle_mcmc.jl | 8 ++++++-- test/mcmc/particle_mcmc.jl | 24 +++++++++++++++++++++++- 4 files changed, 37 insertions(+), 4 deletions(-) diff --git a/HISTORY.md b/HISTORY.md index 580932cb71..f4c1a4aab7 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -1,3 +1,10 @@ +# 0.46.1 + +Fixed a bug that biased `PG` / `CSMC` posteriors. +Particles forked from the reference kept replaying the retained trajectory instead of sampling afresh, so every descendant of the reference was a copy of it and the sweep lost the diversity that makes particle Gibbs valid. +On a two-state HMM with ten observations, `PG(16)` state marginals sat up to seven Monte Carlo standard errors away from the exact forward-backward values; they are now within one. +The bug dates back to v0.41.0, and affects any model sampled with `PG` / `CSMC`, whether on its own or as a Gibbs component. + # 0.46.0 ## Breaking changes diff --git a/Project.toml b/Project.toml index e96c2dd443..0e09c7ad6c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "Turing" uuid = "fce5fe82-541a-59a6-adf8-730c64b5f9a0" -version = "0.46.0" +version = "0.46.1" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" diff --git a/src/mcmc/particle_mcmc.jl b/src/mcmc/particle_mcmc.jl index 2bef555edc..8cc05de5db 100644 --- a/src/mcmc/particle_mcmc.jl +++ b/src/mcmc/particle_mcmc.jl @@ -66,9 +66,13 @@ function AdvancedPS.delete_retained!(trace::TracedModel) # This method is called if, during a CSMC update, we perform a resampling # and choose the reference particle as the trajectory to carry on from. # In such a case, we need to ensure that when we continue sampling (i.e. - # the next time we hit tilde_assume!!), we don't use the values in the + # the next time we hit tilde_assume!!), we don't use the values in the # reference particle but rather sample new values. - return TracedModel(trace.model, trace.varinfo, true, trace.fargs, trace.kwargs) + # + # This has to mutate: `AdvancedPS.fork` calls it for its side effect and discards the + # return value. + trace.resample = true + return trace end function AdvancedPS.reset_model(trace::TracedModel) diff --git a/test/mcmc/particle_mcmc.jl b/test/mcmc/particle_mcmc.jl index 2cbc4dbece..a94c7c1b58 100644 --- a/test/mcmc/particle_mcmc.jl +++ b/test/mcmc/particle_mcmc.jl @@ -3,7 +3,7 @@ module ParticleMCMCTests using ..Models: gdemo_default using ..SamplerTestUtils: test_chain_logp_metadata using AdvancedPS: ResampleWithESSThreshold, resample_systematic, resample_multinomial -using Distributions: Bernoulli, Beta, Gamma, Normal, sample +using Distributions: Bernoulli, Beta, Gamma, Normal, Poisson, sample using FlexiChains: VNChain using Random: Random using StableRNGs: StableRNG @@ -168,6 +168,28 @@ end @test length(unique(c[:s])) == 1 end + @testset "conditional sweeps keep the population diverse" begin + # `k[t] ~ Poisson(1)` reweighted by `c^k[t]` is exactly `Poisson(c)`, since `e⁻¹cᵏ/k!` + # normalises to `e⁻ᶜcᵏ/k!`, and the reweighting is the only term carrying information: + # the sweep alone has to produce `E[k[t]] = c`. A descendant of the reference that + # replays the retained trajectory rather than branching off it is a copy of the + # reference, and the chain then over-visits that (high weight, hence high `k`) + # trajectory: `E[k]` used to come out between 2.4 and 2.7 across eight seeds, against + # within 0.05 of 2 once fixed. + c = 2.0 + @model function tilted_poisson(T, c) + k = Vector{Int}(undef, T) + for t in 1:T + k[t] ~ Poisson(1.0) + @addlogprob! k[t] * log(c) + end + end + chn = sample(StableRNG(468), tilted_poisson(4, c), PG(16), 2_000) + ks = reduce(vcat, collect(chn[@varname(k)])) + @test mean(ks) ≈ c atol = 0.15 + @test mean(iszero, ks) ≈ exp(-c) atol = 0.025 + end + @testset "addlogprob leads to reweighting" begin # Make sure that PG takes @addlogprob! into account. It didn't use to: # https://github.com/TuringLang/Turing.jl/issues/1996 From 30208a54cecac2244f50f4fcf373555dfd46093f Mon Sep 17 00:00:00 2001 From: Hong Ge Date: Wed, 19 Aug 2026 20:47:00 +0100 Subject: [PATCH 2/5] particle_mcmc: Test the reference against an exact posterior The previous test asserted on the `resample` flag through `AdvancedPS.fork`, which only a trace-based implementation has. This one goes through `sample` only: a Gibbs run over a switching model, compared against the exact posterior from enumerating all `2 * 2^8` configurations weighted by the model's own log density. It is also sensitive to the other way a reference can be wrong. Rebuilding it by replaying random numbers rather than reusing values drifts off the retained trajectory once the other Gibbs component re-conditions the latents, which needs the chain structure and the parameter placement this model has. Mean absolute error over the state marginals, across four seeds: 0.002 to 0.005 for a correct sweep, 0.020 to 0.028 for descendants that copy the reference, 0.017 to 0.022 for a replayed reference. The threshold is 0.01, and the testset costs about 50 seconds. Co-Authored-By: Claude Code --- HISTORY.md | 4 +-- src/mcmc/particle_mcmc.jl | 8 ++---- test/mcmc/particle_mcmc.jl | 58 +++++++++++++++++++++++++------------- 3 files changed, 43 insertions(+), 27 deletions(-) diff --git a/HISTORY.md b/HISTORY.md index f4c1a4aab7..2d9a0503d5 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -1,9 +1,7 @@ # 0.46.1 -Fixed a bug that biased `PG` / `CSMC` posteriors. -Particles forked from the reference kept replaying the retained trajectory instead of sampling afresh, so every descendant of the reference was a copy of it and the sweep lost the diversity that makes particle Gibbs valid. +Fixed a bug, present since v0.41.0, that biased `PG` / `CSMC` posteriors for any model, whether sampled on its own or as a Gibbs component. On a two-state HMM with ten observations, `PG(16)` state marginals sat up to seven Monte Carlo standard errors away from the exact forward-backward values; they are now within one. -The bug dates back to v0.41.0, and affects any model sampled with `PG` / `CSMC`, whether on its own or as a Gibbs component. # 0.46.0 diff --git a/src/mcmc/particle_mcmc.jl b/src/mcmc/particle_mcmc.jl index 8cc05de5db..ee4f86c325 100644 --- a/src/mcmc/particle_mcmc.jl +++ b/src/mcmc/particle_mcmc.jl @@ -66,11 +66,9 @@ function AdvancedPS.delete_retained!(trace::TracedModel) # This method is called if, during a CSMC update, we perform a resampling # and choose the reference particle as the trajectory to carry on from. # In such a case, we need to ensure that when we continue sampling (i.e. - # the next time we hit tilde_assume!!), we don't use the values in the - # reference particle but rather sample new values. - # - # This has to mutate: `AdvancedPS.fork` calls it for its side effect and discards the - # return value. + # the next time we hit tilde_assume!!), we don't use the values in the + # reference particle but rather sample new values. This has to mutate: + # `AdvancedPS.fork` calls it for its side effect and discards the return value. trace.resample = true return trace end diff --git a/test/mcmc/particle_mcmc.jl b/test/mcmc/particle_mcmc.jl index a94c7c1b58..a3a5c60736 100644 --- a/test/mcmc/particle_mcmc.jl +++ b/test/mcmc/particle_mcmc.jl @@ -3,7 +3,7 @@ module ParticleMCMCTests using ..Models: gdemo_default using ..SamplerTestUtils: test_chain_logp_metadata using AdvancedPS: ResampleWithESSThreshold, resample_systematic, resample_multinomial -using Distributions: Bernoulli, Beta, Gamma, Normal, Poisson, sample +using Distributions: Bernoulli, Beta, Categorical, Gamma, Normal, sample using FlexiChains: VNChain using Random: Random using StableRNGs: StableRNG @@ -168,26 +168,46 @@ end @test length(unique(c[:s])) == 1 end - @testset "conditional sweeps keep the population diverse" begin - # `k[t] ~ Poisson(1)` reweighted by `c^k[t]` is exactly `Poisson(c)`, since `e⁻¹cᵏ/k!` - # normalises to `e⁻ᶜcᵏ/k!`, and the reweighting is the only term carrying information: - # the sweep alone has to produce `E[k[t]] = c`. A descendant of the reference that - # replays the retained trajectory rather than branching off it is a copy of the - # reference, and the chain then over-visits that (high weight, hence high `k`) - # trajectory: `E[k]` used to come out between 2.4 and 2.7 across eight seeds, against - # within 0.05 of 2 once fixed. - c = 2.0 - @model function tilted_poisson(T, c) - k = Vector{Int}(undef, T) - for t in 1:T - k[t] ~ Poisson(1.0) - @addlogprob! k[t] * log(c) + @testset "conditional sweeps target the exact posterior" begin + # A reference particle that is not exactly the retained trajectory shows up here. The + # chain's stay probability is the other Gibbs component, so every `z[t]` is + # re-conditioned when `i` moves, and the observations are sharp enough to make the + # weights uneven; either is enough to bias the marginals. Measured over four seeds, a + # correct sweep keeps the mean error under 0.005, descendants that copy the reference + # rather than branching off it give 0.020 to 0.028, and a reference rebuilt by + # replaying random numbers instead of reusing values gives 0.017 to 0.022. + # + # Enumerating all `2 * 2^8` configurations and weighting them by the model's own log + # density keeps the target out of the hands of a reimplementation. + means, sd, stay = (-1.0, 1.0), 0.8, (0.35, 0.65) + @model function switching(y) + i ~ Categorical(2) + z = Vector{Int}(undef, length(y)) + z[1] ~ Categorical([0.5, 0.5]) + y[1] ~ Normal(means[z[1]], sd) + for t in 2:length(y) + p = stay[i] + z[t] ~ Categorical(z[t - 1] == 1 ? [p, 1 - p] : [1 - p, p]) + y[t] ~ Normal(means[z[t]], sd) end end - chn = sample(StableRNG(468), tilted_poisson(4, c), PG(16), 2_000) - ks = reduce(vcat, collect(chn[@varname(k)])) - @test mean(ks) ≈ c atol = 0.15 - @test mean(iszero, ks) ≈ exp(-c) atol = 0.025 + y = [-0.9163, -2.4106, -2.1881, 0.3716, 1.3404, -1.2046, -1.8294, -0.3521] + model = switching(y) + + confs = [ + (i, collect(z)) for i in 1:2 for z in Iterators.product(fill(1:2, length(y))...) + ] + w = exp.([logjoint(model, (; i=i, z=z)) for (i, z) in confs]) + w ./= sum(w) + exact = [ + sum(w[k] * (confs[k][2][t] - 1) for k in eachindex(w)) for t in eachindex(y) + ] + + alg = Gibbs(@varname(i) => MH(), @varname(z) => CSMC(8)) + chn = sample(StableRNG(468), model, alg, 6_000) + zs = stack(collect(z) for z in chn[@varname(z)]) + marginals = [mean(view(zs, t, :) .== 2) for t in eachindex(y)] + @test mean(abs, marginals .- exact) < 0.01 end @testset "addlogprob leads to reweighting" begin From 618c49799e62fb5b636f544c279e0aeb0487af32 Mon Sep 17 00:00:00 2001 From: Hong Ge Date: Wed, 19 Aug 2026 20:48:31 +0100 Subject: [PATCH 3/5] particle_mcmc: Cut the HISTORY entry to one line Co-Authored-By: Claude Code --- HISTORY.md | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/HISTORY.md b/HISTORY.md index 2d9a0503d5..861dc78de1 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -1,7 +1,6 @@ # 0.46.1 -Fixed a bug, present since v0.41.0, that biased `PG` / `CSMC` posteriors for any model, whether sampled on its own or as a Gibbs component. -On a two-state HMM with ten observations, `PG(16)` state marginals sat up to seven Monte Carlo standard errors away from the exact forward-backward values; they are now within one. +Fixed a bug, present since v0.41.0, that biased `PG` / `CSMC` posteriors, whether sampled on their own or as a Gibbs component. # 0.46.0 From 0af57e2562725552b1dfcba85fcaa870f0371d7a Mon Sep 17 00:00:00 2001 From: Hong Ge Date: Wed, 19 Aug 2026 20:59:29 +0100 Subject: [PATCH 4/5] Note the one-line HISTORY.md convention Co-Authored-By: Claude Code --- CLAUDE.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CLAUDE.md b/CLAUDE.md index 4922775aff..f8a12b64ac 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -98,3 +98,4 @@ User-facing functions accept `initial_params` as a convenience. `_convert_initia - Non-breaking changes target `main`; breaking changes target the `breaking` branch. - Julia ≥ 1.10.8 required (see `[compat]` in `Project.toml`). + - `HISTORY.md`: one line for a bugfix or internal change. Only a breaking change or new feature earns more, and only what a user needs to act on it: what broke, and the old → new form. The mechanism and any measurements belong in the commit and the PR. From cbe4e82fa8c82ea43304a56ce25413a209446b61 Mon Sep 17 00:00:00 2001 From: Hong Ge Date: Thu, 20 Aug 2026 12:59:31 +0100 Subject: [PATCH 5/5] particle_mcmc: Clarify the exact-posterior test The comment now states the property the test rests on -- conditional SMC is invariant only if the reference is exactly the retained path -- and the replay rule that makes the parameter's placement in the transition a requirement rather than a choice. In the model, the stay probability becomes a transition matrix indexed by the previous state, replacing a ternary over two probability vectors. The exact posterior now sums the parameter out in its own step, and reads the marginals off a flat vector of paths instead of indexing into tuples. Same numbers: the unfixed sampler still gives 0.0354. Co-Authored-By: Claude Code --- test/mcmc/particle_mcmc.jl | 43 +++++++++++++++++++------------------- 1 file changed, 22 insertions(+), 21 deletions(-) diff --git a/test/mcmc/particle_mcmc.jl b/test/mcmc/particle_mcmc.jl index a3a5c60736..2cecc39d75 100644 --- a/test/mcmc/particle_mcmc.jl +++ b/test/mcmc/particle_mcmc.jl @@ -169,44 +169,45 @@ end end @testset "conditional sweeps target the exact posterior" begin - # A reference particle that is not exactly the retained trajectory shows up here. The - # chain's stay probability is the other Gibbs component, so every `z[t]` is - # re-conditioned when `i` moves, and the observations are sharp enough to make the - # weights uneven; either is enough to bias the marginals. Measured over four seeds, a - # correct sweep keeps the mean error under 0.005, descendants that copy the reference - # rather than branching off it give 0.020 to 0.028, and a reference rebuilt by - # replaying random numbers instead of reusing values gives 0.017 to 0.022. + # Conditional SMC is invariant only if the reference is exactly the retained path, and + # this model is shaped so that either way of getting that wrong biases the marginals. + # The observations are sharp enough to keep the weights uneven, so resampling fires + # and a descendant of the reference that copies it rather than branching off shows up. + # The stay probability is the other Gibbs component, so a reference rebuilt by + # replaying random numbers -- `z[t] = z[t-1]` exactly when `u[t] < p` -- lands on a + # different path as soon as `i` moves. Mean absolute error over the marginals, across + # four seeds: under 0.005 for a correct sweep, 0.020 to 0.028 for the first failure, + # 0.017 to 0.022 for the second. # - # Enumerating all `2 * 2^8` configurations and weighting them by the model's own log - # density keeps the target out of the hands of a reimplementation. + # All `2 * 2^8` configurations enumerate the exact posterior, weighted by the model's + # own log density rather than by a reimplementation of it. means, sd, stay = (-1.0, 1.0), 0.8, (0.35, 0.65) @model function switching(y) i ~ Categorical(2) + p = stay[i] + transition = [p 1-p; 1-p p] z = Vector{Int}(undef, length(y)) z[1] ~ Categorical([0.5, 0.5]) y[1] ~ Normal(means[z[1]], sd) for t in 2:length(y) - p = stay[i] - z[t] ~ Categorical(z[t - 1] == 1 ? [p, 1 - p] : [1 - p, p]) + z[t] ~ Categorical(transition[z[t - 1], :]) y[t] ~ Normal(means[z[t]], sd) end end y = [-0.9163, -2.4106, -2.1881, 0.3716, 1.3404, -1.2046, -1.8294, -0.3521] model = switching(y) + T = length(y) - confs = [ - (i, collect(z)) for i in 1:2 for z in Iterators.product(fill(1:2, length(y))...) - ] - w = exp.([logjoint(model, (; i=i, z=z)) for (i, z) in confs]) - w ./= sum(w) - exact = [ - sum(w[k] * (confs[k][2][t] - 1) for k in eachindex(w)) for t in eachindex(y) - ] + paths = vec([collect(z) for z in Iterators.product(fill(1:2, T)...)]) + logws = [logjoint(model, (; i=i, z=path)) for i in 1:2, path in paths] + ws = exp.(logws .- maximum(logws)) + path_probs = vec(sum(ws; dims=1)) ./ sum(ws) # posterior over paths, `i` summed out + exact = [path_probs' * [path[t] == 2 for path in paths] for t in 1:T] alg = Gibbs(@varname(i) => MH(), @varname(z) => CSMC(8)) chn = sample(StableRNG(468), model, alg, 6_000) - zs = stack(collect(z) for z in chn[@varname(z)]) - marginals = [mean(view(zs, t, :) .== 2) for t in eachindex(y)] + draws = stack(collect(z) for z in chn[@varname(z)]) # T x ndraws + marginals = [mean(==(2), view(draws, t, :)) for t in 1:T] @test mean(abs, marginals .- exact) < 0.01 end