Paper deep dive
Uncovering and Fixing Collider Bias in Bayesian PINNs
Michael Obermayr, Robert Peharz
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 10/10/2026, 2:23:52 AM
Summary
This paper identifies and analyzes 'collider bias' in Bayesian Physics-Informed Neural Networks (B-PINNs), where the standard modeling assumption of independent physics and trajectory parameters coupled via residuals induces systematic posterior bias. The authors propose a hierarchical chain model that correctly models the generative process (physics -> trajectories -> observations) but introduces a doubly intractable inference problem. They resolve this by discretizing the stochastic dynamics, allowing exact sampling via particle MCMC, and provide analytical approximations for the bias mechanisms.
Entities (10)
Relation Signals (8)
Bayesian Physics-Informed Neural Networks → suffersfrom → Collider Bias
confidence 97% · We show that this modeling choice can induce severe systematic bias in the posterior over physical parameters... This is the so-called collider bias.
Chain Model → avoids → Collider Bias
confidence 96% · The chain model does not suffer from this posterior bias... Experiments confirm the predicted bias and show that the chain formulation successfully avoids it.
Collider Bias → causedby → Collider Structure
confidence 95% · They are commonly formulated via a collider structure... We show that this modeling choice can induce severe systematic bias
Discretization → enables → Particle MCMC
confidence 94% · This challenge can be resolved by discretizing the underlying stochastic dynamics, after which the chain posterior can be sampled exactly with particle MCMC.
Chain Model → poses → Doubly Intractable
confidence 93% · The chain model does not suffer from this posterior bias, but it poses a harder, so-called doubly intractable, inference problem
Residual Metric Bias → partof → Collider Bias
confidence 92% · We identify two distinct mechanisms characterizing the collider bias... Mechanism 1: Residual metric bias.
Dissipation Bias → partof → Collider Bias
confidence 92% · We identify two distinct mechanisms characterizing the collider bias... Mechanism 2: Dissipation bias.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Bayesian physics-informed neural networks (B-PINNs) are a popular framework for parameter and state inference from sparse or noisy observations. They are commonly formulated via a collider structure, in which physical and trajectory parameters are assumed to be a priori independent and become coupled through virtual likelihoods on differential-equation residuals that enforce physical consistency. We show that this modeling choice can induce severe systematic bias in the posterior over physical parameters: even when the prior is favorably centered on the ground-truth parameters, the resulting posterior can drift away and concentrate far from them. As a remedy, we advocate a hierarchical chain model in which physics generates trajectories, which in turn generate observations. The chain model does not suffer from this posterior bias, but it poses a harder, so-called doubly intractable, inference problem due to a physics-dependent normalization constant. This challenge can be resolved by discretizing the underlying stochastic dynamics, after which the chain posterior can be sampled exactly with particle MCMC. We identify two distinct mechanisms characterizing the collider bias, derive analytical approximations of their magnitudes, and establish diagnostic criteria for predicting when standard B-PINNs remain reliable. Experiments confirm the predicted bias and show that the chain formulation successfully avoids it.
Tags
Links
- Source: https://arxiv.org/abs/2610.11737v1
- Canonical: https://arxiv.org/abs/2610.11737v1
Trouble viewing inline? Open PDF directly →
Full Text
114,819 characters extracted from source content.
Expand or collapse full text
Uncovering and Fixing Collider Bias in Bayesian PINNs Michael Obermayr Robert Peharz TU Graz TU Graz Abstract Bayesian physics-informed neural networks (B-PINNs) are a popular framework for parameter and state inference from sparse or noisy observations. They are commonly formulated via a collider structure, in which physical and trajectory parameters are assumed to be a priori independent and become coupled through virtual likelihoods on differential-equation residuals that enforce physical consistency. We show that this modeling choice can induce severe systematic bias in the posterior over physical parameters: even when the prior is favorably centered on the ground-truth parameters, the resulting posterior can drift away and concentrate far from them. As a remedy, we advocate a hierarchical chain model in which physics generates trajectories, which in turn generate observations. The chain model does not suffer from this posterior bias, but it poses a harder, so-called doubly intractable, inference problem due to a physics-dependent normalization constant. This challenge can be resolved by discretizing the underlying stochastic dynamics, after which the chain posterior can be sampled exactly with particle MCMC. We identify two distinct mechanisms characterizing the collider bias, derive analytical approximations of their magnitudes, and establish diagnostic criteria for predicting when standard B-PINNs remain reliable. Experiments confirm the predicted bias and show that the chain formulation successfully avoids it. 1 Introduction Inferring unknown physical parameters from dynamical systems is a fundamental challenge across science and engineering. Physics-informed neural networks (PINNs) [22] tackle this inverse problem by training neural trajectory surrogates to jointly minimize data mismatch and physical law violations. Under limited and noisy data, a Bayesian treatment naturally quantifies uncertainty over trajectories and parameters. Yet, their assumed dependency structure requires careful consideration, which is the focus of this paper. (a) Collider (B-PINN) (b) Chain (ours) Figure 1: Modeling choices for physics-informed Bayesian inference. (a) Collider (B-PINN): Physics parameters λ and trajectory parameters θ are a priori independent, interacting only when conditioned on equation residuals r. (b) Chain (ours): Physics λ directly generates trajectories θ, which in turn generate data y. Figure 2: Posterior marginals for the damped linear oscillator. System governed by x¨=−cmx˙−kmx+σrξ(t) x=- cm x- kmx+ _rξ(t) with white noise ξ(t)ξ(t), mass m, damping c and spring constant k, ground-truth parameters m,c,k=1.0,1.2,4.0\m,c,k\=\1.0,1.2,4.0\, residual scale fixed at σr=0.3 _r=0.3, and 40 noisy observations (σobs=0.02 _obs=0.02). Posteriors are sampled via NUTS (24,000 samples) for B-PINN and PMMH [2] (14,000 samples) for our proposed chain model. Although all three physical parameters are inferred jointly, only the ratios c/mc/m and k/mk/m are identifiable; we therefore display their corresponding marginal posteriors. Despite initialization at the ground truth and identical priors centered at the ground truth, only the chain model accurately recovers the true values. In contrast, B-PINN exhibits a systematic collider bias, underestimating dissipation (c/mc/m) and stiffness (k/mk/m). The bias is reproducible over resampled datasets. Full experimental details are provided in Section 5. Bayesian PINNs (B-PINNs) [29] assume the factorization pcollider(λ,θ,r,y)=p(λ)⏟PhysicsPriorp(θ)⏟Traj.Priorp(r|λ,θ)⏟ResidualLikelihoodp(y|θ)⏟Observ.Likelihood.p_collider(λ,θ,r,y)= p(λ)_ subarraycPhysics\\[-1.0pt] Prior subarray\, p(θ)_ subarraycTraj.\\[-1.0pt] Prior subarray\, p(r\,|\,λ,θ)_ subarraycResidual\\[-1.0pt] Likelihood subarray\, p(y\,|\,θ)_ subarraycObserv.\\[-1.0pt] Likelihood subarray. (1) Here p(λ)p(λ) and p(θ)p(θ) are the priors over physics and trajectory parameters, respectively, and p(y|θ)p(y\,|\,θ) is the observation likelihood. Residuals r serve as “virtual observations” of zero measuring physics violations at collocation points, typically under a narrow zero-mean Gaussian likelihood p(r|λ,θ)p(r\,|\,λ,θ). Thus, λ and θ are a priori independent and coupled only upon conditioning on r=r=0. In graphical-model terms [13], this forms a collider at r (Fig. 1(a)), a structure widely adopted in the literature [26, 14, 25, 20]. Alternatively, one might use a hierarchical model, in which physics λ is equipped with a prior, trajectories θ depend on λ and observations y depend on θ: pchain(λ,θ,y)∝p(λ)⏟PhysicsPriorp(θ|λ)⏟TrajectoryLikelihoodp(y|θ)⏟ObservationLikelihood.p_chain(λ,θ,y) ~ p(λ)_ subarraycPhysics\\[-1.0pt] Prior subarray~~ p(θ\,|\,λ)_ subarraycTrajectory\\[-1.0pt] Likelihood subarray~~ p(y\,|\,θ)_ subarraycObservation\\[-1.0pt] Likelihood subarray. (2) Here the dependency structure between λ and θ forms a chain (head-to-tail node) [13], as illustrated in Fig. 1(b). Unlike the collider (1), θ depends directly on λ through the trajectory likelihood p(θ|λ)p(θ\,|\,λ) rather than indirectly through residual likelihood p(r|λ,θ)p(r\,|\,λ,θ). Which of the two formulations is preferable? The collider offers practical convenience: its factors are straightforward Gaussian priors and likelihoods. By contrast, specifying p(θ|λ)p(θ\,|\,λ) in the chain is less straightforward; it requires describing a distribution over trajectory parameters θ that may vary substantially with the physics λ. An intuitive construction is to start from the prior p(θ)p(θ), similar as in the collider formulation, reweight it according to the exponential of the negative squared residual error, and renormalize the result. As detailed in Section 2, this renders the chain model (2) nearly identical to the collider model (1), up to a physics-dependent normalization constant ZλZ_λ whose calculation requires an intractable marginalization over θ. Despite the computational simplicity of the collider, the chain should be preferred. Conceptually, the residuals r in the collider formulation act merely as auxiliary variables to couple physics λ and trajectories θ, rather than quantities of substantive interest. The chain avoids this artificial construction by modeling trajectory dependence directly via p(θ|λ)p(θ\,|\,λ). Furthermore, forward-sampling the collider yields an implausible generative process—typically generating incompatible pairs of λ and θ that demand unlikely residuals. By contrast, the chain follows the natural hierarchy in which physics λ generates trajectories θ, which in turn generate observations y. Crucially, the collider systematically biases the posterior over physical parameters. To demonstrate this, we sample synthetic data from a damped linear oscillator with known ground truth, center prior modes at the true values, and sample the collider posterior with NUTS [12], using Fourier features to prevent spectral bias [27]. Despite these favorable conditions, the collider marginals drift substantially from the ground truth (Fig. 2), by about two posterior standard deviations in c/mc/m. This is alarming, as the collider model can return confidently wrong estimates of the underlying physics. By contrast, our proposed chain inference exhibits the expected behavior: its posterior concentrates at the ground truth. As shown in Section 2, this bias stems directly from the physics-dependent normalization constant ZλZ_λ, which systematically favors less restrictive physics in the collider setting. While the chain formulation eliminates this bias, it yields a challenging, doubly-intractable problem [17]. Despite the prevalence of collider-based approaches in physics-informed neural inference, this structural bias remains largely unaddressed. To our knowledge, the only exceptions are Alberts, Hao, and Bilionis [1, 10, 11], who employ a similar chain structure and trace the posterior bias to ZλZ_λ; however, they address the doubly-intractable problem via computationally demanding nested inference without further analyzing the bias. Still, their formulation has gained little traction in the community and the collider formulation remains the default in B-PINN practice [16, 24]. 2 Problem Formulation and Probabilistic Modeling We consider parameterized second-order ordinary differential equations (ODEs) governing a trajectory x(t)∈ℝdx(t) ^d: Mλ(x,x˙)x¨=Fλ(x,x˙,t),M_λ(x, x)\, x=F_λ(x, x,t), (3) where x˙ x and x¨ x denote the first and second time derivatives, Mλ(x,x˙)∈ℝd×dM_λ(x, x) ^d× d is an invertible generalized mass matrix, and Fλ(x,x˙,t)∈ℝdF_λ(x, x,t) ^d is a force field driving the system, both parameterized by unknown physical parameters λ∈ℝpλ ^p. Equation (3) is the general form of a broad class of second-order dynamical systems arising in classical mechanics, structural dynamics, and robotics [8, 18]. In practice, the trajectory x(t)x(t) is rarely observed directly and is therefore treated as latent. Instead, we only have access to sparse, noisy observations y=ykk=1Ky=\y_k\_k=1^K at times tkt_k, governed by a Gaussian measurement model yk=x(tk)+ϵk,ϵk∼(0,Σobs),y_k=x(t_k)+ _k, _k (0, _obs), (4) with known observation covariance Σobs _obs. Beyond the measurement, the dynamics themselves are typically uncertain: simplifications and unmodeled degrees of freedom inevitably perturb the ideal dynamics (3). We account for these discrepancies by introducing process noise on the velocity v:=x˙v:= x. With this, the phase-space state z:=(x⊤,v⊤)⊤∈ℝ2dz:=(x ,v ) ^2d evolves according to an Itô stochastic differential equation (SDE) dzt=[vtaλ(zt,t)]⏟=:fλ(zt,t)dt+[0Aλ(zt)−1Σr1/2]dWt,dz_t= bmatrixv_t\\ a_λ(z_t,t) bmatrix_=:\,f_λ(z_t,t)dt+ bmatrix0\\ A_λ(z_t)^-1 _r^1/2 bmatrixdW_t, (5) driven by standard Brownian motion Wt∈ℝdW_t ^d. Here, aλ(z,t):=Mλ(x,v)−1Fλ(x,v,t)a_λ(z,t):=M_λ(x,v)^-1F_λ(x,v,t) is the acceleration field and fλf_λ the drift. The physical disturbances are modeled as constant white noise with covariance Σr _r, and the invertible, possibly state- and parameter-dependent matrix Aλ(z)∈ℝd×dA_λ(z) ^d× d maps them into velocity increments. The choice of AλA_λ determines how process noise enters the dynamics: setting Aλ=IdA_λ=I_d injects white noise directly into the accelerations, whereas setting Aλ=MλA_λ=M_λ models the perturbations as stochastic forces. The latter is more realistic for many mechanical systems, where physical disturbances act as fluctuating generalized forces or torques. As we show in Section 3, parameter dependence in AλA_λ is one of the factors contributing to collider bias in the inferred physical parameters. With process noise present (Σr≻0 _r 0), particular physics parameters λ induce a distribution p(x|λ)p(x\,|\,λ) over trajectories rather than a single deterministic trajectory. Our primary goal is to infer the marginal posterior distribution over the physical parameters given sparse observations of the trajectory: p(λ|y)∝p(λ)∫p(y|x)p(x|λ)x.p(λ\,|\,y) p(λ) p(y\,|\,x)\,p(x\,|\,λ)\,dx. (6) Following standard B-PINN practice, we fix the residual scale Σr _r as a hyperparameter rather than inferring it alongside λ [29, 14]. 2.1 Physics-Informed Residuals and Energy To connect the trajectory distribution induced by the stochastic dynamics (5) with physics-informed inference, we express deviations from the governing dynamics through a residual energy. In physics-informed machine learning, the continuous trajectory is represented by a differentiable surrogate xθ(t)x_θ(t) parameterized by θ (e.g., neural network weights). At a set of collocation points t1<⋯<tNt_1<…<t_N, we define the physics residual rn(θ,λ):=Aλ(zθ(tn))(x¨θ(tn)−aλ(zθ(tn),tn)),r_n(θ,λ):=A_λ(z_θ(t_n)) ( x_θ(t_n)-a_λ(z_θ(t_n),t_n) ), (7) where zθ:=(xθ⊤,x˙θ⊤)⊤z_θ:=(x_θ , x_θ ) . To relate these pointwise residuals to the process-noise model, suppose the collocation points are uniformly spaced with step size h. Over an interval of width h, the Brownian increment driving (5) has covariance hΣrh _r. Equivalently, the corresponding average stochastic forcing over that interval has covariance Σr/h _r/h. This motivates modeling the residuals at the collocation points as Gaussian deviations with covariance Σr/h _r/h, yielding the quadratic residual energy E[θ,λ]=h2∑n=1Nrn(θ,λ)⊤Σr−1rn(θ,λ).E[θ,λ]= h2 _n=1^Nr_n(θ,λ) _r^-1r_n(θ,λ). (8) Thus, the stochastic dynamics induce a natural quadratic energy measuring the compatibility between a trajectory xθx_θ and physical parameters λ. This residual energy is the common ingredient underlying the probabilistic formulations considered next. 2.2 Probabilistic Models: Collider vs. Chain We can now define the two probabilistic models—collider (1) and chain (2)—in a parallel manner and using similar components. Both factorizations share the same prior p(λ)p(λ) over physical parameters and the same observational likelihood, p(y|θ)=∏k=1K(yk|xθ(tk),Σobs).p(y\,|\,θ)= _k=1^KN (y_k\,|\,x_θ(t_k), _obs ). (9) The key difference between them is how the residual energy (8) is used to couple the trajectory parameters θ and the physical parameters λ. Collider formulation. Standard Bayesian PINNs [29, 26, 14, 25] place an independent prior p(θ)p(θ) on the trajectory parameters and enforce the physical law through auxiliary residuals treated as pseudo-observations equipped with a narrow zero-mean Gaussian: p(r=|λ,θ)=∏n=1N(|rn(θ,λ),1hΣr)=1Zrexp(−E[θ,λ]), splitp(r=0\,|\,λ,θ)&= _n=1^NN (0 |r_n(θ,λ), 1h _r )\\ &= 1Z_r (-E[θ,λ] ), split (10) where the Gaussian normalization constant Zr=det(2πΣr/h)N/2Z_r= (2π _r/h)^N/2 depends on Σr _r, but is invariant to λ and θ. Because λ and θ are assumed a priori independent and interact only via conditioning on r=r=0, the residual acts as a collider node in the underlying directed acyclic graph (Fig. 1(a)). Chain formulation. In contrast, the hierarchical chain formulation follows the generative direction λ→θ→yλ→θ→ y (Fig. 1(b)) by defining a conditional distribution over trajectories, p(θ|λ)=1Zλexp(−E[θ,λ])p(θ).p(θ\,|\,λ)= 1Z_λ (-E[θ,λ] )p(θ). (11) Here, p(θ)p(θ) is interpreted as a reference density over the trajectory parametrization, which can encode preferences for certain trajectories. Unlike in the collider formulation, it is not a prior over θ; instead, the exponential factor tilts this reference density toward trajectories that are compatible with the λ. Normalizing the resulting conditional distribution requires the normalization constant Zλ=∫exp(−E[θ,λ])p(θ)θ,Z_λ= (-E[θ,λ] )p(θ)\,dθ, (12) which generally depends on the physical parameters λ. Model Difference and The Collider Bias. The two formulations differ only through their normalization factors and we can write: pcollider(λ,θ,r=,y)=ZλZrpchain(λ,θ,y).p_collider(λ,θ,r=0,y)= Z_λZ_rp_chain(λ,θ,y). (13) Since ZrZ_r is constant with respect to λ, this difference carries directly to the posterior over the physical parameters: Δ(λ):=logpcollider(λ|y,r=)pchain(λ|y)=logZλ+C, (λ):= p_collider(λ\,|\,y,r=0)p_chain(λ\,|\,y)= Z_λ+C, (14) where we call the model difference Δ(λ) (λ) the collider bias. It is the mechanism underlying the posterior drift observed in Fig. 2. The term bias is deliberate: Δ displaces the posterior from the data-generating parameters, so that the collider’s credible intervals miss the ground truth, not merely the chain posterior (Section 5). The partition function ZλZ_λ measures the amount of trajectory space, weighted by the reference density p(θ)p(θ), that is compatible with the dynamics under λ. The collider formulation therefore reweights the desired posterior by ZλZ_λ, favoring λ that admit a larger volume of low-residual trajectories, that is, less restrictive physics. 3 Analysis of the Collider Bias For B-PINNs with a neural trajectory surrogate, the collider bias Δ(λ) (λ) decomposes into two distinct mechanisms. In Section 4.2, we show that adopting a Markovian state-space latent model eliminates one mechanism and renders the other analytically tractable. Laplace approximation of the collider bias. Because the partition function (12) is analytically intractable for neural surrogates, we place a Gaussian prior p(θ)=(0,P−1)p(θ)=N(0,P^-1) on the weights and apply a Laplace expansion around the data-free mode θ∗(λ)θ^*(λ). The collider bias (14) evaluates to Δ(λ)≈−12logdet(Gλ+P)+C, (λ)≈- 12 (G_λ+P )+C, (15) where GλG_λ is the Gauss–Newton curvature of the residual energy (8) at θ∗θ^*. The residual constrains neff:=tr[Gλ(Gλ+P)−1]≤rankGλ≤min(Nd,dimθ)n_eff:=tr\! [G_λ(G_λ+P)^-1 ] _λ≤ (Nd, θ) directions. Evaluating this curvature along the nominal trajectory z(t):=zθ∗(t)z(t):=z_θ^*(t) over t∈[0,]t∈[0,T] decomposes (15) into: Δ(λ)≈−neffNd∑n=1Nlog|detAλ(z(tn))|⏟(1) Residual metric bias−12∫0∑k|Resk(t)|dt⏟(2) Dissipation bias+C, split (λ)≈& - n_effNd _n=1^N A_λ(z(t_n)) _(1) Residual metric bias\\ & - 12 _0^T _k _k(t) \,dt_(2) Dissipation bias\;+\;C, split (16) where sk(t)\s_k(t)\ are the eigenvalues of the drift Jacobian DfλDf_λ along z(t)z(t) (full derivation in Appendix A.1). Eq. (16) explains the origin and direction of the bias. Derivation Assumptions. This decomposition relies on four simplifying assumptions (Appendix A.1): (i) the surrogate closely satisfies the ODE at θ∗θ^*, with the prior penalty depending only weakly on λ; (i) GλG_λ is evaluated along a single nominal trajectory; (i) the residual loss dominates the prior along constrained directions (neff≈rankGλn_eff _λ locally constant); and (iv) the residual metric and drift decouple additively. We also treat AλA_λ and DfλDf_λ as piecewise constant in time, and assume the surrogate bandwidth Ω exceeds the system’s characteristic dynamical rates. Mechanism 1: Residual metric bias. Because ZrZ_r normalizes over residual space, it omits the Jacobian determinant |detAλ| A_λ of the acceleration-to-residual mapping, whereas ZλZ_λ accounts for it along all neffn_eff constrained directions. Consequently, acceleration residuals (Aλ=IdA_λ=I_d) avoid this bias, whereas force residuals (Aλ=MλA_λ=M_λ) add −neffNd∑n=1Nlog|detMλ(z(tn))|- n_effNd _n=1^N M_λ(z(t_n)) , shifting posterior mass toward vanishing mass determinants. The residual metric is dictated by where the noise enters the dynamics (Section 2.1); switching to acceleration residuals for convenience changes the noise model and is therefore not a free choice. This pathology worsens along drift-invariant parameter directions where jointly scaling masses and forces preserves drift. Because the bias grows with neffn_eff while the data only constrain scale via finite observations, the artifact can overwhelm the likelihood and collapse the posterior toward vanishing masses (Appendix A.2). Mechanism 2: Dissipation bias. Lacking Markovian temporal structure, the neural surrogate couples perturbations through the operator dt−Dfλ ddt-Df_λ nonlocally over [0,][0,T], generating the second term in (16). When all perturbations decay (Resk≤0Res_k≤ 0), the spectral sum reduces to phase-space divergence: ∑k|Resk|=−∇⋅fλ=−tr(∂aλ/∂v) _k _k =-∇· f_λ=-tr(∂ a_λ/∂ v). The dissipation bias simplifies to Δdiss(λ)=12∫0∇⋅fλ(z(t),t)t+(/Ω)+C, _diss(λ)= 12 _0^T∇· f_λ(z(t),t)\,dt+O(T/ )+C, (17) whose leading term matches the Onsager–Machlup path-action correction [19]. In dissipative systems (∇⋅fλ<0∇· f_λ<0), eΔdisse _diss spuriously penalizes dissipation, underestimating friction, damping, or drag. Crucially, this penalty scales as Θ() (T), so that the B-PINN’s error grows with the record length relative to its own posterior width (Section 5.2). In the presence of growing modes, it analogously penalizes stronger instability (Appendix A.1). 4 Tractable Inference in the Chain Model With a global surrogate, the chain model yields a doubly intractable posterior [17], usually requiring expensive inner-loop re-estimation of ZλZ_λ via variational inference or Langevin dynamics that often diverges (Appendix A.4). Instead of attempting global integration, we discretize the generative SDE [23]. This factorizes ZλZ_λ into analytic per-step normalizers (21), enabling asymptotically exact inference via particle MCMC. 4.1 The Chain Model as a State-Space Model Discretizing [0,][0,T] into T intervals of width h defines states zi=(xi⊤,vi⊤)⊤z_i=(x_i ,v_i ) at ti=iht_i=ih and applying a semi-implicit Euler–Maruyama scheme to (5) yields: vi=vi−1+haλ(zi−1,ti−1)+hAλ(zi−1)−1Σr1/2ξi,xi=xi−1+hvi, splitv_i&=v_i-1+h\,a_λ(z_i-1,t_i-1)+ h\,A_λ(z_i-1)^-1 _r^1/2 _i,\\ x_i&=x_i-1+h\,v_i, split (18) with independent increments ξi∼(0,Id) _i (0,I_d). This discretization induces a Markovian trajectory prior over θ=z0:Tθ=z_0:T: p(z0:T|λ)=p(z0)∏i=1TKλ(zi|zi−1),Kλ(zi|zi−1)=δ(xi−xi−1−hvi)⋅(vi|v¯i,hAλ−1ΣrAλ−⊤), splitp(z_0:T\,|\,λ)&=p(z_0) _i=1^TK_λ(z_i\,|\,z_i-1),\\ K_λ(z_i\,|\,z_i-1)&=δ\! (x_i-x_i-1-hv_i )\\ & ·N\! (v_i\, |\, v_i,\;h\,A_λ^-1 _rA_λ^- ), split (19) where v¯i:=vi−1+haλ(zi−1,ti−1) v_i:=v_i-1+h\,a_λ(z_i-1,t_i-1) and AλA_λ is evaluated at zi−1z_i-1. Unlike (11), this prior admits analytic, per-step normalization. Combined with the prior p(λ)p(λ) and observations at tk⊆tii=1T\t_k\ \t_i\_i=1^T, this completes the chain model (2) over (λ,z0:T)(λ,z_0:T) to target the exact marginal posterior p(λ|y)p(λ\,|\,y). 4.2 Exact Collider Bias on the State-Space Latent On the Markovian latent state θ=z0:Tθ=z_0:T, the collider bias can be evaluated in closed form. The velocity update in (18) corresponds to the discrete residual ri=Aλ(zi−1)(vi−vi−1h−aλ(zi−1,ti−1)),r_i=A_λ(z_i-1) ( v_i-v_i-1h-a_λ(z_i-1,t_i-1) ), (20) making the transition kernel (19) proportional to exp(−h2ri⊤Σr−1ri) (- h2r_i _r^-1r_i), which matches the summand of (8) for N=TN=T. Because rir_i is affine in viv_i, each one-step transition normalizes analytically: Zλ(zi−1):=∫exp(−h2ri⊤Σr−1ri)dvi=(2πh)d/2det(Σr)1/2|detAλ(zi−1)|. splitZ_λ(z_i-1)&:= (- h2r_i _r^-1r_i )dv_i\\ & := (2π h)^d/2 ( _r)^1/2 A_λ(z_i-1) . split (21) The drift aλa_λ acts merely as a mean shift and drops out of the integral, causing the dissipation bias to vanish. Substituting ∏i=1TZλ(zi−1) _i=1^TZ_λ(z_i-1) for ZλZ_λ in (13) and subtracting logZr=T2logdet(2πΣr/h) Z_r= T2 (2π _r/h) gives the exact collider bias: ΔMarkov(λ,z0:T−1)=−∑i=1Tlog|detAλ(zi−1)|+C. _Markov(λ,z_0:T-1)=- _i=1^T A_λ(z_i-1) +C. (22) Only the residual metric bias remains: under acceleration residuals (Aλ=IdA_λ=I_d), ΔMarkov _Markov is independent of λ and both posteriors coincide; under force residuals (Aλ=MλA_λ=M_λ), ΔMarkov _Markov biases the collider posterior toward vanishing mass determinants. 4.3 Exact Parameter Inference via Particle MCMC For nonlinear dynamics, the marginal likelihood p(y|λ)=∫p(z0:T|λ)p(y|z0:T)dz0:Tp(y\,|\,λ)= p(z_0:T\,|\,λ)\,p(y\,|\,z_0:T)\,dz_0:T lacks an analytic form. Rather than approximating it variationally [3], we exploit its Markovian structure, which admits sequential Monte Carlo estimation via a bootstrap particle filter [9]. We propagate P particles forward under (18), reweight by (yk|xtk,Σobs)N(y_k\,|\,x_t_k, _obs) at observation times tkt_k, and adaptively resample [15]. The resulting estimator p^(y|λ) p(y\,|\,λ) costs (PT)O(PT) and is strictly unbiased. By the pseudo-marginal property [2], substituting p^(y|λ) p(y\,|\,λ) into the Metropolis–Hastings acceptance ratio guarantees Particle Marginal Metropolis–Hastings (PMMH) targets the exact posterior p(λ|y)p(λ\,|\,y) for any P≥1P≥ 1, where P governs only chain mixing. Sampling a trajectory alongside each accepted λ′λ additionally yields exact joint samples from p(λ,z0:T|y)p(λ,z_0:T\,|\,y). 5 Experiments We test the analysis of Section 3 on three systems that isolate its bias mechanisms, and compare our inference scheme with nested inference. Table 1: The chain model recovers the truth in every system; the collider fails where Section 3 predicts. logα α: posterior mean of the log scale factor along the weakly identified drift-invariant direction (Eq. (37); correct: ≈0≈ 0). Other columns: posterior mean of an identified quantity divided by its true value (correct: 1.001.00). Bold: off by more than a factor of two. † : chains did not converge or stalled; median of the draws, for both columns of the system. ‡ : exact Kalman likelihood, since PMMH cannot reach this posterior. All parameters in Appendix A.3.3. Oscillator Double Pendulum Relativistic Particle Model (latent) logα α c/mc/m logα α m2/m1m_2/m_1 logα α c Acceleration residuals (Aλ=IdA_λ=I_d): metric bias inactive Chain (state space) −0.01-0.01 0.99 −0.03-0.03 0.99 −0.02-0.02 0.89 Collider (state space) −0.01-0.01 0.99 −0.03-0.03 0.99 −0.02-0.02 0.89 Chain (Fourier MLP) −0.00-0.00 1.00 — — — — Collider B-PINN (Fourier MLP) −0.00-0.00 0.61 −0.03†-0.03 0.94 +0.00+0.00 1.78 Force / torque residuals (Aλ=MλA_λ=M_λ): metric bias active Chain (state space) −0.01-0.01 1.01 +0.04+0.04 1.00 −0.18-0.18 0.89 Collider (state space) −42.42‡-42.42 27 −2.37†-2.37 0.30 −8.17†-8.17 4.1⋅4.1·10^3 Chain (Fourier MLP) +0.29+0.29 0.74 — — — — Collider B-PINN (Fourier MLP) −0.93-0.93 2.05 −2.79-2.79 0.54 −1.88-1.88 1.61 5.1 Systems and Setup Data are generated from the SDE (5) with isotropic process noise Σr=σr2Id _r= _r^2I_d fixed at ground truth; physics priors are log-normal with medians at the truth. Jointly scaling masses and forces by α leaves the drift invariant; this drift-invariant direction is only weakly identifiable via process noise. Damped oscillator (d=1d=1): mass m, stiffness k, linear damping c, drift mx¨=−cx˙−kxm x=-c x-kx, truth (m,c,k)=(1,1.2,4)(m,c,k)=(1,1.2,4), and σr=0.3 _r=0.3. Because only c/mc/m and k/mk/m enter the drift, the invariant direction is (m,c,k)→α(m,c,k)(m,c,k)→α(m,c,k). We observe position K=40K=40 times over =24T=24 (σobs=0.02 _obs=0.02). Dissipation bias is active (∇⋅fλ=−c/m<0∇\!·\!f_λ=-c/m<0). Double pendulum (d=2d=2): masses m1,m2m_1,m_2 on massless rods of lengths ℓ1,ℓ2 _1, _2, with vertical angles x=(ϑ1,ϑ2)x=( _1, _2). The drift is Mλ(x)x¨=Fλ(x,x˙)M_λ(x)\, x=F_λ(x, x) with state-dependent mass matrix MλM_λ and forces FλF_λ: Mλ=[(m1+m2)ℓ12μcosδμcosδm2ℓ22], M_λ= bmatrix(m_1+m_2) _1^2&μ δ\\ μ δ&m_2 _2^2 bmatrix, Fλ=−[μϑ˙22sinδ+(m1+m2)gℓ1sinϑ1−μϑ˙12sinδ+m2gℓ2sinϑ2], F_λ=- bmatrixμ _2^2 δ+(m_1+m_2)g _1 _1\\ -μ _1^2 δ+m_2g _2 _2 bmatrix, with δ:=ϑ1−ϑ2δ:= _1- _2, μ:=m2ℓ1ℓ2μ:=m_2 _1 _2, and g=9.81g=9.81. Truth: (m1,m2,ℓ1,ℓ2)=(1,0.5,1,0.8)(m_1,m_2, _1, _2)=(1,0.5,1,0.8), σr=0.02 _r=0.02. Linearity in masses yields the invariant direction (m1,m2)→α(m1,m2)(m_1,m_2)→α(m_1,m_2). Both angles are observed K=20K=20 times over =12T=12 (σobs=0.02 _obs=0.02 rad). The conservative divergence integrates to a bounded boundary term (38), but chaos yields a Θ() (T) surrogate determinant term (Appendix A.3.6). Relativistic particle (d=1d=1): rest mass m0m_0, charge q, accelerated by E(t)=1.81cos(πt)E(t)=1.81 (π t) up to 0.5c0.5c. Drift: m0γ3x¨=qE(t)m_0γ^3 x=qE(t) with Lorentz factor γ=(1−x˙2/c2)−1/2γ=(1- x^2/c^2)^-1/2; truth: (m0,q,c)=(1,1,1)(m_0,q,c)=(1,1,1), σr=0.02 _r=0.02. As only q/m0q/m_0 governs the drift, the invariant direction is (m0,q)→α(m0,q)(m_0,q)→α(m_0,q). Position is observed K=20K=20 times over =12T=12 (σobs=0.005 _obs=0.005). The divergence integrates to a bounded boundary term but alternates sign, so ∑k|Resk|=|∇⋅fλ| _k|Re\,s_k|=|∇\!·\!f_λ| grows with T (Appendix A.3.6). Evaluation grid. We evaluate two residual metrics—acceleration (Aλ=IdA_λ=I_d, metric bias off) and force/torque (Aλ=MλA_λ=M_λ, metric bias on)—across four inference schemes: Chain (state space): our method (Section 4; PMMH, 1,0001,000–2,0002,000 particles, 24,00024,000 iterations). Collider (state space): differs from Chain solely by omitting the normalizer. Collider B-PINN: standard model (Section 2.2) with a Fourier-feature MLP (dimθ≈1.4 θ≈ 1.4k, Nd≤720Nd≤ 720), jointly sampled over (θ,λ)(θ,λ) via NUTS (44 chains, 9,0009,000 steps). Chain (Fourier MLP): augments the B-PINN with the closed-form correction −Δ- (16) to restore ZλZ_λ. We run this only on the oscillator, where state-independent AλA_λ and ∇⋅fλ∇\!·\!f_λ allow closed-form evaluation at each NUTS step; in the other systems, state dependence makes it prohibitive (Appendix A.4). Details and code: Appendix A.3 and https://github.com/obermm/ColliderBias. Figure 3: Longer records amplify the dissipation bias (oscillator, acceleration residuals). Posterior mean of log(c/m) (c/m) minus its true value, with ±1± 1 posterior standard deviation. The chain stays near the truth on both latents, and its posterior contracts as T grows; on the state-space latent the collider coincides with it (open rings). The B-PINN drifts further from the truth while its posterior width stays constant, so its error grows from 1.551.55 to 2.672.67 standard deviations. 5.2 Collider Bias and Its Removal In Table 1, the chain model recovers the truth throughout: |logα|≤0.18| α|≤ 0.18 (≤1.1≤ 1.1 posterior SDs) and identified parameters within 1.31.3 SDs (Appendix A.3.3). The collider fails as predicted. Null condition. With acceleration residuals (Aλ=IdA_λ=I_d) on the Markovian state-space latent, both biases are inactive: Δ(λ) (λ) in Eq. (22) is constant, yielding bit-identical chain and collider draws. Mechanism 1: residual metric bias. Under force residuals, the metric term in Eq. (16) rewards small masses by adding −Tdlogα-Td α to the state-space log-density and ≈−nefflogα≈-n_eff α to the B-PINN, opposed along the invariant ray only by process noise (Appendix A.2). On state-space latents (T=1,248T=1,248 for the oscillator), this reward dominates: the Kalman posterior collapses to logα=−42.4 α=-42.4, and PMMH drifts toward zero mass until Euler integration fails. Measured Δ(λ) (λ) values match Eq. (22) to machine precision (Appendix A.3.8). On the B-PINN, a weaker reward (neff≈110n_eff≈ 110) is checked by noise, yet mass scales are still underestimated by factors of 2.52.5 to 1616. The damage is not confined to the drift-invariant direction: since detMλ∝ℓ12ℓ22 M_λ _1^2 _2^2 also rewards short rods, the pendulum B-PINN shrinks ℓ1 _1 and ℓ2 _2 to 0.140.14 and 0.270.27 of their true values, so removing the scale ambiguity by reparameterization would not remove the bias. Mechanism 2: dissipation bias. Omitting divergence (17) biases B-PINNs toward underestimating dissipation. On the oscillator with acceleration residuals, B-PINN underestimates damping (c/m=0.73c/m=0.73 vs. true 1.201.20; Fig. 2). This is systematic: across ten dataset resamples, B-PINN underestimates log(c/m) (c/m) by 0.42±0.030.42± 0.03 and lies below the truth on every dataset, with 90%90\% credible intervals that cover the truth in only 44 of 1010 cases, whereas our chain’s intervals always do (Appendix A.3.4). With σr _r matched to the data, the bias persists for four tested σr∈[0.05,1] _r∈[0.05,1] (Appendix A.3.5). For the relativistic particle, B-PINN overestimates the speed of light (c=1.78c=1.78 vs. 1.001.00): although the divergence integral vanishes, its absolute path integral does not, penalizing smaller c by 9.59.5 nats per e-fold (Appendix A.3.6). The pendulum is unassessed, as its B-PINN does not converge. Longer records amplify dissipation bias. As the penalty scales as Θ() (T), longer horizons (=12→48T=12→ 48) push B-PINN damping estimates further from the truth without posterior contraction (Fig. 3). Reweighting the chain posterior by eΔdisse _diss predicts the bias’s direction, growth, and 60–75% of its size (Appendix A.3.8). Restoring the normalizer on the surrogate. Chain (Fourier MLP) adds −Δ- to an otherwise identical B-PINN. Under acceleration residuals on the oscillator, it removes 97.5%97.5\% of the log(c/m) (c/m) bias. Under force residuals, where metric correction uses approximate neffn_eff, it eliminates two-thirds of the mass bias (overshooting slightly; Appendix A.3.7). Thus, the missing normalizer, not the neural surrogate, causes the pathology. Comparison with nested inference. PMMH is up to 13×13× faster than the B-PINN, depending on the system state size (Appendix A.3.1). We further benchmark PMMH against unmodified public reimplementations of NSVI [10] and NPSGLD [11] on Euler–Maruyama latents with closed-form ZλZ_λ (Appendix A.4). PMMH converges across all settings (0.510.51–6.716.71 s/ESS). NPSGLD fails to converge, requiring 21×21×–2,900×2,900× more compute per ESS while underestimating posterior SDs by up to 30×30×. NSVI converges to a biased damping ratio on the oscillator (c/m≈0.80c/m≈ 0.80 vs. 1.201.20) and raising its compute by a factor of 100 does not remove this bias. 6 Conclusions B-PINNs suffer from an overlooked structural defect: conditioning on virtual residual observations turns physical parameters λ and trajectory surrogates into colliders. This implicitly neglects the parameter-dependent normalizing constant ZλZ_λ, systematically biasing posterior inference. Even when initialized at the ground truth, standard B-PINNs underestimate damping and mass, and grow increasingly overconfident with longer observation records. We approximated two bias mechanisms in closed form, confirmed them across three systems, yet showed that analytical corrections can remain insufficient. To eliminate collider bias entirely, we reformulate the trajectory as a Markovian state-space model. This enables an asymptotically exact particle MCMC scheme that avoids both the computational overhead of nested Langevin dynamics [11] and the variational gap of nested VI [10]. When are B-PINNs safe? Collider B-PINNs yield unbiased posteriors if two null conditions hold simultaneously (Section 3; Appendix Table 2): (i) the residual metric AλA_λ is parameter-independent (e.g., measuring acceleration rather than force); and (i) the dissipation bias (17) is negligible, with Δdiss _diss remaining roughly constant across posterior support. Otherwise, practitioners must adopt a Markovian chain formulation or, when AλA_λ and ∇⋅fλ∇\!·\!f_λ are state-independent, apply the closed-form correction −Δ- directly. Limitations and Future Work. Our particle MCMC scheme targets only stochastic ODEs (Σr≻0 _r 0) in low-to-moderate dimensions; particle filters degenerate as σr→0 _r→ 0 and face weight collapse on finely discretized PDEs. Because collider bias persists in PDE systems, extending our framework to spatial domains will require scalable alternatives such as ensemble Kalman or score-based filtering. Furthermore, selecting Σr _r to prevent residual misspecification remains open. Ideally Σr _r would be inferred jointly with λ, but this introduces a third collider-bias mechanism that we leave to future work. Implications. Although standard B-PINNs rely on collider structures, the resulting bias has gone largely unnoticed, as they appear well-behaved when the null conditions approximately hold. When violated, MCMC still converges cleanly to confidently wrong posteriors whose error is undetectable without ground truth. Crucially, this distortion peaks in weakly identified regimes and long horizons with sparse data—precisely when Bayesian rigor matters most. Practitioners should therefore check these null conditions and prefer a chain formulation. Acknowledgements This research was funded by the Austrian Research Promotion Agency (FFG) via the project VENTUS (910263) under the FFG AI for Green program. References [1] A. Alberts and I. Bilionis (2023) Physics-informed information field theory for modeling physical systems with uncertainty quantification. Journal of Computational Physics 486, p. 112100. External Links: Document Cited by: §A.4, §1. [2] C. Andrieu, A. Doucet, and R. Holenstein (2010) Particle Markov chain Monte Carlo methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72 (3), p. 269–342. External Links: Document Cited by: Figure 2, §4.3. [3] C. Archambeau, M. Opper, Y. Shen, D. Cornford, and J. Shawe-Taylor (2007) Variational inference for diffusion processes. In Advances in Neural Information Processing Systems, J. Platt, D. Koller, Y. Singer, and S. Roweis (Eds.), Vol. 20. External Links: Link Cited by: §4.3. [4] E. Bingham, J. P. Chen, M. Jankowiak, F. Obermeyer, N. Pradhan, T. Karaletsos, R. Singh, P. A. Szerlip, P. Horsfall, and N. D. Goodman (2019) Pyro: deep universal probabilistic programming. Journal of Machine Learning Research 20 (28), p. 1–6. Cited by: §A.3.1. [5] M. Bolt (2026) Niff-replication: a JAX reimplementation of the Neural Information Field Filter. Note: https://github.com/BoltMaxwell/niff-replicationCommit dab4403 Cited by: §A.3.1, §A.4. [6] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang (2018) JAX: composable transformations of Python+NumPy programs. Note: https://github.com/jax-ml/jax Cited by: §A.3.1. [7] D. Dür and A. Bach (1978) The Onsager–Machlup function as Lagrangian for the most probable path of a diffusion process. Communications in Mathematical Physics 60 (2), p. 153–170. External Links: Document Cited by: §A.1. [8] H. Goldstein, C. Poole, and J. Safko (2002) Classical mechanics. Third edition, Addison Wesley. Cited by: §2. [9] N. J. Gordon, D. J. Salmond, and A. F. M. Smith (1993) Novel approach to nonlinear/non-Gaussian Bayesian state estimation. IEE Proceedings F (Radar and Signal Processing) 140 (2), p. 107–113. External Links: Document Cited by: §4.3. [10] K. Hao and I. Bilionis (2024) An information field theory approach to Bayesian state and parameter estimation in dynamical systems. Journal of Computational Physics 512, p. 113139. External Links: Document Cited by: 1st item, §A.4, §1, §5.2, §6. [11] K. Hao and I. Bilionis (2025) Neural information field filter. Mechanical Systems and Signal Processing 226, p. 112253. External Links: Document Cited by: 2nd item, §A.4, §1, §5.2, §6. [12] M. D. Hoffman and A. Gelman (2014) The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research 15 (47), p. 1593–1623. External Links: Link Cited by: §1. [13] D. Koller and N. Friedman (2009) Probabilistic graphical models: principles and techniques. MIT Press, Cambridge, MA. Cited by: §1, §1. [14] K. Linka, A. Schäfer, X. Meng, Z. Zou, G. E. Karniadakis, and E. Kuhl (2022) Bayesian physics-informed neural networks for real-world nonlinear dynamical systems. Computer Methods in Applied Mechanics and Engineering 402, p. 115346. External Links: Document Cited by: §1, §2.2, §2. [15] J. S. Liu and R. Chen (1998) Sequential Monte Carlo methods for dynamic systems. Journal of the American Statistical Association 93 (443), p. 1032–1044. External Links: Link Cited by: §4.3. [16] H. Moon, H. Cho, W. Demeke, B. Ryu, and S. Ryu (2025) Thermal conductivity estimation of thermoelectric materials with uncertainty quantification using Bayesian physics-informed neural networks. Note: arXiv:2510.16723 Cited by: §1. [17] I. Murray, Z. Ghahramani, and D. J. C. MacKay (2006) MCMC for doubly-intractable distributions. In Proceedings of the Twenty-Second Conference on Uncertainty in Artificial Intelligence, p. 359–366. Cited by: §1, §4. [18] R. M. Murray, Z. Li, and S. S. Sastry (1994) A mathematical introduction to robotic manipulation. CRC Press. Cited by: §2. [19] L. Onsager and S. Machlup (1953) Fluctuations and irreversible processes. Physical Review 91 (6), p. 1505–1512. External Links: Document Cited by: §A.1, §3. [20] A. Pensoneault and X. Zhu (2024) Efficient Bayesian physics informed neural networks for inverse problems via ensemble Kalman inversion. Journal of Computational Physics 508, p. 113006. External Links: Document Cited by: §1. [21] D. Phan, N. Pradhan, and M. Jankowiak (2019) Composable effects for flexible and accelerated probabilistic programming in NumPyro. Note: arXiv:1912.11554 Cited by: §A.3.1. [22] M. Raissi, P. Perdikaris, and G. E. Karniadakis (2019) Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics 378, p. 686–707. External Links: Document Cited by: §1. [23] S. Särkkä and A. Solin (2019) Applied stochastic differential equations. Institute of Mathematical Statistics Textbooks, Vol. 10, Cambridge University Press. External Links: Document Cited by: §4. [24] K. Shukla, Z. Zou, T. Kaeufer, M. Triantafyllou, and G. E. Karniadakis (2026) Uncertainty quantification in PINNs for turbulent flows: Bayesian inference and repulsive ensembles. Note: arXiv:2604.17156 Cited by: §1. [25] S. Stock, J. Stiasny, D. Babazadeh, C. Becker, and S. Chatzivasileiadis (2023) Bayesian physics-informed neural networks for robust system identification of power systems. In 2023 IEEE Belgrade PowerTech, p. 1–6. External Links: Document Cited by: §1, §2.2. [26] L. Sun and J. Wang (2020) Physics-constrained Bayesian neural network for fluid flow reconstruction with sparse and noisy data. Theoretical and Applied Mechanics Letters 10 (3), p. 161–169. External Links: Document Cited by: §1, §2.2. [27] S. Wang, H. Wang, and P. Perdikaris (2021) On the eigenvector bias of Fourier feature networks: from regression to solving multi-scale PDEs with physics-informed neural networks. Computer Methods in Applied Mechanics and Engineering 384, p. 113938. External Links: Document Cited by: §1. [28] E. Wong and M. Zakai (1965) On the convergence of ordinary integrals to stochastic integrals. The Annals of Mathematical Statistics 36 (5), p. 1560–1564. External Links: Document Cited by: §A.1. [29] L. Yang, X. Meng, and G. E. Karniadakis (2021) B-PINNs: Bayesian physics-informed neural networks for forward and inverse PDE problems with noisy data. Journal of Computational Physics 425, p. 109913. External Links: Document Cited by: §1, §2.2, §2. Appendix A SUPPLEMENTARY MATERIAL A.1 Laplace Approximation of the Collider Bias for Global Surrogates This appendix derives the two-term decomposition (16) for a global surrogate such as a neural network. Step 1 approximates ZλZ_λ by Laplace’s method, Step 2 identifies the two factors through which λ enters the resulting curvature, Steps 3 and 4 differentiate with respect to these factors, and Step 5 collects the terms; a worked example for the oscillator follows. Along drift-invariant directions, the metric bias admits exact statements, derived in Appendix A.2. Setup. We write the residual (7) in continuous time, rθ(t,λ):=Aλ(zθ(t))(v˙θ(t)−aλ(zθ(t),t)),vθ:=x˙θ,r_θ(t,λ):=A_λ(z_θ(t)) ( v_θ(t)-a_λ(z_θ(t),t) ), v_θ:= x_θ, (23) so that rn(θ,λ)=rθ(tn,λ)r_n(θ,λ)=r_θ(t_n,λ); since vθv_θ is the derivative of xθx_θ, only the acceleration equation can be violated. The surrogate carries the weight prior p(θ)=(0,P−1)p(θ)=N(0,P^-1) of Section 3, where P→0P→ 0 corresponds to a flat reference measure. Chain and collider share this prior and the energy (8) and differ only in their normalizers, ZλZ_λ in (12) and ZrZ_r in (10), so that Δ=logZλ−logZr+const = Z_λ- Z_r+const by (14). Only ZλZ_λ depends on λ; since Σr _r is fixed, logZr Z_r is a constant. Step 1: Laplace approximation. Let θ∗(λ)θ^*(λ) minimize E[θ,λ]−logp(θ)E[θ,λ]- p(θ), the negative logarithm of the data-free integrand of ZλZ_λ. A second-order expansion around θ∗θ^* gives logZλ≈−E[θ∗,λ]−12θ∗⊤Pθ∗−12logdetHλ+const,Hλ:=∇θ2E|θ∗+P, Z_λ≈-E[θ^*,λ]- 12θ^* Pθ^*- 12 H_λ+const, H_λ:=∇^2_θE |_θ^*+P, (24) where the constant (12logdetP 12 P) does not depend on λ. We drop the first two terms, assuming that the surrogate (nearly) solves the ODE, E[θ∗,λ]≈0E[θ^*,λ]≈ 0, and that θ∗θ^* follows an ODE solution rather than the prior mode, so that 12θ∗⊤Pθ∗ 12θ^* Pθ^* depends only weakly on λ. This is exact for linear, unforced dynamics and a surrogate with xθ≡0x_θ≡ 0 at θ=0θ=0, where θ∗=0θ^*=0 for every λ. Since the residuals vanish at θ∗θ^*, so does the part of ∇θ2E∇^2_θE that is proportional to them, which leaves Hλ≈Gλ+PH_λ≈ G_λ+P with the Gauss–Newton curvature Gλ:=h∑n=1N(∂rθ(tn,λ)∂θ)⊤Σr−1∂rθ(tn,λ)∂θ,logZλ≈−12logdet(Gλ+P)+const.G_λ:=h _n=1^N ( ∂ r_θ(t_n,λ)∂θ )^\! _r^-1\, ∂ r_θ(t_n,λ)∂θ, Z_λ≈- 12 (G_λ+P)+const. (25) Since logZr Z_r is constant, this gives Eq. (15); the remaining steps evaluate its first term. Step 2: Two factors of the curvature. Let B(t):=∂zθ(t)/∂θB(t):=∂ z_θ(t)/∂θ be the tangent map of the surrogate at θ∗θ^*, and Dfλ(t):=∂fλ/∂zDf_λ(t):=∂ f_λ/∂ z the drift Jacobian along the fitted trajectory z(t):=zθ∗(t)z(t):=z_θ^*(t). By the chain rule, ∂rθ(t,λ)∂θ≈Aλ(z(t))(LλB)v(t),Lλ:=dt−Dfλ(t), ∂ r_θ(t,λ)∂θ≈ A_λ(z(t))\,(L_λB)_v(t), L_λ:= ddt-Df_λ(t), (26) where (⋅)v(·)_v denotes the velocity block. The position block of LλBL_λB vanishes because vθ=x˙θv_θ= x_θ, and the term containing the derivative of AλA_λ multiplies the residual, which vanishes at θ∗θ^*. Stacking over the collocation points gives Gλ=⊤λ⊤λ⊤−1λλ,λ=diag(Aλ(z(t1)),…,Aλ(z(tN))),=h−1diag(Σr,…,Σr).G_λ=B L_λ A_λ \,S^-1\,A_λL_λB, _λ=diag (A_λ(z(t_1)),…,A_λ(z(t_N)) ), =h^-1diag( _r,…, _r). (27) The stacked tangent map B depends on λ only through θ∗(λ)θ^*(λ), which we hold fixed. Since S is fixed, λ enters GλG_λ through two factors: the residual metric λA_λ and the drift operator λL_λ. Since these factors are rectangular, logdet(Gλ+P) (G_λ+P) does not split into a sum over them; we therefore differentiate it with respect to each factor in turn, holding the other fixed. Step 3: Residual metric (Mechanism 1). Let Aλ(z)=κλA¯(z)A_λ(z)= _λ A(z) with a scalar κλ>0 _λ>0 and a λ-independent A¯ A, e.g., Aλ=mA_λ=m for the oscillator. Then the metric scale factors out of (27), Gλ=κλ2G~,G~:=⊤λ⊤¯⊤−1¯λ,G_λ= _λ^2\, G, G:=B L_λ A S^-1 AL_λB, (28) where ¯:=diag(A¯(z(t1)),…,A¯(z(tN))) A:=diag ( A(z(t_1)),…, A(z(t_N)) ), and G~ G does not depend on κλ _λ. Since ∂Gλ/∂logκλ=2Gλ∂ G_λ/∂ _λ=2G_λ, Jacobi’s formula gives ∂logκλ(−12logdet(Gλ+P))=−tr[(Gλ+P)−1Gλ]=−neff. ∂ _λ (- 12 (G_λ+P) )=-tr\! [(G_λ+P)^-1G_λ ]=-\,n_eff. (29) When is neffn_eff constant? For P≻0P 0, let gk\g_k\ be the eigenvalues of P−1/2G~P−1/2⪰0P^-1/2 GP^-1/2 0. Then det(Gλ+P)=det(P)∏k(1+gkκλ2) (G_λ+P)= (P) _k(1+g_k _λ^2), so that −12logdet(Gλ+P)=−12∑klog(1+gkκλ2)+const,neff=∑kgkκλ21+gkκλ2.- 12 (G_λ+P)=- 12 _k \! (1+g_k _λ^2 )+const, n_eff= _k g_k _λ^21+g_k _λ^2. (30) Each summand lies in [0,1)[0,1) and is nonzero only for gk>0g_k>0, so neff<rankGλ≤min(Nd,dimθ)n_eff<rankG_λ≤ (Nd, θ) whenever Gλ≠0G_λ≠ 0. Once the physics dominates the prior on every direction it constrains (gkκλ2≫1g_k _λ^2 1), neffn_eff saturates at rankGλrankG_λ and no longer depends on κλ _λ; we call this the saturated regime. For P=0P=0 and full-rank GλG_λ (which requires dimθ≤Nd θ≤ Nd), neff=dimθn_eff= θ for every κλ _λ. Mechanism 1 (residual metric). Since ZrZ_r does not depend on κλ _λ, Eq. (29) gives ∂Δ/∂logκλ=−neff∂ /∂ _λ=-n_eff. With ∑nlog|detAλ(z(tn))|=Ndlogκλ+const _n | A_λ(z(t_n))|=Nd _λ+const, this integrates, for constant neffn_eff, to the metric term −neffNd∑nlog|detAλ(z(tn))|- n_effNd _n | A_λ(z(t_n))| of (16). For a state-dependent Aλ(z(t))A_λ(z(t)), the term holds locally: on time windows that are short compared with the variation of Aλ(z(t))A_λ(z(t)) but long compared with the resolution of the surrogate, AλA_λ is nearly constant, and if the constrained directions are spread uniformly over time and components, each window contributes its share neff/Ndn_eff/Nd of the local sum of log|detAλ| | A_λ|. This neglects the coupling between windows and the anisotropy of AλA_λ relative to Σr _r, both absent in the scalar case. State-space latent. On the Markovian latent of Section 4, θ consists of the TdTd velocities under a flat prior (P=0P=0), and ∂r/∂v∂ r/∂ v is block lower-triangular with invertible diagonal blocks Aλ(zi−1)/hA_λ(z_i-1)/h. Hence neff=Td=Ndn_eff=Td=Nd exactly, and the metric term becomes the exact collider bias (22). Step 4: Drift operator (Mechanism 2). We set Aλ=IdA_λ=I_d, so that λ enters (27) only through λL_λ. We assume that the span of the surrogate’s tangent map is approximately invariant under LλL_λ, with a λ-independent Gram determinant, and that the physics dominates P on the resolved modes. Then 12logdet(Gλ+P) 12 (G_λ+P) equals log|detLλ| | L_λ| restricted to the resolved modes, up to λ-independent terms. On the position component, LλL_λ acts as d2dt2−∂aλ∂vdt−∂aλ∂x d^2dt^2- ∂ a_λ∂ v ddt- ∂ a_λ∂ x, whose matrix symbol at frequency ω has, by the Schur complement, the same determinant as L^λ(ω):=iωI2d−Dfλ L_λ(ω):=iω I_2d-Df_λ. For constant DfλDf_λ (relaxed below), a surrogate with effective bandwidth Ω on [0,][0,T] resolves real Fourier modes at density /πT/π, so that 12logdet(Gλ+P)≈π∫0Ωlog|detL^λ(ω)|ω. 12 (G_λ+P)\;≈\; Tπ _0 | L_λ(ω) |dω. (31) The eigenvalue integral. Let skk=12d\s_k\_k=1^2d be the eigenvalues of DfλDf_λ, so that log|detL^λ(ω)|=∑klog|iω−sk| | L_λ(ω)|= _k |iω-s_k|. The integral (31) diverges as Ω→∞ →∞, but at a rate independent of λ; differentiating with respect to the eigenvalues isolates the λ-dependent part. For a complex pair s=−ρ±iβs=-ρ± iβ with β>0β>0 (ρ>0ρ>0 damped, ρ<0ρ<0 growing), I(ρ,β):=∫0Ωlog|(iω−s)(iω−s¯)|ωI(ρ,β):= _0 |(iω-s)(iω- s)|\,dω satisfies ∂I∂ρ ∂ I∂ρ =∫0Ω(ρ(ω−β)2+ρ2+ρ(ω+β)2+ρ2)ω=arctanΩ−βρ+arctanΩ+βρ=πsign(ρ)+(1/Ω), = _0 ( ρ(ω-β)^2+ρ^2+ ρ(ω+β)^2+ρ^2 )dω= -βρ+ +βρ\;=\;π\,sign(ρ)+O(1/ ), (32) ∂I∂β ∂ I∂β =12log(Ω+β)2+ρ2(Ω−β)2+ρ2=(β/Ω), = 12 ( +β)^2+ρ^2( -β)^2+ρ^2=O(β/ ), (33) where the boundary terms at ω=0ω=0 cancel. Up to (1/Ω)O(1/ ), the damping rate thus enters with weight π and the frequency β not at all, so each pair contributes π|ρ|π|ρ|. A real eigenvalue s=−ρs=-ρ gives ∂ρ∫0Ωlog|iω+ρ|ω=arctan(Ω/ρ)=π2sign(ρ)+(1/Ω) _ρ _0 |iω+ρ|\,dω= ( /ρ)= π2sign(ρ)+O(1/ ), i.e., a contribution π2|ρ| π2|ρ|. Every eigenvalue therefore contributes π2|Resk| π2|Re\,s_k|, and multiplying by the mode density /πT/π gives 12logdet(Gλ+P)=2∑k=12d|Resk|+const+ℛ,ℛ=(/Ω). 12 (G_λ+P)= T2 _k=1^2d|Re\,s_k|+const+R, =O(T/ ). (34) Without growing modes (Resk≤0Re\,s_k≤ 0 for all k), ∑k|Resk|=−tr(Dfλ)=−∇⋅fλ _k|Re\,s_k|=-tr(Df_λ)=-∇\!·\!f_λ, and inserting (34) into (15) gives Δdiss(λ)=2∇⋅fλ+const+ℛ. _diss(λ)= T2\,∇\!·\!f_λ+const+R. (35) Otherwise, Δ(λ) (λ) is not determined by the divergence alone: growing modes are penalized like damped ones, and a divergence that changes sign along the orbit enters through its absolute value. Appendix A.3.6 discusses both cases for conservative systems. Nonlinear drift. For nonlinear dynamics, we treat Dfλ(z(t))Df_λ(z(t)) as constant on windows that are short compared with T but long compared with the correlation time Ω−1 ^-1 of the surrogate and the relaxation times |Resk|−1|Re\,s_k|^-1. Summing over windows replaces ∇⋅fλT\,∇\!·\!f_λ by ∫0∇⋅fλ(z(t),t)t _0^T∇\!·\!f_λ(z(t),t)\,dt and gives Eq. (17). Its leading term is the Onsager–Machlup correction to the path action [19, 7], obtained here as the fluctuation determinant of a band-limited surrogate. The remainder ℛR collects the (1/Ω)O(1/ ) tails of (32) and (33), the residual floor E[θ∗,λ]E[θ^*,λ] of a finite-capacity fit, and the windowing. Since all of these grow with T at the same rate as the leading term, their relative size is set by Ω compared with the time scales of the system, not by the length of the record. State-space latent. Under the Euler–Maruyama transition (18), viv_i depends only on the noise increments ΔWj:=hξj W_j:= h\, _j up to step i, so the Jacobian of (ΔW1,…,ΔWT)↦(v1,…,vT)( W_1,…, W_T) (v_1,…,v_T) is block lower-triangular with diagonal blocks Aλ(zi−1)−1Σr1/2A_λ(z_i-1)^-1 _r^1/2, which do not involve the drift. This is why the exact collider bias (22) has no dissipation term: the Markovian latent normalizes the path density step by step (the Itô pre-point form), whereas a global surrogate, which imposes no ordering, yields the midpoint (Onsager–Machlup) form, consistent with the Wong–Zakai limit of smooth approximations to Brownian paths [28]. Step 5: Collecting the terms. Adding the metric term of Step 3 and the dissipation term of Step 4 gives Eq. (16); without growing modes, its dissipation term is the divergence integral (17). Of the four simplifications listed in Section 3, three need a caveat in practice. Simplification (i) fails when the data-free integrand of ZλZ_λ is dominated by the prior rather than the physics; the terms dropped from (24) then depend on λ. For (i), when x≡0x≡ 0 solves the ODE, θ∗=0θ^*=0, where the tangent map of a neural surrogate degenerates; we therefore evaluate the curvature at a representative state on the solution manifold (Appendix A.3.7). Since the factors in (27) are rectangular, the additivity (iv) holds only to first order. Eq. (16) thus captures the dominant dependence of Δ on λ and shows where the collider bias is large, but adding −Δ- to the log-posterior of a B-PINN need not recover the chain posterior exactly. Table 2 summarizes both mechanisms. Table 2: Taxonomy of structural biases in collider B-PINNs. Scaling denotes the asymptotic rate of the collider bias Δ=log(pcollider/pchain) = (p_collider/p_chain). Null conditions indicate when Δ is invariant to inferred parameters; sk\s_k\ are eigenvalues of DfλDf_λ. Mechanism Physical Origin Scaling Null Condition Residual Metric State-space map AλA_λ depends on λ Θ(neff) (n_eff) or Θ() (T) AλA_λ invariant to λ Dissipation Phase-space volume contraction Θ() (T) Markovian latent, or ∑k|Resk| _k _k invariant to λ Worked example: the damped linear oscillator. For mx¨+cx˙+kx=σrξm x+c x+kx= _rξ with force residuals (Aλ=mA_λ=m, d=1d=1), the drift fλ(x,v)=(v,−kmx−cmv)⊤f_λ(x,v)=(v,\,- kmx- cmv) has the constant Jacobian Dfλ=[01−k/m−c/m]Df_λ= bmatrix0&1\\ -k/m&-c/m bmatrix with divergence ∇⋅fλ=−c/m∇\!·\!f_λ=-c/m. Its eigenvalues are −ρ±iβ-ρ± iβ with ρ=c/(2m)ρ=c/(2m) and β=k/m−c2/(4m2)β= k/m-c^2/(4m^2) if underdamped, and real and negative if overdamped; in both cases ∑k|Resk|=c/m _k|Re\,s_k|=c/m. No windowing is needed, and (16) reads Δ(λ)≈−nefflogm⏟metric−c2m⏟dissipation+const. (λ)≈ -\,n_eff m_metric\; -\; c\,T2m_dissipation\;+\;const. (36) The dissipation term penalizes large c/mc/m in proportion to T, so the collider underestimates damping. The stiffness enters only through β, which by (33) contributes no leading-order term; the (β/Ω)O(β/ ) remainder, however, still grows with T and rewards smaller k/mk/m, so restoring only the leading-order terms leaves part of the stiffness bias (Appendix A.3.7). Along the drift-invariant rescaling (m,c,k)→α(m,c,k)(m,c,k)→α(m,c,k), only the metric term varies, ∂Δ/∂logα=−neff∂ /∂ α=-n_eff: the collider is rewarded without bound for shrinking the mass scale, opposed only by the prior and by the process-noise level, which identifies σr/m _r/m. On the state-space latent, the same holds exactly with neff=N=Tn_eff=N=T, since (22) gives ΔMarkov=−Tlogm+const _Markov=-T m+const; Appendix A.2 discusses its effect on the collider posterior. A.2 Exact Collider Bias along Drift-Invariant Directions While the decomposition derived in Appendix A.1 is generally approximate, the metric bias (Mechanism 1) admits exact statements along drift-invariant directions. Consider a drift-invariant direction parameterized by jointly scaling masses and forces by a factor α. Suppose the scaled components follow independent log-normal priors with log-scale standard deviations sjs_j. Measuring logα α relative to the prior medians λ~j λ_j via the precision-weighted mean logα=∑jsj−2log(λj/λ~j)∑jsj−2, α= _js_j^-2 ( _j/ λ_j) _js_j^-2, (37) the prior factorizes into independent components along logα α and the identified directions, with a Gaussian marginal of standard deviation (∑jsj−2)−1/2( _js_j^-2)^-1/2 along logα α. Under acceleration residuals, the residuals themselves remain invariant under this scaling. Consequently, neither the residual energy (8) nor the observational likelihood depends on α, and ZλZ_λ is independent of α. Both the chain and the collider joint then depend on α solely through the prior, so both posteriors leave the distribution of logα α exactly at its prior, for any latent representation. Under force residuals, Aλ=MλA_λ=M_λ scales with α while Σr _r is fixed, so the process noise on the velocities scales with α−1α^-1: the data constrain α through the noise level, and the chain posterior of logα α differs from its prior. On top of this, the collider joint acquires an additional explicit α-dependence: on the state-space latent, ΔMarkov=−Tdlogα+const _Markov=-Td α+const exactly by Eq. (22), and on a global surrogate, Δ≈−nefflogα+const ≈-n_eff α+const by the metric term of Eq. (16). The collider posterior of logα α is therefore the chain posterior reweighted by α−Tdα^-Td (respectively α−neffα^-n_eff). This reward grows with the number of time steps (respectively constrained directions), whereas the constraint from the noise level is set by the finitely many observations; on the oscillator, where the exact Kalman likelihood is available, it places the state-space collider at logα=−42.4 α=-42.4 (Table 1). A.3 Experimental Details and Full Results This appendix supports Section 5 and follows its order. Appendix A.3.1 specifies the setup, and Appendix A.3.2 verifies that our method samples the chain posterior correctly, since every bias in this paper is measured against it. Appendix A.3.3 reports all parameters of the runs summarized in Table 1, and Appendix A.3.4 shows that the bias on the oscillator is reproducible across datasets, and Appendix A.3.5 how it depends on the residual scale. The next two subsections support the discussion of individual mechanisms in Section 5.2: the dissipation term on the two conservative systems (Appendix A.3.6) and the closed-form correction of the neural surrogate (Appendix A.3.7). Appendix A.3.8 details the scaling-law experiments of Figs. 3 and 4. A.3.1 Setup Evaluation grid. Each of the three systems is combined with two residual metrics (acceleration and force/torque residuals) and with the inference variants listed below. In every configuration, the process noise is isotropic, Σr=σr2Id _r= _r^2I_d, with σr _r fixed at its true value. Within a system, all runs share the dataset, ground truth and prior, the time discretization, number of particles, sampler budget and random seed, and the network and collocation grid of the surrogate. Two runs that differ in a single factor of the grid therefore isolate the effect of that factor, for instance of the residual metric. Systems and data. Table 3 lists the settings that differ between systems. The data are simulated with the same semi-implicit Euler–Maruyama scheme (18) that the inference uses, but on a finer time grid: 2×2× finer for the oscillator and about 50×50× finer for the pendulum and the relativistic particle (Appendix A.3.2 shows that the results do not rely on this match). The number of Euler–Maruyama steps T used for inference and the number of particles were chosen per system by a convergence study, and are identical across all state-space runs of that system, so that chain and collider remain comparable. Table 3: Per-system settings. The MLP surrogate maps a Fourier-feature embedding of dimension 3333 through two hidden layers of width 2424 to the d output coordinates. NdNd is the number of residual components, i.e., collocation points times dimension. Oscillator Double pendulum Relativistic particle ground truth λ (m,c,k)=(1, 1.2, 4)(m,c,k)=(1,\,1.2,\,4) (m1,m2,ℓ1,ℓ2)=(1, 0.5, 1, 0.8)(m_1,m_2, _1, _2)=(1,\,0.5,\,1,\,0.8) (m0,q,c)=(1, 1, 1)(m_0,q,c)=(1,\,1,\,1) residual scale σr _r (fixed) 0.300.30 0.020.02 0.020.02 horizon T 2424 1212 1212 observations K 4040 2020 2020 observation noise σobs _obs 0.020.02 0.020.02 rad 0.0050.005 Euler–Maruyama steps T 12481248 12351235 1216012160 particles 10001000 10001000 20002000 MLP layer widths 33-24-24-1 33-24-24-2 33-24-24-1 dimθ θ 14411441 14661466 14411441 residual components NdNd 128128 720720 300300 Inference variants. • Chain (state space) is our method (Section 4). Its latent is the Euler–Maruyama path z0:Tz_0:T with the normalized transition kernel (19); a bootstrap particle filter estimates the marginal likelihood p^(y|λ) p(y\,|\,λ), and PMMH samples λ. • Collider (state space) uses the same latent, particle filter and sampler, but omits ZλZ_λ. It is a control that differs from our method only in the normalizer. • Collider B-PINN is the standard formulation of Section 2.2: a Fourier-feature MLP surrogate xθx_θ with a Gaussian weight prior, sampled with NUTS jointly over (θ,λ)(θ,λ). • Chain (Fourier MLP) is the Collider B-PINN with the closed-form correction −Δ- for the missing normalizer added to its log-density (Appendix A.3.7). It is run on the oscillator only (Section 5.1). Because the oscillator is linear, the marginal likelihood of its state-space latent is also available exactly from a Kalman filter. We use this in one place. Under force residuals, PMMH cannot reach the posterior of the state-space collider (see below), so on the oscillator we sample it with an adaptive random-walk Metropolis sampler instead, targeting the exact chain likelihood multiplied by eΔ(λ)=m−Te (λ)=m^-T from Eq. (22) (44 chains, R^≤1.001 R≤ 1.001). These entries are marked with ‡ . Applied to the chain model, the same sampler reproduces the PMMH posterior (e.g., logm m: 0.014±0.1390.014± 0.139 against 0.006±0.1380.006± 0.138). Priors. All physical parameters have log-normal priors with median at the ground truth. The standard deviations of the log-parameters are 0.40.4 for masses, charges and the damping coefficient c of the oscillator, and 0.250.25 for lengths, stiffnesses and the light speed c of the relativistic particle. The Gaussian weight prior of the B-PINN follows LeCun scaling at unit gain: each weight of a layer with ninn_in inputs has variance 1/nin1/n_in, and each bias has variance 11. Samplers. PMMH runs for 24,00024,000 iterations. The first 4,0004,000 adapt the proposal, and the first 30%30\% of the remaining 20,00020,000 are discarded as burn-in, which leaves 14,00014,000 draws. NUTS runs 44 chains, each with 3,0003,000 warm-up iterations and 6,0006,000 draws at a maximum tree depth of 1111, which gives 24,00024,000 draws in total. All samplers use random seed 00. Each run uses a single GPU (NVIDIA Quadro RTX 8000 or NVIDIA L40) on an internal compute cluster. On the oscillator, for instance, PMMH takes 442442 s and the B-PINN 5,8345,834 s (NUTS) for these budgets. Software and licenses. All methods are implemented in JAX [6]; NUTS and the autocorrelation-based effective sample sizes use NumPyro [21, 4]. The nested baselines of Appendix A.4 use the public reimplementation of Bolt [5], which additionally depends on Optax. JAX, NumPyro and Optax are released under the Apache License 2.0; the reimplementation does not state a license. We release our code under the MIT license. Reporting conventions. Unless stated otherwise, tables report the posterior mean divided by the ground truth, so that 1.001.00 means no error. There are two exceptions. Table 1 reports the drift-invariant direction as the precision-weighted log scale factor logα α of Eq. (37). Tables 7 and 9 report displacements in log space, i.e., the posterior mean of the log parameter minus the log of its true value. Since the mean of a logarithm is not the logarithm of the mean, these displacements are not exactly the logarithms of the ratios reported elsewhere; for instance, the damping ratio of the Collider B-PINN under acceleration residuals has a displacement of −0.526-0.526 (Table 9), and exp(−0.526)=0.59 (-0.526)=0.59, whereas Table 1 reports the ratio 0.610.61. Runs without valid posterior samples († ). Two kinds of runs do not yield valid posterior samples. The first is the state-space collider under force residuals. Along the drift-invariant direction, its log-density gains −Tdlogα-Td α, which rewards small masses. Together with the log-normal prior, this still defines a proper posterior, but its mode lies far outside the range the sampler can represent; for the oscillator, where the exact Kalman likelihood locates it, at logα≈−42 α≈-42. PMMH moves toward it until the Euler–Maruyama recursion becomes numerically unstable; from then on, the particle filter returns non-finite likelihoods for a large fraction of the proposals, and only a handful of the 20,00020,000 post-adaptation proposals are accepted. The second kind are B-PINN runs whose chains have not mixed (R^>1.05 R>1.05). For both kinds, we report the median of the draws, which describes where the sampler stopped rather than the posterior. Marginal plots show such runs as an arrow at the edge of the axis, labelled with the value reached. A.3.2 Validation of the Reference Posterior All biases in this paper are measured against the chain posterior that our method samples. Three checks confirm that this reference is correct: the first tests the likelihood estimator, the second the simulation of the data, and the third the sampler. Particle filter against the exact likelihood. For the linear oscillator, the Kalman filter gives the exact log-likelihood logp(y|λ) p(y\,|\,λ). At three points along the drift-invariant direction, the particle-filter estimate deviates from it by +0.14+0.14, −0.10-0.10 and +0.27+0.27 nats, against Monte Carlo standard deviations of 0.270.27, 0.380.38 and 0.380.38 across five independent filter runs. The likelihood estimator underlying PMMH thus shows no bias at the resolution relevant for our comparisons. Simulation step. Because the data are simulated with the same scheme that the inference uses, the chain model could benefit from a matched discretization. To rule this out, we regenerated ten datasets of the oscillator configuration at a 32×32× finer step and computed the exact (Kalman) chain posterior on each. Its displacement from the truth is −0.01±0.07-0.01± 0.07 in log(c/m) (c/m) and +0.04±0.02+0.04± 0.02 in log(k/m) (k/m) (mean ± standard error over datasets), and its 90%90\% credible intervals contain the truth in 88 or 99 of the 1010 datasets for every identified parameter. A different sampler. On the double pendulum and the relativistic particle, we also sampled both models on the state-space latent with NUTS instead of PMMH. NUTS runs jointly over (z0,ξ,λ)(z_0,ξ,λ), where ξ=(ξ1,…,ξT)ξ=( _1,…, _T) are the standardized noise increments of Eq. (18). The path z0:Tz_0:T is a deterministic function of these variables, and the standard normal prior on ξ does not depend on λ. In this non-centred parameterization, the chain model therefore needs no explicit normalizer, while the collider differs from it by Δ(λ) (λ) in (22). The resulting posteriors agree with those obtained with PMMH, and the collapse of the collider under torque residuals is reproduced; neither result is therefore an artefact of the particle filter. NUTS is thus a valid alternative to PMMH for the chain model, but a more expensive one on fine time grids: per effective sample of the drift-invariant scale, it takes 4.94.9 s against 2.72.7 s for PMMH on the pendulum (1,2351,235 steps), and about 1,8001,800 s against 66 s on the relativistic particle (12,16012,160 steps; 5151 h of NUTS in total). A.3.3 Full Posterior Results Tables 4–6 report the posterior of every parameter for every inference variant, including the parameters and parameter combinations that Table 1 omits. They show the same pattern as Table 1. The chain model keeps the drift-invariant scale near the truth and recovers every identified quantity to within 1.31.3 posterior standard deviations. The errors of the collider, in contrast, concentrate where Eq. (16) predicts: under force residuals on the parameters that enter the residual metric, and on the parameters that enter the dissipation term (the damping of the oscillator, the light speed of the relativistic particle). Table 4: Oscillator. The three inference legs under both residual metrics, on one dataset and one prior, with σr _r held at its true value. Posterior mean ± standard deviation over 14000 draws (state space) and 24000 draws (B-PINN). † : a walk at the Euler map’s stability boundary, not a posterior; its median is quoted; ‡ : chains that have not merged (R^>1.05 R>1.05), likewise a median rather than a posterior. residual leg m (1.00) c (1.20) k (4.00) c/mc/m (1.20) k/mk/m (4.00) force residual Chain (state space) 1.015±0.141.015± 0.14 1.209±0.201.209± 0.20 3.981±0.513.981± 0.51 1.209±0.241.209± 0.24 3.948±0.433.948± 0.43 Collider (state space) 9.3⋅10−69.3·10^-6 ‡ 3.7⋅10−43.7·10^-4 ‡ 0.033870.03387 ‡ 39.7439.74 ‡ 36563656 ‡ Collider (B-PINN) 0.2391±0.0500.2391± 0.050 0.5663±0.130.5663± 0.13 1.832±0.251.832± 0.25 2.462±0.732.462± 0.73 7.935±1.77.935± 1.7 acceleration residual Chain (state space) 1.035±0.221.035± 0.22 1.219±0.281.219± 0.28 4.030±0.794.030± 0.79 1.191±0.201.191± 0.20 3.927±0.413.927± 0.41 Collider (state space) 1.035±0.221.035± 0.22 1.219±0.281.219± 0.28 4.030±0.794.030± 0.79 1.191±0.201.191± 0.20 3.927±0.413.927± 0.41 Collider (B-PINN) 1.314±0.281.314± 0.28 0.9419±0.240.9419± 0.24 4.115±0.814.115± 0.81 0.7311±0.180.7311± 0.18 3.162±0.353.162± 0.35 Table 5: Double pendulum. The three inference legs under both residual metrics, on one dataset and one prior, with σr _r held at its true value. Posterior mean ± standard deviation over 14000 draws (state space) and 24000 draws (B-PINN). † : a walk at the Euler map’s stability boundary, not a posterior; its median is quoted; ‡ : chains that have not merged (R^>1.05 R>1.05), likewise a median rather than a posterior. residual leg m1m_1 (1.00) m2m_2 (0.500) ℓ1 _1 (1.00) ℓ2 _2 (0.800) m1+m2m_1+m_2 (1.50) m2/m1m_2/m_1 (0.500) force residual Chain (state space) 1.078±0.261.078± 0.26 0.5367±0.130.5367± 0.13 1.011±0.0111.011± 0.011 0.7890±0.0100.7890± 0.010 1.615±0.391.615± 0.39 0.4976±0.0120.4976± 0.012 Collider (state space) 0.17130.1713 ‡ 0.025380.02538 ‡ 1.3861.386 ‡ 0.075890.07589 ‡ 0.19670.1967 ‡ 0.14810.1481 ‡ Collider (B-PINN) 0.08769±0.0200.08769± 0.020 0.02269±0.00540.02269± 0.0054 0.1441±0.0190.1441± 0.019 0.2137±0.0340.2137± 0.034 0.1104±0.0220.1104± 0.022 0.2709±0.0850.2709± 0.085 acceleration residual Chain (state space) 1.021±0.301.021± 0.30 0.5062±0.150.5062± 0.15 1.011±0.0121.011± 0.012 0.7926±0.0110.7926± 0.011 1.527±0.451.527± 0.45 0.4958±0.0110.4958± 0.011 Collider (state space) 1.021±0.301.021± 0.30 0.5062±0.150.5062± 0.15 1.011±0.0121.011± 0.012 0.7926±0.0110.7926± 0.011 1.527±0.451.527± 0.45 0.4958±0.0110.4958± 0.011 Collider (B-PINN) 1.0041.004 ‡ 0.47150.4715 ‡ 1.0231.023 ‡ 0.78080.7808 ‡ 1.4761.476 ‡ 0.46990.4699 ‡ Table 6: Relativistic particle. The three inference legs under both residual metrics, on one dataset and one prior, with σr _r held at its true value. Posterior mean ± standard deviation over 14000 draws (state space) and 24000 draws (B-PINN). † : a walk at the Euler map’s stability boundary, not a posterior; its median is quoted; ‡ : chains that have not merged (R^>1.05 R>1.05), likewise a median rather than a posterior. residual leg m0m_0 (1.00) q (1.00) c (1.00) q/m0q/m_0 (1.00) force residual Chain (state space) 0.8233±0.150.8233± 0.15 0.8681±0.150.8681± 0.15 0.8939±0.120.8939± 0.12 1.057±0.0551.057± 0.055 Collider (state space) 2.1⋅10−52.1·10^-5 ‡ 0.003830.00383 ‡ 41334133 ‡ 182.8182.8 ‡ Collider (B-PINN) 0.1630±0.0360.1630± 0.036 0.1485±0.0320.1485± 0.032 1.608±0.271.608± 0.27 0.9127±0.0360.9127± 0.036 acceleration residual Chain (state space) 0.9941±0.290.9941± 0.29 1.049±0.311.049± 0.31 0.8881±0.0870.8881± 0.087 1.056±0.0461.056± 0.046 Collider (state space) 0.9941±0.290.9941± 0.29 1.049±0.311.049± 0.31 0.8881±0.0870.8881± 0.087 1.056±0.0461.056± 0.046 Collider (B-PINN) 1.101±0.331.101± 0.33 0.9872±0.290.9872± 0.29 1.783±0.251.783± 0.25 0.8969±0.0150.8969± 0.015 A.3.4 Replication over Datasets All runs above use a single dataset per system. On one dataset, a posterior can lie away from the truth by chance; a bias means that it does so systematically, on average over datasets. To distinguish the two, we repeat the oscillator configuration of Fig. 2 (acceleration residuals) on independent datasets, each with a new SDE path and new observation noise. We do this in two complementary ways. Sampled posteriors. We rerun the B-PINN (NUTS) and our chain model (state space, PMMH) with the settings of Appendix A.3.1 on ten datasets, including the one of Fig. 2. All B-PINN runs reach R^≤1.001 R≤ 1.001. Exact posteriors. To remove sampling error and convergence issues altogether, we use ten further datasets and replace the neural network by a surrogate on which all posteriors can be computed exactly. The collider is a B-PINN on a linear cosine basis (6565 coefficients with prior standard deviation 0.30.3, Nd=128Nd=128 collocation residuals). Because the residual is linear in the coefficients, they can be integrated out in closed form, and ZλZ_λ can be restored exactly on the same basis, which gives a chain model on the surrogate. The exact Kalman likelihood of the Euler–Maruyama latent gives the state-space chain. We evaluate all three posteriors on a grid over the identified coordinates (log(c/m),log(k/m))( (c/m), (k/m)), with the log-normal prior of Appendix A.3.1; less than 10−310^-3 of each posterior’s mass lies at the boundary of the grid. Table 7: Replication of the oscillator configuration over datasets (acceleration residuals). Top: the neural B-PINN and our chain model, sampled as in the main text, on ten datasets. Bottom: exactly computed posteriors on a linear surrogate, on ten further datasets. displ.: mean over datasets of the posterior mean of the log parameter minus the log truth, ± its standard error. sd: mean posterior standard deviation. 90%90\% cov.: number of datasets whose central 90%90\% credible interval contains the truth. log(c/m) (c/m) log(k/m) (k/m) Model (latent) displ. sd 90%90\% cov. displ. sd 90%90\% cov. Sampled posteriors, ten datasets Chain (state space) +0.02±0.05+0.02± 0.05 0.170.17 10/1010/10 +0.02±0.02+0.02± 0.02 0.110.11 10/1010/10 Collider B-PINN (MLP) −0.42±0.03-0.42± 0.03 0.240.24 /4/10 −0.12±0.03-0.12± 0.03 0.110.11 /7/10 Exact posteriors on a linear surrogate, ten datasets Chain (state space) −0.03±0.05-0.03± 0.05 0.170.17 8/108/10 +0.02±0.04+0.02± 0.04 0.110.11 9/109/10 Chain (linear surrogate) −0.05±0.05-0.05± 0.05 0.180.18 9/109/10 +0.02±0.03+0.02± 0.03 0.090.09 9/109/10 Collider (linear surrogate) −0.46±0.04-0.46± 0.04 0.230.23 /1/10 −0.08±0.03-0.08± 0.03 0.100.10 /7/10 Results. Table 7 confirms that the displacements reported in the main text are a systematic bias. Over ten datasets, the B-PINN underestimates c/mc/m by 0.42±0.030.42± 0.03 log units and k/mk/m by 0.12±0.030.12± 0.03, and it lies below the truth on every dataset in both coordinates. Measured relative to the chain model on the same dataset, which removes the variation between datasets, the displacements are −0.45±0.02-0.45± 0.02 and −0.14±0.02-0.14± 0.02. The B-PINN posterior is not narrower than the chain’s (mean standard deviation 0.240.24 against 0.170.17 in log(c/m) (c/m)), but it is displaced by almost two of its own standard deviations, so its 90%90\% credible intervals contain the truth in only 44 and 77 of the ten datasets, whereas those of the chain model contain it in all ten. The dataset of Fig. 2 lies at the upper end: its displacement in log(c/m) (c/m) (−0.53-0.53) is the third largest of the ten, against a median of −0.39-0.39. The exact posteriors on the linear surrogate show the same picture, with displacements of −0.46±0.04-0.46± 0.04 and −0.08±0.03-0.08± 0.03 and coverage in 11 and 77 of ten datasets. Restoring ZλZ_λ exactly on the same surrogate removes the bias in both coordinates and restores nominal coverage, as does the state-space chain. The bias is therefore reproducible across datasets, and it is caused by the missing normalizer rather than by the particular dataset or by the limited expressiveness of the surrogate. A.3.5 Sensitivity to the Residual Scale All experiments above fix σr _r at the value that generated the data. To test how the bias depends on this value, we repeat the exact-posterior comparison of Appendix A.3.4 on the same ten oscillator datasets (acceleration residuals), varying σr _r in the data and in the models. For σr≤0.1 _r≤ 0.1, we evaluate the posteriors on a finer grid, since they are narrower than its original spacing. Table 8: Sensitivity to the residual scale (oscillator, acceleration residuals, exact posteriors, ten datasets). Displacement of log(c/m) (c/m) from the truth (mean over datasets ± standard error) and number of datasets whose central 90%90\% credible interval contains the truth. Top: σr _r matched to the data. Bottom: σr _r misspecified. σr _r Chain (state space) Chain (linear surrogate) Collider (linear surrogate) data model displ. 90%90\% cov. displ. 90%90\% cov. displ. 90%90\% cov. 0.050.05 0.050.05 −0.01±0.02-0.01± 0.02 9/109/10 −0.03±0.02-0.03± 0.02 9/109/10 −0.07±0.02-0.07± 0.02 /7/10 0.10.1 0.10.1 −0.01±0.03-0.01± 0.03 9/109/10 −0.02±0.03-0.02± 0.03 9/109/10 −0.12±0.03-0.12± 0.03 /7/10 0.30.3 0.30.3 −0.03±0.05-0.03± 0.05 8/108/10 −0.05±0.05-0.05± 0.05 9/109/10 −0.46±0.04-0.46± 0.04 /1/10 1.01.0 1.01.0 −0.04±0.07-0.04± 0.07 8/108/10 −0.16±0.07-0.16± 0.07 9/109/10 −0.74±0.03-0.74± 0.03 /0/10 00 0.30.3 +0.40±0.01+0.40± 0.01 /0/10 +0.37±0.01+0.37± 0.01 /1/10 −0.15±0.01-0.15± 0.01 10/1010/10 0.30.3 0.90.9 +0.97±0.05+0.97± 0.05 /0/10 +0.96±0.05+0.96± 0.05 /0/10 −0.30±0.02-0.30± 0.02 10/1010/10 With σr _r matched to the data, the collider bias is present at every scale and grows with it, from −0.07-0.07 to −0.74-0.74 in log(c/m) (c/m), while both chain models remain close to the truth (Table 8, top). At small σr _r, the displacement is small in absolute terms but still about one posterior standard deviation, and coverage remains below nominal. A misspecified residual scale introduces a different bias, which affects the chain model as well. For noise-free data, or for a residual scale three times too large, the chain overestimates the damping, and its intervals miss the truth (Table 8, bottom). The downward collider bias then partly offsets this error, so that the collider’s wider intervals cover the truth; this reflects two biases of opposite sign cancelling, not a correct model. Ideally, the residual scale is therefore matched to the data. When it is unknown, a natural remedy is to infer σr _r jointly with λ, as noted in Section 6. For the chain model, this is a direct extension of PMMH, since the per-step normalizers (21) remain analytic in σr _r. For the collider, however, the bias Δ then also depends on σr _r, which gives rise to a further bias mechanism; we leave its analysis to future work. A.3.6 The Dissipation Term on the Conservative Systems The double pendulum and the relativistic particle are both conservative, yet the dissipation term of a global surrogate need not vanish for them. We first explain why, and then evaluate the term for each system. Why conservative systems can retain a dissipation term. Suppose that the forces derive from a (possibly time-dependent) potential and that Mλ=∂2ℒ/∂v∂v⊤M_λ=∂^2L/∂ v\,∂ v is the velocity Hessian of the Lagrangian ℒL. Then the flow is Hamiltonian in the canonical coordinates (x,p)(x,p) with p=∂ℒ/∂vp= /∂ v, and divergence-free in these coordinates (Liouville’s theorem). Since the change of variables (x,v)↦(x,p)(x,v) (x,p) has Jacobian determinant detMλ M_λ, ∇⋅fλ(z(t),t)=−dtlogdetMλ(z(t)),∫0∇⋅fλdt=logdetMλ(z(0))detMλ(z()),∇\!·\!f_λ(z(t),t)=- ddt M_λ(z(t)), _0^T∇\!·\!f_λ\,dt= M_λ(z(0)) M_λ(z(T)), (38) a boundary term that does not grow with T. On a global surrogate, however, Δ(λ) (λ) depends on ∑k|Resk| _k|Re\,s_k| (Eq. (34)) rather than on the divergence, and the two can differ in two ways. First, the eigenvalues of a Hamiltonian linearization come in pairs ±s± s: oscillatory pairs (Res=0Re\,s=0) contribute nothing, but hyperbolic pairs (real ±s± s) contribute |s|T|s| although their divergence vanishes, so a chaotic system can retain a λ-dependent term that grows with T. Second, if DfλDf_λ has a single eigenvalue with nonzero real part, as for a one-dimensional system with velocity-dependent inertia, then ∑k|Resk|=|∇⋅fλ| _k|Re\,s_k|=|∇\!·\!f_λ| pointwise, and the windowed dissipation term becomes −12∫0|∇⋅fλ|dt- 12 _0^T|∇\!·\!f_λ|\,dt. This grows linearly in T whenever the divergence changes sign along the orbit, although ∫0∇⋅fλt _0^T∇\!·\!f_λ\,dt stays bounded. The pendulum is an instance of the first case and the relativistic particle of the second. We evaluate both quantities along the noise-free trajectory at the ground truth; this explains why the dissipation term leaves the pendulum results in Table 1 unaffected but biases the relativistic particle. Double pendulum. The time average of ∑k|Resk| _k|Re\,s_k| over [0,][0,T] is 0.5080.508, so the surrogate retains a term of 2⟨∑k|Resk|⟩=3.05 T2 _k|Re\,s_k| =3.05 nats at =12T=12, although the mean divergence is only 1.2⋅10−21.2· 10^-2. This term does not affect the shift logα α along the drift-invariant direction reported in Table 1: scaling (m1,m2)→α(m1,m2)(m_1,m_2)→α(m_1,m_2) leaves the acceleration aλ=Mλ−1Fλa_λ=M_λ^-1F_λ, and hence DfλDf_λ, unchanged, and indeed the time average agrees to six digits at α=0.5α=0.5, 11 and 22. Off the direction, it is negligible as well: it varies by less than 2⋅10−32· 10^-3 over ±0.05± 0.05 in logℓ1 _1, i.e., |∂Δ/∂logℓ1|≈0.04|∂ /∂ _1|≈ 0.04 nats per unit of logℓ1 _1, which is small against a prior with log-standard deviation 0.250.25. Relativistic particle. The drive returns the particle to rest at both ends of the record (γ(0)=γ()=1γ(0)=γ(T)=1), so ∫0∇⋅fλt=0 _0^T∇\!·\!f_λ\,dt=0 to machine precision: the divergence term is not merely bounded but absent. The drift Jacobian Dfλ=[010∇⋅fλ]Df_λ= bmatrix0&1\\ 0&∇· f_λ bmatrix, however, has the eigenvalues 00 and ∇⋅fλ=−3qE(t)v/(m0c2γ)∇\!·\!f_λ=-3qE(t)v/(m_0c^2γ), which changes sign twice per drive period. Hence ∑k|Resk|=|∇⋅fλ| _k|Re\,s_k|=|∇\!·\!f_λ|, and ∫0|∇⋅fλ|t=10.36 _0^T|∇\!·\!f_λ|\,dt=10.36, so that Δdiss=−5.18 _diss=-5.18 nats at the truth, a term that grows linearly in T. Differentiating with respect to the light speed at a fixed trajectory gives ∂Δdiss/∂logc=+9.5∂ _diss/∂ c=+9.5 nats per unit of logc c. The collider is thus rewarded for a larger c, because a larger c keeps γ closer to one and thereby flattens the divergence. This reward has the sign and the order of magnitude of the observed error, and it is too strong for a prior with log-standard deviation 0.250.25 to counteract. The B-PINN reports c=1.78c=1.78 (logc c displaced by +0.569+0.569), against 0.890.89 for the chain. We do not derive a quantitative prediction from these numbers, because Δdiss _diss is strongly nonlinear in c (roughly ∝−c−2 -c^-2) over the range that the posterior covers. A.3.7 Restoring the Normalizer on the Neural Surrogate In principle, the bias could stem from the probabilistic model (collider or chain), the latent (state space or neural surrogate), or the estimator (PMMH or NUTS). On the state-space latent, the model alone decides: chain and collider share latent and estimator and coincide exactly whenever Δ(λ) (λ) is constant (Section 5.2), and NUTS reproduces the PMMH results (Appendix A.3.2). Here, we hold the neural surrogate fixed and change only the model. Construction. Chain (Fourier MLP) adds the term −Δ=logZr−logZλ- = Z_r- Z_λ to the log-density of the oscillator B-PINN. By Eq. (13), this turns the collider into the chain, up to the approximations of Appendix A.1. For the oscillator, Eq. (16) at d=1d=1, with detAλ=m A_λ=m under force residuals and detAλ=1 A_λ=1 under acceleration residuals, gives logZr−logZλ=c2m⏟dissipation+nefflogm⏟metric, force residuals only+const. Z_r- Z_λ= c\,T2m_dissipation\;+\; n_eff m_metric, force residuals only+const. (39) Everything else (architecture, weight prior, data, collocation grid, sampler and seed) is identical to the B-PINN. Effective number of constraints. The correction requires neff=tr[Gλ(Gλ+P)−1]n_eff=tr[G_λ(G_λ+P)^-1]. We compute it once, from the Gauss–Newton curvature at λ=λtrueλ= _true, and keep it fixed during sampling. This is the idealization under which the closed forms are derived (assumption (i) in Appendix A.1): by Eq. (30), neffn_eff depends on λ only through κλ _λ and is constant in the saturated regime. We evaluate the curvature at the final network weights of the corresponding B-PINN run rather than at the data-free mode θ∗=0θ^*=0, where xθ≡0x_θ≡ 0 solves the ODE but the network Jacobian is degenerate (assumption (i)). With Nd=128Nd=128 and dimθ=1441 θ=1441, this gives neff=110.2n_eff=110.2 for force residuals. For acceleration residuals, only the dissipation term is active and no neffn_eff is needed. This value depends only moderately on λ but strongly on the network weights. Evaluating the curvature at the chain’s posterior mean instead of λtrue _true changes it by less than 0.5%0.5\%, and over draws from the chain posterior it varies by about ±10%± 10\% (9898–116116 around 110.2110.2). At the final state of the acceleration-residual B-PINN, by contrast, it is 99.999.9, and at a MAP fit of the network it drops to 4444. Moreover, neffn_eff is not constant along the B-PINN’s own posterior (5555–7878, where the B-PINN has collapsed the mass scale), so the saturated regime of assumption (i) holds only approximately. This sensitivity to the expansion point is the main source of the overcorrection discussed below. Table 9: Restoring the normalizer on the neural surrogate (oscillator). Displacements from the truth in log space (posterior mean of the log parameter minus the log truth) for the Collider B-PINN and for Chain (Fourier MLP), and the fraction of the displacement removed by the correction, 1−|chain|/|collider|1-|chain|/|collider|. Active terms are the terms of Eq. (39) that depend on the inferred parameters in that configuration. The last configuration is the =24T=24 run of the horizon study in Fig. 3. Configuration Active terms Direction Collider Chain (MLP) Removed force residuals metric, dissipation logm m −1.453-1.453 +0.497+0.497 65.8%65.8\% acceleration residuals dissipation log(c/m) (c/m) −0.526-0.526 −0.013-0.013 97.5%97.5\% Results. Table 9 reports two configurations, each of which activates a different subset of the terms in (39). All runs converge (R^=1.000 R=1.000 over 24,00024,000 draws), so these are valid posteriors. How much of the bias the correction removes follows how accurate the closed forms are in each configuration: • With acceleration residuals, only the dissipation term is active and no neffn_eff enters. The correction removes 97.5%97.5\% of the displacement of log(c/m) (c/m). • With force residuals, the approximate neffn_eff matters, and the correction can overshoot. The collider’s logm=−1.45 m=-1.45 becomes +0.50+0.50 (m≈1.64m≈ 1.64 against a truth of 1.001.00): the correction crosses the truth and leaves a displacement of about one third of the original on the other side. A term linear in logm m with a fixed coefficient neffn_eff replaces a function of κλ _λ and can therefore overcorrect. As cautioned at the end of Appendix A.1, the closed forms thus locate the bias and remove most of it, but they do not calibrate it exactly. Remaining stiffness bias. The correction also leaves part of the stiffness bias: with acceleration residuals, Chain (Fourier MLP) reports k/m=3.48k/m=3.48, against 3.163.16 for the B-PINN and 3.933.93 for the chain on the state-space latent (truth 4.004.00). This remainder is still collider bias rather than a representation error of the surrogate. On the linear cosine surrogate of Appendix A.3.4, where ZλZ_λ can be restored exactly, the collider displaces log(k/m) (k/m) from the truth by −0.08±0.03-0.08± 0.03 (mean ± standard error over ten datasets), whereas the exactly restored chain (+0.02±0.03+0.02± 0.03) and the state-space chain (+0.02±0.04+0.02± 0.04) are unbiased. The leading-order correction misses this part because k enters logZλ Z_λ only through the (β/Ω)O(β/ ) remainder of Eq. (33). A.3.8 Details on the Scaling-Law Experiments Figs. 3 and 4 use the oscillator configuration of Appendix A.3.1, in each case with the mechanism that is not under test switched off. Fig. 4 checks the metric term of Δ(λ) (λ) in a setting where its closed form is exact; Fig. 3 tests the approximate dissipation term on the neural surrogate. Figure 4: Exact closed form of the residual metric bias against measurements on the oscillator (Eq. (22)). Closed form in grey, the collider on the Markovian state-space latent in dark red, the chain in blue. On the state-space latent, Δ(λ) (λ) along the drift-invariant direction matches −Tlogα-T α to within 2.8⋅10−122.8· 10^-12 nats for each number of time steps shown (T=312T=312, 1,2481,248, 2,4962,496; T=624T=624, omitted for clarity, agrees equally), whereas the chain’s profile stays flat. Residual metric bias (Fig. 4). The figure shows Δ(λ) (λ) along the drift-invariant direction (m,c,k)→α(m,c,k)(m,c,k)→α(m,c,k) on the state-space latent, for force residuals; the Markovian latent switches off the dissipation term. We use four discretizations, T∈312,624,1248,2496T∈\312,624,1248,2496\. The measured Δ(λ) (λ) matches −Tlogα-T α to within 2.8⋅10−122.8· 10^-12 nats, i.e., to machine precision, as expected since Eq. (22) is exact on this latent. Over the swept range of α, Δ(λ) (λ) varies by 225225, 449449, 899899 and 17971797 nats, doubling with the number of steps, whereas the chain’s log-posterior varies by only 6.06.0 nats over the same range. Dissipation bias (Fig. 3). Unlike the metric term, the dissipation term is not a property of the latent alone, so we measure its effect on the posterior rather than Δ(λ) (λ) itself. We use acceleration residuals, which switch off the metric bias, and vary the horizon over ∈12,24,48T∈\12,24,48\; the run at =24T=24 is the standard oscillator configuration. Across horizons, we keep the surrogate’s bandwidth per unit time and the Euler step fixed and scale K, N and T with T, so that only the amount of data changes. The panel shows the posterior mean of log(c/m) (c/m) minus its true value for Chain (state space), Chain (Fourier MLP) and the Collider B-PINN, each with a bar of ±1± 1 posterior standard deviation. The state-space collider is not plotted separately, since under acceleration residuals Δ(λ) (λ) is constant and its draws are identical to the chain’s. The B-PINN lies at −0.377-0.377, −0.526-0.526 and −0.669-0.669, i.e., at 1.551.55, 2.082.08 and 2.672.67 of its own posterior standard deviations, which barely change (0.2430.243, 0.2530.253, 0.2510.251). The chain stays near the truth, both on the state-space latent (−0.086-0.086, −0.022-0.022, +0.041+0.041, with standard deviations 0.1980.198, 0.1760.176, 0.1430.143) and on the surrogate with the correction (−0.095-0.095, −0.013-0.013, +0.141+0.141). The chain’s posterior thus contracts as the record grows, whereas the B-PINN’s does not, so the B-PINN’s error grows relative to its own stated uncertainty. To compare with the closed form (36), we reweight the chain posterior of log(c/m) (c/m), approximated as Gaussian, by eΔdiss=e−(c/m)/2e _diss=e^-(c/m)T/2. This predicts shifts relative to the chain of −0.21-0.21, −0.32-0.32 and −0.41-0.41, whereas the B-PINN is shifted by −0.29-0.29, −0.50-0.50 and −0.71-0.71. The leading-order term thus predicts the direction of the bias and its growth with the record, but only between three quarters and three fifths of its size; it also predicts a posterior that contracts with T, which the B-PINN’s does not. A.4 Comparison with Nested Inference To our knowledge, the only prior work that adopts the chain model, by Alberts, Hao, and Bilionis [1, 10, 11], handles the unknown normalizer ZλZ_λ by nested inference: an inner loop estimates the gradient ∂λlogZλ _λ Z_λ anew at every step of an outer loop over the joint state. We compare our method, the particle filter with PMMH (Section 4.3), against the two nested schemes that target the same chain posterior. • NSVI [10] (nested stochastic variational inference) fits an outer Gaussian guide over path coefficients, physical parameters and initial state by maximizing an ELBO. The normalizer gradient in this ELBO is supplied by a second, persistent inner guide, which is re-fitted at every outer step. • NPSGLD [11] (nested preconditioned stochastic-gradient Langevin dynamics) replaces both guides by preconditioned Langevin chains: a persistent auxiliary chain over the path supplies the normalizer gradient to a main chain over the joint state. We use the public JAX reimplementation of both methods by Bolt [5] (commit dab4403), which is not by the original authors; we run it unmodified and supply only the problem definition and the evaluation. To isolate the inference algorithm, all methods use the Euler–Maruyama path as latent, on which ZλZ_λ is in fact available in closed form. The comparison therefore measures how expensive and how accurate nested estimation of ZλZ_λ is where it could be avoided. The baselines were designed for latents without a closed-form normalizer, a setting we do not test. Setup. We use two systems in the configurations of Section 5.1: the damped oscillator under both residual metrics, and the double pendulum under torque residuals. All methods share the model, the latent (T=1248T=1248 Euler–Maruyama steps for the oscillator and T=1235T=1235 for the pendulum), the observations, and the log-normal priors with median at the truth. As in all experiments, σr _r is fixed at its true value. Instead of tuning a single run per baseline, we sweep their hyperparameters. For NSVI, we vary the outer learning rate in [5⋅10−5,10−3][5· 10^-5,10^-3] and the number of inner iterations from 88 to 256256, with budgets of up to 1.2⋅1061.2· 10^6 iterations. For NPSGLD, we vary the step size in [10−8,10−3][10^-8,10^-3], use 44 or 1616 chains and both of its preconditioners (rmsprop and dense_fisher), and switch step-size annealing on or off, with budgets of up to 10710^7 iterations. This amounts to 2323 runs on the oscillator and 2121 on the pendulum. All runs, including those of our method, use the same Quadro RTX 8000 GPU, so that run times are directly comparable. Evaluation metrics. Wall-clock time alone favours samplers that have not mixed, since they produce many draws that carry little information. We therefore measure cost in seconds per effective sample, i.e., wall-clock time divided by the effective sample size (ESS), and report the split-R R alongside it. As in Table 1, runs with R^>1.05 R>1.05 do not yield a valid posterior; we report their median and mark them with a dagger. For our method, the ESS is the smallest over the components of λ; for NPSGLD, it is the ESS of the drift-invariant scale. This choice favours the baseline: on the drift-invariant scale alone, the ESS of our method on the pendulum is about 870870 rather than 343343, which would lower its cost from 6.716.71 to 2.62.6 s per effective sample. The cost and R R of our method refer to a single PMMH chain. To check its mixing, we additionally ran four independent chains (sampler seeds 00–33) per configuration; they give split-R^≤1.001 R≤ 1.001 for all parameters, and their posterior means of the drift-invariant scale agree to within 0.020.02 (oscillator) and 0.0030.003 (pendulum). NSVI draws independent samples from a fitted guide, so R R and ESS do not apply to it; we assess it by its bias instead. Finally, the width of a run is the ratio of its posterior standard deviation to that of our method. It measures what matters for uncertainty quantification: a method can locate the posterior correctly and still misrepresent its uncertainty. Table 10: Nested inference compared with our particle filter with PMMH. Ranges are over the hyperparameter sweep of each method. Location is the drift-invariant scale (m for the oscillator, m1+m2m_1+m_2 for the pendulum), reported as posterior mean ± standard deviation where R^≤1.05 R≤ 1.05 and as a daggered median otherwise. Width is the posterior standard deviation divided by that of our method; 1.001.00 is correct, and 0.030.03 means thirty-fold overconfidence. Cost is in seconds per effective sample. A run counts as diverged if its location is off by more than a factor of 100100 or its spread exceeds 100100 reference standard deviations. System Method Runs Location R R ESS Width Cost Oscillator PF + PMMH (ours) 11 1.015±0.1401.015± 0.140 1.001.00 604604 1.001.00 0.810.81 force residuals NSVI 88 0.780.78–1.431.43 — — 0.150.15–0.280.28 — NPSGLD 1010 1.001.00–1.121.12 † 1.31.3–3.03.0 33–3232 0.030.03–1.251.25 5656–23852385 Oscillator PF + PMMH (ours) 11 1.035±0.2201.035± 0.220 1.001.00 875875 1.001.00 0.510.51 accel. residuals NSVI 22 1.281.28 — — 1.361.36–1.381.38 — NPSGLD 33 1.001.00–1.121.12 † 1.11.1–1.51.5 66–8686 1.091.09–1.901.90 2121–190190 Double pendulum PF + PMMH (ours) 11 1.615±0.3941.615± 0.394 1.001.00 343343 1.001.00 6.716.71 torque residuals NSVI 1616 0.170.17–1.581.58 — — 0.030.03–0.100.10 — NPSGLD 55 1.501.50–1.921.92 † , 22 diverged 2.52.5–277277 33–44 0.070.07–0.940.94 141141–495495 Cost. Our method produces one effective sample every 0.50.5–6.76.7 s. NPSGLD, the only baseline that produces MCMC samples, needs 2121 to 2,3852,385 s per effective sample, i.e., 21×21× to 2,900×2,900× more, depending on the configuration. Its cheapest pendulum run finishes about five times sooner than our method, but it is 21×21× more expensive per effective sample, since its 423423 s yield only three effective samples. A larger budget does not improve its mixing either: the longest NPSGLD run, with 10710^7 iterations over 11,44411,444 s, yields five effective samples. NSVI converges to a biased posterior. On the oscillator, NSVI is fast, and most of its settings agree on the drift-invariant scale; but the posterior they agree on is the wrong one. In six of its eight force-residual runs, the damping ratio lies at c/m=0.78c/m=0.78–0.810.81, whereas the truth is 1.201.20 and our method gives 1.209±0.2361.209± 0.236, an underestimate of about 1.81.8 reference standard deviations. The other two runs, both with 10510^5 outer iterations, end at c/m=0.69c/m=0.69 and 0.940.94. More computation does not reduce the bias: raising the outer iterations from 3⋅1043· 10^4 to 3⋅1053· 10^5 and the inner iterations from 3232 to 256256 increases the run time 113113-fold, from 128128 s to 14,47414,474 s, but moves c/mc/m only from 0.7900.790 to 0.8000.800. The bias is therefore a property of the converged approximation, not of insufficient computation. On the chaotic pendulum, NSVI does not converge at all: across 1616 settings, the mass scale ranges from 0.170.17 to 1.581.58, so its final estimate reflects where the budget ran out rather than the posterior. It is also strongly overconfident there: its posterior standard deviations are only 0.030.03–0.100.10 times that of our method (Table 10), i.e., ten to thirty times too small, so that its scattered estimates come with misleadingly narrow uncertainty. NPSGLD misrepresents the uncertainty. None of the 1818 NPSGLD runs, on either system and in any configuration, reaches R^≤1.05 R≤ 1.05. Because its chains do not mix, its narrow posteriors reflect insufficient exploration rather than precision. On the pendulum, NPSGLD reports m1+m2=1.495m_1+m_2=1.495 with a standard deviation fourteen times too small; on the oscillator with force residuals, it reports m=1.000m=1.000 with a standard deviation thirty times too small. At the same time, its posteriors over the identified ratios are up to 2.5×2.5× too wide. The NPSGLD posterior is therefore not only mislocated but also misshapen. Our method, in contrast, converges in every configuration (R^≤1.001 R≤ 1.001 over four independent chains) and is asymptotically exact. Why the two systems differ. Both baselines perform considerably better on the oscillator than on the pendulum, for structural reasons. The oscillator is linear, so its path posterior is close to Gaussian, and the Gaussian outer guide of NSVI is nearly well specified. Moreover, with force residuals, the comparatively large process noise (σr=0.30 _r=0.30) determines the drift-invariant scale, which is the direction that nested schemes find most difficult. The pendulum, in contrast, is chaotic, its mass matrix depends on the state, and its path posterior over 1,2351,235 steps and two degrees of freedom is strongly non-Gaussian. Nested inference therefore does not fail everywhere, but it succeeds only where the problem is close to what its approximations assume, and even there NPSGLD costs more than 40×40× as much per effective sample as our method.