arXiv CSopen access

Structured Inference with Large Language Gibbs

Structured Inference with Large Language Gibbs Sanghyeok Choi1

Henry Gouk1

Esmeralda S. Whitammer1,2

1 University of Edinburgh, School of Informatics 2 CIFAR Fellow, Learning in Machines and Brains

{sanghyeok.choi, henry.gouk, esmeralda.whitammer}@ed.ac.uk

arXiv:2606.19264v1 [cs.LG] 17 Jun 2026

Abstract

left-to-right pass, but sequential decoding can introduce order-dependent artifacts such as recency effects (Lu et al., 2022) and ignorance of context (Liu et al., 2024), which may distort the intended joint distribution. In this work, we ask whether using the conditional beliefs of LLMs over individual variables allows extracting more information than autoregressive generation permits. Concretely, we study Gibbs-style sampling procedures that iteratively resample one variable conditioned on the current values of the others (Fig. 1), motivated by recent work on the use of masked language models as Gibbs transition operators (Wang and Cho, 2019; Yamakoshi et al., 2022; Torroba Hennigen and Kim, 2023) and the prior use of Gibbs sampling for simple structured prediction tasks with LLMs (Zhu et al., 2024; Zhu and Griffiths, 2025) and inference of reasoning chains (Xu et al., 2024); see related work in §A. This procedure turns the model into a transition operator over structured states whose stationary distribution reflects a compromise among local conditionals, removing the bias of any fixed generation order. We propose several ways to define the transition operator using the LLM’s next-token conditionals (§3.1 and 3.2) and study their interaction with the chosen model. We demonstrate this large language Gibbs approach on three illustrative applications. First, for synthetic distributions, we show that sequential prompting can fail to produce conditionally independent and marginally unbiased samples, consistent with recent observations (Hopkins et al., 2023; Hu et al., 2024; Gu et al., 2026), while Gibbs updates can substantially correct this behaviour (§4). Second, for unsupervised consistent reasoning on benchmarks such as GSM8K and TruthfulQA, we use iterative conditional updates to find coherent answers across sets of related questions (§5.1) and find that this offers a more principled and better-performing alternative to

The knowledge encoded in large language models (LLMs) can serve as a substrate for structured reasoning over variables describing a complex world, but accessing this knowledge in a probabilistically coherent manner poses a difficult inference problem. We propose Large Language Gibbs, a scheme for structured probabilistic inference that uses conditional distributions of an LLM as transition operators. Rather than sampling structured objects through single-pass autoregressive generation, we iteratively resample individual variables conditioned on others using an LLM’s nexttoken conditionals. This approach avoids orderdependent biases and produces a stationary distribution that reflects a compromise between all local conditionals. We apply this approach to sampling from synthetic distributions, consistent reasoning tasks, and Bayesian structure learning. The results suggest that the use of LLM conditionals in MCMC is a practical alternative to one-pass generation for structured probabilistic inference under a world prior accessible through noisy LLM conditionals. Code: https://github.com/hyeok9855/ large-language-gibbs.

1

Introduction

Large language models (LLMs) are increasingly used to represent distributions over structured objects, such as label assignments over large input sets (Gadetsky et al., 2025; Wen et al., 2025), tabular data records (Borisov et al., 2023), and other combinatorial objects. In these settings, where the target sample space is high-dimensional and the joint distribution is constrained by statistical dependencies among variables, the estimation of the joint distribution matters as much as the quality of marginal predictions on a single variable, which is the focus of most LLM applications. A common approach is to generate the full object in a single 1

Metadata Breed Coat Age (years) Sleep (h/day)

enum(...) enum(...) int[0,30] int[0,24]

LM (autoreg. init.)

Sample 1

Sample 2

Sample 3

Sample 𝑛

shorthair tuxedo 3 12

shorthair calico 3 12

shorthair calico 3 14

abyssinian tabby 10 15

LM

LM

LM

···

LM

resampling prompt + constraints What is the brightest star in the night sky? The North Star is the brightest star in the night sky. Sirius is the brightest star in the night sky. Polaris is the brightest star in the night sky. The sun is the brightest star in the night sky.

bool bool bool bool

LM (autoreg. init.)

Sample 1

Sample 2

True True False True

False True False True

LM

Sample 3 LM

False True False False

LM

···

resampling prompt + constraints

Figure 1: Above: Several iterations of large language Gibbs sampling for four variables jointly describing a cat. The initial values are sampled autoregressively, then updates are performed by resampling one variable at a time from the LLM’s conditional distribution given the others, serialised in random order (§3.1). Metadata describing the variables can be used to impose constraints on the resampling. Below: Application to consistent reasoning, where each variable is a binary answer to a question (§5.1). Resampling encourages consistency among the answers. For example, the LM’s knowledge that ‘Polaris’ and ‘North Star’ refer to the same star may be being imposed in the first resampling step and the knowledge that at most one star is the brightest in the second step.

recent algorithms (Wen et al., 2025). Third, for Bayesian structure learning from tabular data, we use Gibbs sampling to generate synthetic data and define a world-knowledge-driven prior over causal structures, finding that this can improve structure learning performance on real datasets (§5.2).

2

𝑝, and define a Markov chain by the resampling operator 𝑖 ∼ Unif ({1, . . . , 𝑛}), (𝑡 −1) 𝑋𝑖(𝑡 ) ∼ 𝑞 𝑖 (𝑋𝑖(𝑡 ) | 𝑋−𝑖 ).

In contrast to (1), the estimate 𝑞 𝑖 is used for resampling instead of the true conditional of 𝑝. If each 𝑞 𝑖 has full support, then the kernel is still given by a stochastic matrix whose 𝑛-th power is positive, so it has a unique stationary distribution, to which iterative resampling using (2) converges, even if the unary conditionals 𝑞 𝑖 are not compatible with any joint distribution. Computing the density of this stationary distribution, or its marginals, amounts Î to a 𝑖 |X𝑖 |-dimensional eigenvector computation, but sampling from it can be easily approximated simply by simulating the Markov chain. Note also that the stationary distribution of (2) is only one notion among many of a joint distribution that is optimally compatible with a collection of unary conditionals that do not necessarily arise from the same joint. For example, the ArnoldGokhale construction (Arnold and Gokhale, 1998) reconstructs a joint distribution from conditionals 𝑞 𝑖 as the one whose unary conditionals are KLnearest to the 𝑞 𝑖 .

Background: Gibbs Sampling

Consider a random variable 𝑋 = (𝑋1 , . . . , 𝑋𝑛 ) taking values in X = X1 × · · · × X𝑛 , where each X𝑖 is a finite (or perhaps infinite discrete) space. It is well-known (Geman and Geman, 1984; MacKay, 2003) that if the joint probability mass function 𝑝 of 𝑋 has full support, then it is the unique
∞stationary distribution of the Markov chain (𝑡 ) 𝑋 𝑡=1 defined by the Gibbs kernel, whose transition rule is given by 𝑖 ∼ Unif ({1, . . . , 𝑛}), (𝑡 −1) 𝑋𝑖(𝑡 ) ∼ 𝑝(𝑋𝑖(𝑡 ) | 𝑋−𝑖 ),

(2)

(1)

where 𝑋−𝑖 indicates that 𝑖-th variable is omitted from 𝑋. That is, the kernel randomly resamples one of the variables from its conditional given all other variables. Variants may select 𝑖 non-uniformly (for example, introducing temporal dependence by resampling variables in a fixed order) or resample multiple variables at a time, known as block Gibbs sampling. See §3.3. Suppose now that one has a collection of approximations 𝑞 𝑖 (𝑋𝑖 | 𝑋−𝑖 ) to the unary conditionals of

3

Large Language Gibbs

We consider approximate Gibbs kernels (2) in which the 𝑞 𝑖 are defined by an autoregressive 2

Algorithm. Unravelling the approximate Gibbs procedure with the kernel described above, the basic large language Gibbs procedure is the following: (1) Initialise variables 𝑋1 , . . . , 𝑋𝑛 (e.g., using autoregressive generation with 𝑝 LM ). (2) For 𝑡 = 1, . . . , 𝑇: (a) Sample an index 𝑖 uniformly from 1, . . . , 𝑛. (b) Sample a bijection 𝜎−𝑖 uniformly from the set of all bijections {1, . . . , 𝑛 − 1} → {1, . . . , 𝑛} \ {𝑖}. (c) Sample 𝑋𝑖 from 𝑝 LM (𝑋𝑖 | seq(𝑋−𝑖 , 𝜎−𝑖 ) ⊕ [𝑐 𝑖 ]). (d) Yield current value of 𝑋 as 𝑋 (𝑡 ) . After simulating the Gibbs chain, we collect a certain number of samples using typical MCMC procedures, such as burn-in (forgetting a fixed number of samples from the start of the chain) and thinning (outputting every 𝐾-th element of the chain). On the stationary distribution. Because the transition kernels have full support, the large language Gibbs chain has a unique stationary distribution, which we denote by 𝑞 ∗ . As mentioned before, the conditionals defined in (3) are not necessarily compatible with any joint, and in general ∗ 𝑞 LM 𝑖 (𝑋𝑖 | 𝑋 −𝑖 ) ≠ 𝑞 (𝑋𝑖 | 𝑋 −𝑖 ). Below, we construct a natural candidate joint 𝑞 sym and identify a sufficient condition under which 𝑞 ∗ = 𝑞 sym , which also exposes why the random shuffling step in (3) is essential. We start by rewriting (3) as

(large) language model. Assume that we have a fixed language model 𝑝 LM , with next-token distribution 𝑝 LM (𝑤 𝑗 | 𝑤 1: 𝑗 −1 ) for any token sequence 𝑤 1: 𝑗 −1 in a fixed vocabulary. This language model defines a conditional distribution over continuations of any token sequence. We are interested in the joint distribution over variables 𝑋1 , . . . , 𝑋𝑛 that 𝑝 LM models. Each variable 𝑋𝑖 is assumed to be associated with a natural language context 𝑐 𝑖 that describes it, and its sample space X𝑖 is considered to be in bijection with finite-length sequences of tokens under specific constraints (for example, 𝑋𝑖 may be an integer in a certain range, which can be represented by numeric tokens with appropriate constraints, or may be forced to belong to a certain set of strings). We assume that the variables possess no inherent semantic ordering. The main question is how to define a resampling kernel 𝑞 LM 𝑖 (𝑋𝑖 | 𝑋 −𝑖 ) using the language model, which only gives the conditional distribution over next token(s) given a previous token sequence. 3.1

Large Language Gibbs: Basic Version

We can define a kernel as the expectation of autoregressive kernels under random permutations, 𝑞 LM 𝑖 (𝑋𝑖 | 𝑋 −𝑖 )   (3) ≜ E 𝜎−𝑖 𝑝 LM (𝑋𝑖 | seq(𝑋−𝑖 , 𝜎−𝑖 ) ⊕ [𝑐 𝑖 ]) , where 𝜎−𝑖 is drawn uniformly from the set of all bijections {1, . . . , 𝑛 − 1} → {1, . . . , 𝑛} \ {𝑖}, ⊕ denotes sequence concatenation, and seq(𝑋, 𝜎) serialises the variables in 𝑋 into a token sequence [𝑐 𝜎 (1) , 𝑋 𝜎 (1) , . . . , 𝑐 𝜎 (𝑛−1) , 𝑋 𝜎 (𝑛−1) ] in the order 𝜎 (with appropriate seperation tokens), possibly prepended with some global task metadata or an instruction prompt. For example, at the first step of resampling in the example at the top of Fig. 1, the variable index to be resampled is 𝑖 = 2 (“Coat”). With a certain 𝜎−2 , we may resample 𝑋2 by continuing the following prompt:

𝑞 LM 𝑖 (𝑋𝑖 | 𝑋 −𝑖 )   (4) ∝ E 𝜎−𝑖 𝑝 LM (seq(𝑋−𝑖 , 𝜎−𝑖 ) ⊕ [𝑐 𝑖 , 𝑋𝑖 ]) . Next, we define a target joint distribution   𝑞 sym (𝑋) ≜ E 𝜎 𝑝 LM (seq(𝑋, 𝜎)) ,

(5)

where 𝜎 is uniform on the set of all bijections from {1, . . . , 𝑛} to itself. This is a natural choice when no variable ordering is preferred over others, and since each 𝑝 LM (seq(·, 𝜎)) is a probability density, so is 𝑞 sym . A bijection 𝜎−𝑖 in (4) corresponds to a 𝜎 in (5) satisfying 𝜎(𝑛) = 𝑖. Using this correspondence, we define the target with fixed last-position index 𝑖 as   sym 𝑞 𝑖 (𝑋) ≜ E 𝜎 : 𝜎 (𝑛)=𝑖 𝑝 LM (seq(𝑋, 𝜎)) , (6) Í sym so 𝑞 sym (𝑋) = 𝑛1 𝑛𝑗=1 𝑞 𝑗 (𝑋). With this definition, (4) can be rewritten as

Here are some characteristics of a cat. Sleep (h/day, int[0,24]): 12 Breed (enum(...)): shorthair Age (years, int[0,30]): 3 Coat (enum(...)):

Note that the variables are shuffled. Practically, to enforce constraints, we can use “structured output” generation in modern language modelling software, e.g., vLLM (Kwon et al., 2023).

sym

𝑞 LM 𝑖 (𝑋𝑖 | 𝑋 −𝑖 ) ∝ 𝑞 𝑖 3

(𝑋)

(7)

3.2

(proportionality in 𝑋𝑖 ). On the other hand, the unary conditional of 𝑞 sym is 𝑞 sym (𝑋𝑖 | 𝑋−𝑖 ) ∝ Í sym The two systems of condition𝑗 𝑞 𝑗 (𝑋). sym als therefore agree precisely when 𝑞 𝑖 (𝑋) ∝ 𝑞 sym (𝑋) as a function of 𝑋𝑖 , for each 𝑖. Formally, a sufficient condition is the following.

The basic large language Gibbs procedure in §3.1 directly uses the language model’s generative conditional 𝑝 LM (𝑋𝑖 | seq(𝑋−𝑖 , 𝜎−𝑖 ) ⊕ [𝑐 𝑖 ]) as the resampling kernel. Depending on the language model type (e.g., instruction-tuned vs. base model, see §4 for motivating examples), the LM’s discriminative prediction may be more reliable than generation, as reinforcement learning fine-tuning can lead to miscalibrated probabilities (Swamy et al., 2026). This motivates Metropolis-within-Gibbs kernels in which a candidate is drawn from a simple proposal (e.g., uniform) and accepted via the LM’s discriminative output. The algorithm is exactly the same as before (§2) except that we replace the resampling step (2c) with the following: (i) Draw a candidate 𝑋𝑖′ ∼ 𝑞(·) from a proposal 𝑞 on X𝑖 (taken to be uniform throughout). (ii) Accept or reject 𝑋𝑖′ by querying the language model. We consider two ways to implement the acceptance step (ii). Barker Gibbs. Present the language model with the current value 𝑋𝑖 and the candidate value 𝑋𝑖′ , and ask it to choose the one that is more plausible given 𝑋−𝑖 (serialised in a random order). Let
LM ′ 𝑝 𝑋𝑖 ≻ 𝑋𝑖 | 𝑋−𝑖 be the probability of 𝑋𝑖′ being preferred by the language model. If we assume that the language model is calibrated in the sense of Barker’s rule (Barker, 1965), i.e.,

Assumption 3.1. For all 𝑖 ∈ {1, . . . , 𝑛} and all 𝑋, sym

𝑞𝑖

(𝑋) = 𝑞 sym (𝑋).

(8)

This assumption states that the joint likelihood of the variables as estimated by the LM is invariant to the last variable in the order in which they are serialised. Under this assumption, we have the following proposition. Proposition 3.2. Under Assumption 3.1, the conditionals (3) coincide with the unary conditionals of 𝑞 sym : sym 𝑞 LM (𝑋𝑖 | 𝑋−𝑖 ), 𝑖 (𝑋𝑖 | 𝑋 −𝑖 ) = 𝑞

∀ 𝑖, 𝑋.

Consequently, 𝑞 sym is invariant under the large language Gibbs kernel; if in addition each 𝑞 LM 𝑖 (· | 𝑋−𝑖 ) has full support on X𝑖 , the chain is irreducible and 𝑞 ∗ = 𝑞 sym . sym

Proof. By (7), 𝑞 LM 𝑖 (𝑋𝑖 | 𝑋 −𝑖 ) ∝ 𝑞 𝑖 (𝑋) in 𝑋𝑖 . sym Under Assumption 3.1, 𝑞 𝑖 (𝑋) = 𝑞 sym (𝑋), so sym (𝑋) ∝ 𝑞 sym (𝑋 | 𝑋 ) in 𝑞 LM 𝑖 −𝑖 𝑖 (𝑋𝑖 | 𝑋 −𝑖 ) ∝ 𝑞 𝑋𝑖 ; both sides being normalised in 𝑋𝑖 , they are equal. Compatibility of the unary conditionals with 𝑞 sym then makes 𝑞 sym stationary by the standard Gibbs argument; positivity gives irreducibility, hence uniqueness, i.e., 𝑞 ∗ = 𝑞 sym . □

𝑝 LM 𝑋𝑖′ ≻ 𝑋𝑖 | 𝑋−𝑖 𝑞 ∗𝑖 (𝑋−𝑖 , 𝑋𝑖′ ) = ∗ , 𝑞 𝑖 (𝑋−𝑖 , 𝑋𝑖 ) + 𝑞 ∗𝑖 (𝑋−𝑖 , 𝑋𝑖′ )

(9)

where 𝑞 ∗𝑖 is the stationary distribution of the basic large language Gibbs, then the algorithm (combined with the uniform proposal) also has 𝑞 ∗ as the stationary distribution. The usage of Barker’s rule is also proposed for MCMC with LLMs in Zhu et al. (2024), and we extend it by combining it with the Gibbs sampling framework. Gambling Gibbs. Here, we cast the decision of accepting 𝑋𝑖′ or not as a gamble. Along with the random candidate 𝑋𝑖′ , we sample 𝑉 ∼ 𝑈 [0, 100] and ask the language model whether it would place a bet $𝑉 that 𝑋𝑖′ is more plausible than 𝑋𝑖 , which will pay out $100 if correct. If the model places a bet, then we accept 𝑋𝑖′ . Again, if we assume the language model’s internal belief of plausibility of

Now, consider the augmented distribution 1 LM 𝑝 (seq(𝑋, 𝜎)), in which 𝜎 is a 𝑃(𝑋, 𝜎) = 𝑛! priori uniform. Then (8) is equivalent to 𝑃 (𝜎(𝑛) = 𝑖 | 𝑋) = 𝑛1

Kernels with Discriminative Prediction

∀ 𝑖, 𝑋.

A simple sufficient condition for this is full permutation invariance – 𝑝 LM (seq(𝑋, 𝜎)) being constant in 𝜎 – but Assumption 3.1 is strictly weaker, because it only constrains the average over the ordering of 𝑋−𝑖 (6). In practice, while the assumption of permutation invariance or Assumption 3.1 is unlikely to hold exactly, 𝑞 LM removes the effect of 𝑖 ordering biases in 𝑋−𝑖 , which gives a better approximation to 𝑞 sym than any fixed-order conditionals. 4

Figure 2: Empirical distribution of generated samples from Llama-3.1-8B (base) for target distribution as 𝑈 [{0, 1, . . . , 99}] (top) and N (0, 1) (bottom). Results are aggregated over 25 seeds, with 256 samples generated per seed. Maximum autocorrelation over 1-128 lags and over 25 runs is also reported (max |ACF|). Independent sampling has high bias; batch sampling has high autocorrelation. Gibbs substantially corrects both issues.

Figure 3: Empirical distribution of generated samples from Llama-3.1-8B-Instruct, same format as Fig. 2.

𝑋𝑖 is proportional to 𝑞 ∗𝑖 (𝑋−𝑖 , 𝑋𝑖 ), it should place

stationary distribution may be different than with random choices of variables to resample.

𝑞𝑖∗ (𝑋−𝑖 ,𝑋𝑖′ ) 𝑉 the bet only if 100 < 𝑞∗ (𝑋−𝑖 ,𝑋 ∗ ′ , i.e., 𝑖 )+𝑞𝑖 (𝑋−𝑖 ,𝑋𝑖 ) 𝑖

𝑝

Block Gibbs. One can resample sets of 𝐵 > 1 variables at a time from an approximation to their conditional given the others. The large language Gibbs kernel is readily generalised to this case: one serialises the remaining 𝑛 − 𝐵 variables in random order and places them in the context, then autoregressively resamples the 𝐵 variables to be updated. Such resampling can help the chain escape local modes in the case of highly correlated variables. Similar results to Proposition 3.2, characterising the stationary distribution, can be obtained for the block Gibbs case. The corresponding assumptions generalising Assumption 3.1 require that the joint likelihood symmetrised over all permutations equals the joint likelihood symmetrised only over permutations where the sets of the first 𝑛 − 𝐵 and last 𝐵 indices are fixed.

𝑋𝑖′ ≻ 𝑋𝑖 | 𝑋−𝑖

LM

 (10) 𝑞 ∗𝑖 (𝑋−𝑖 , 𝑋𝑖′ ) 𝑉 . =P < ∗ ∗ ′ 100 𝑞 𝑖 (𝑋−𝑖 , 𝑋𝑖 ) + 𝑞 𝑖 (𝑋−𝑖 , 𝑋𝑖 ) 

which recovers the Barker acceptance probability (9). The advantage of Gambling Gibbs over Barker Gibbs is that the language model is required only to provide a deterministic threshold judgment, which can be obtained by greedy decoding and does not rely upon calibrated probabilities. 3.3

Practical Variants

Variations on the basic Gibbs sampling from §2 are possible. We describe two that we consider in combination with the large language Gibbs kernels. Sweeping Gibbs. Instead of randomly selecting a variable to resample at each step, we can cycle through the variables in a fixed order, a common practice in Gibbs sampling that can lead to faster convergence. If the resampling kernels are exact conditionals of some target distribution, then this procedure has that target as its stationary distribution. However, in the case of approximate conditionals that are not compatible with any joint, the

4

Motivating Example: Sampling from Simple Distributions

It has been reported multiple times in the literature that LLMs are not able to sample from simple distributions over integers and real numbers, such as uniform and Gaussian (Hopkins et al., 2023; Hu et al., 2024; Gu et al., 2026). Notably, Gu 5

Table 1: Task accuracies on the consistent reasoning tasks (mean±std over three seeds).

et al. (2026) showed that batched sampling with LLMs (generating a batch of samples sequentially in one shot) gives better marginal distributions than independent sampling (generating one sample at a time). However, the batch samples sometimes have high autocorrelation, meaning that the samples are not independent despite having more accurate marginals. Here, we use large language Gibbs and its variants to sample from simple distributions: the uniform distribution over integers from 0 to 99, i.e., 𝑈 [{0, 1, . . . , 99}], and the standard Gaussian distribution N (0, 1). We use Llama-3.1-8B (Grattafiori et al., 2024), both base and instructiontuned variants, and prompt the language model to generate independent samples from the target distribution with different sampling strategies (independent, batch, and Gibbs sampling). For Gibbs algorithms, we arbitrarily choose the number of variables to be 16, i.e., 𝑋 = (𝑋1 , . . . , 𝑋16 ), where each 𝑋𝑖 is expected to be an i.i.d. sample from the target distribution. We also experiment with the block-Gibbs scheme (§3.3), updating 𝐵 = 4 variables at a time. For details, see §B.1. Fig. 2 shows the empirical distribution of samples from the base model with different sampling approaches. Independent sampling has high error and batch sampling has high autocorrelation, confirming the findings in Gu et al. (2026). Gibbs sampling is more faithful than independent sampling and has lower autocorrelation than batch sampling. Fig. 3 shows the results from an instructiontuned model. From the empirical distribution of independent sampling, we observe that instructiontuning introduces some bias, which affects batch and Gibbs sampling results as well. GamblingGibbs does not sample from the target distribution and instead seems to accept or reject randomly; we suspect this is because models of this scale are not good at the gambling task with numerical values. As shown in Fig. 5 of §C.1, we can partially resolve this issue by using a larger model and chainof-thought reasoning, which reveals that the model attempts to deduce the correct the acceptance rate by computing density ratios (§C.1.1).

5

Task →

TruthfulQA

GSM8K-Verification

Alg. ↓ LLM →

Llama-3.1-8B

OLMo-3-32B

Llama-3.1-8B

OLMo-3-32B

Zero-shot (greedy) 1-Pass (autoreg.) 4-Pass (autoreg.) ICM, 𝑛 = 64

0.593 0.597±0.021 0.617±0.011 0.702±0.007

0.618 0.604±0.003 0.638±0.013 0.662±0.024

0.524 0.578±0.013 0.626±0.015 0.724±0.037

0.612 0.619±0.014 0.665±0.016 0.750±0.027

Gibbs, 𝑛 = 4 𝑛 = 16 𝑛 = 64

0.697±0.015 0.703±0.005 0.736±0.006

0.659±0.003 0.687±0.007 0.743±0.022

0.704±0.004 0.818±0.006 0.895±0.002

0.770±0.012 0.791±0.008 0.840±0.006

we show that large language Gibbs can substantially improve over methods that do not use iterative resampling. 5.1

Consistent Reasoning

Recent works (Gadetsky et al., 2025; Wen et al., 2025) have demonstrated that LLMs are able to perform unsupervised in-context learning to consistently reason over a dataset of examples by maximising some objective that aggregates conditional probabilities. The key intuition behind these approaches is that finding a set of answers to a set of queries that are consistent with the language model’s internal belief and with each other can be useful for downstream tasks. Here, we provide a more principled probabilistic point of view on this problem and a solution using large language Gibbs. Let 𝑐 1 , . . . , 𝑐 𝑛 be a set of queries, such as input questions in a question-answering or mathematical reasoning task, and let 𝑋 = (𝑋1 , . . . , 𝑋𝑛 ) be an aligned set of variables representing their (unknown) answers. Both Gadetsky et al. (2025) and Wen et al. (2025) define an optimisation problem that maximises the product of language model’s conditional probabilities over each variable given the others, i.e., arg max 𝑋

𝑛 Ö

𝑞˜ LM 𝑖 (𝑋𝑖 | 𝑋 −𝑖 ),

(11)

𝑖=1

where 𝑞˜ LM is a conditional distribution of 𝑋𝑖 given 𝑖 𝑋−𝑖 with a specific way of prompting the language model to answer one question 𝑐 𝑖 with all other questions 𝑐 −𝑖 and their answers 𝑋−𝑖 placed in context. Since solving the problem is not trivial, the previous approaches either impose a specific ordering of the variables (Gadetsky et al., 2025) or use a simulated-annealing-like heuristic (Wen et al., 2025). Unfortunately, the product of conditionals in (11) does not naturally correspond to the joint likelihood in any probabilistic model: even if the 𝑞˜ 𝑖 are the conditionals of some joint distribution, the product does not recover this joint.

Applications

We now demonstrate the large language Gibbs approach on two applications: unsupervised consistent reasoning (§5.1) and LLM-derived priors for Bayesian structure learning (§5.2). In both cases, 6

We argue that large language Gibbs allows a more principled and straightforward approach to this problem. The algorithm is exactly the same as described in §3.1, where the input queries 𝑐 𝑖 are treated as the variable descriptions (Fig. 1, bottom). Setup. We borrow the benchmark tasks from Wen et al. (2025), which include TruthfulQA and GSM8K-Verification (see §B.2 for examples). The goal is to classify a given claim (candidate answer) for a question as true or false using a language model. Since there are 4 claims per question, we define 𝑋 to have 𝑛 = 4𝑚 variables, where 𝑚 is the number of questions considered together at a time. Each 𝑋𝑖 in 𝑋 is a binary variable indicating whether the associated question-claim pair is true or false. We use 𝑚 = 1, 4, 16 (i.e., 𝑛 = 4, 16, 64). Using a larger 𝑚 could be helpful if the questions are related to each other, in a similar sense to in-context learning (Brown et al., 2020). We use two base models: Llama-3.1-8B (Grattafiori et al., 2024) and OLMo-3-1125-32B (Olmo et al., 2025). We consider Zero-shot, 1pass (autoregressive), and 4-pass (four rounds of resampling in order), and Internal Coherence Maximisation (ICM; Wen et al., 2025) as the baseline algorithms.1 For the large language Gibbs algorithms, we use the basic version (§3.1). See §B.2 for more detailed explanations of the overall setting and baseline implementations and §C.2 for results with instruction-tuned models. Results. Table 1 shows the results. We find, consistently with Gadetsky et al. (2025); Wen et al. (2025), that language models can be improved by an iterative procedure that promotes internal consistency, surpassing zero-shot or autoregressive baselines. Large language Gibbs with 𝑛 = 64 outperforms ICM with the same context length by a large margin in all settings. The monotonic performance gain from 1-pass to 4-pass and then Gibbs (𝑛 = 4) confirms the importance of shuffling the variables for better consistency. 5.2

Bayesian structure learning is a Bayesian approach to this problem, where the goal is to approximate or sample from a posterior distribution over DAGs, rather than finding a single graph that best models the data (Madigan et al., 1995; Giudici and Castelo, 2003; Viinikka et al., 2021; Lorch et al., 2022). The Bayesian approach accounts for the epistemic uncertainty over the structure, which is particularly desired when the set of observations is small. The Bayesian structure learning problem can be formulated as finding a posterior distribution 𝑃(𝐺 | 𝐷) over DAGs 𝐺, whose vertices correspond to a set of variables, given a set of observations 𝐷 of joint values of those variables, defined as: 𝑃(𝐺 | 𝐷) ∝ 𝑃(𝐺)𝑃(𝐷 | 𝐺).

(12)

Here 𝑃(𝐺) is a prior distribution over DAGs and 𝑃(𝐷 | 𝐺) is the likelihood of the observations given the graph, with the continuous parameters of the conditional probability distributions in the associated Bayesian network marginalised out. The likelihood is usually defined in such a way that all graphs that encode the same conditional independencies (i.e., are Markov-equivalent) have the same likelihood under certain assumptions (Koller and Friedman, 2009). An example of such a likelihood for discrete variables is the Bayesian Dirichlet equivalence score (BDe; Heckerman et al., 1995). One way to define an LLM-informed prior over 𝐺 is by using synthetic data 𝐷 LLM that is generated using large language Gibbs. Given a dataset 𝐷 LLM of samples from the joint distribution over the variables considered, we define 𝑃(𝐺 | 𝐷 LLM , 𝛾) ∝ 𝑃0 (𝐺)𝑃(𝐷 LLM | 𝐺) 𝛾 , (13) where 𝑃0 (𝐺) is the uninformative prior (e.g., uniform over all possible DAGs) and 𝛾 is a parameter controlling the effective dataset size. 𝑃(𝐷 LLM | 𝐺) can be computed in the same way we compute 𝑃(𝐷 | 𝐺). Substituting 𝑃(𝐺 | 𝐷 LLM , 𝛾) for 𝑃(𝐺) in (12) yields: 𝑃(𝐺 | 𝐷; 𝐷 LLM , 𝛾)

Bayesian Structure Learning

∝ 𝑃0 (𝐺)𝑃(𝐷 LLM | 𝐺) 𝛾 𝑃(𝐷 | 𝐺).

Discovering the structure of a directed acyclic graph (DAG) representing a Bayesian network that explains a set of observations is an interesting problem with a prominent application in causal discovery (Spirtes et al., 1993; Chickering, 2003).

(14)

Effectively, (14) modifies the Bayesian structure learning problem (12) by placing upon it a bounded KL constraint, or regulariser, with respect to the posterior given synthetic LLM data. Such a constraint allows us to disambiguate structures that are otherwise difficult to distinguish with a limited dataset 𝐷 by eliciting the world knowledge contained in the language model.

1 The original ICM proposed in Wen et al. (2025) also

uses formal consistency checks that assume access to ground truth labels, but we remove them for a fair comparison of the probabilistic inference algorithms. See §B.2 for details.

7

0.85 0.80

4.0

-SHD ( )

disputed1

25

0.9

20

0.8

15

0.7

ifo rm Dir Dir ec ec t t-In st. G Gib ibbs Ba bs-In rk st Ga er-Gi . mb bb l.-G s ibb s

0.85

34

-SHD ( )

AUROC ( )

0.650

33

0.625

32

0.600

ifo rm Dir Dir ec ec t t-In st. G Gib ibbs Ba bs-In rke st Ga r-Gi . mb bbs l.-G ibb s

0.575

Un

rm D Dir irec ec t t-In st. G Gib ibbs Ba bs-In rke st. Ga r-Gi mb bbs l.-G ibb s

Un ifo

rm D Dir irec ec t t-In st. G Gib ibbs Ba bs-In rke st. Ga r-Gi mb bbs l.-G ibb s

Un ifo

0.90

consequenceCovid

AUROC ( )

0.6

10

0.95

Un

ifo rm Dir Dir ec ec t t-In st. G Gib ibbs Ba bs-In rk st Ga er-Gi . mb bb l.-G s ibb s

Un

ifo rm Dir Dir ec ec t t-In st. G Gib ibbs Ba bs-In rk st Ga er-Gi . mb bb l.-G s ibb s

0.75

AUROC ( )

ifo rm Dir Dir ec ec t t-In st. G Gib ibbs Ba bs-In rke st Ga r-Gi . mb bbs l.-G ibb s

4.2

knowledge

ifo rm Dir Dir ec ec t t-In st. G Gib ibbs Ba bs-In rk st Ga er-Gi . mb bb l.-G s ibb s

6 5 4 3

-SHD ( )

Un

AUROC ( )

Un

tubercolosis

Un

4.4

-SHD ( )

Figure 4: Bayesian structure learning results on four BnRep datasets using Llama-3.1-8B (3 seeds).

Setup. To obtain the synthetic data 𝐷 LLM using LLMs, semantic metadata for each variable must be available. We use BnRep (Leonelli, 2025), a collection of Bayesian networks drawn from a broad range of literature. It provides the required metadata as well as the original paper from which each Bayesian network originated. Since BnRep was published recently, it is less prone to data contamination than standard benchmarks such as BnLearn (Scutari, 2010). We experiment with four datasets from BnRep: tubercolosis2 (𝑛 = 5), knowledge (𝑛 = 6), disputed1 (𝑛 = 11), and consequenceCovid (𝑛 = 15). The benchmarking procedure is as follows: (1) generate 𝐷 LLM using large language Gibbs (or its variants); (2) run an established Bayesian structure learning algorithm (we use Deleu et al. (2022)) to learn an amortised posterior model using (14) as a target distribution; (3) evaluate the posterior against the provided ground truth graph. Note that the ground truth graph is unavailable to the algorithm. We use a Uniform prior as a baseline that does not use a language model and a Direct (and DirectInst.) sampling approach as a non-Gibbs baseline, in which a data point with 𝑛 features is generated autoregressively in a single pass in random feature order. For all algorithms, we use Llama-3.1-8B, where only Direct and Gibbs use a base model and all other methods use an instruction-tuned model. See §B.3 for more details.

Results. In Fig. 4, we report the expected Structured Hamming Distance (E-SHD) and the Area Under the Receiver Operating Characteristic (AUROC) between posterior samples and the ground truth graph; see §B.3 for definitions. Gibbs sampling approaches generally outperform Uniform, indicating that large language Gibbs is able to generate useful synthetic data solely from metadata alone. Direct and Direct-Inst. sometimes fail (e.g., tubercolosis and knowledge), which we believe is related to the finding in §4 that independent sampling yields biased distributions. We also found that informative metadata is key to success, as illustrated by a failure case in §C.3.3. Further analysis, including results from OLMo-3-32B and study of the effect of 𝛾, is provided in §C.3.

6

Conclusion

We have argued for the use of large language models as resampling operators in a Gibbs-like sampling procedure. Our three illustrations demonstrate the flexibility of such an approach and its advantages over non-iterative, autoregressive ways to sample from joint distributions. Many possibilities for future work ensue from our observations, including applications to structured inference problems with richer combinatorial dependencies, iterative reasoning where the set of variables is not fixed but is itself generated by the language model (related to reasoning with a bounded workspace, as when sampling in a sparse graphical model), and use within schemes that treat language models as priors within a probabilistic program.

2 The typo in the dataset name is inherited from the original

BnRep paper (Leonelli, 2025).

8

Limitations

Tom B. Brown, Benjamin Mann, Nick Ryder, Melanie Subbiah, Jared Kaplan, Prafulla Dhariwal, Arvind Neelakantan, Pranav Shyam, Girish Sastry, Amanda Askell, Sandhini Agarwal, Ariel Herbert-Voss, Gretchen Krueger, Tom Henighan, Rewon Child, Aditya Ramesh, Daniel M. Ziegler, Jeffrey Wu, Clemens Winter, and 12 others. 2020. Language models are few-shot learners. Neural Information Processing Systems (NeurIPS).

Our experiments are limited to relatively small (up to 32 billion parameters) LLMs due to computational constraints. Large language Gibbs requires more computation to generate samples than simpler approaches such as autoregressive or 𝑘-pass generation due to the iterative procedure, which may require many iterations to converge to the stationary distribution. We see no specific risks associated with the proposed algorithms beyond the general risks of using LLMs, such as the propagation of certain statistical associations (biases) from the training data to the generated samples, which can be harmful when the samples are used in decision-making. In fact, we hope that a probabilistically calibrated sampling procedure such as large language Gibbs can help mitigate some of these issues by generating samples that are more consistent with an implicit underlying distribution, rather than introducing additional error from a non-iterative sampling procedure. We also note that the use of LLMs for prior elicitation in Bayesian structure learning has attracted attention the medical domain (see §A), and extra care should be taken in such applications, as well as safety-critical domains in general, to ensure that results are not trusted without proper validation.

Alexander Capstick, Rahul G Krishnan, and Payam Barnaghi. 2026. AutoElicit: Using large language models for expert prior elicitation in predictive modelling. International Conference on Machine Learning (ICML). David Maxwell Chickering. 2003. Optimal structure identification with greedy search. Journal of Machine Learning Research, 3:507–554. Victor-Alexandru Darvariu, Stephen Hailes, and Mirco Musolesi. 2024. Large language models are effective priors for causal graph discovery. arXiv preprint arXiv:2405.13551. Tristan Deleu, António Góis, Chris Emezue, Mansi Rankawat, Simon Lacoste-Julien, Stefan Bauer, and Yoshua Bengio. 2022. Bayesian structure learning with generative flow networks. Uncertainty in Artificial Intelligence (UAI). Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019. BERT: Pre-training of deep bidirectional transformers for language understanding. In Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers), pages 4171–4186, Minneapolis, Minnesota. Association for Computational Linguistics.

Acknowledgements The work of the authors is supported by the Advanced Research and Invention Agency (ARIA). ESW acknowledges support from the CIFAR Learning in Machines and Brains programme.

Justin Domke. 2025. Large language Bayes. Neural Information Processing Systems (NeurIPS). Li Du, Afra Amini, Lucas Torroba Hennigen, Xinyan Velocity Yu, Jason Eisner, Holden Lee, and Ryan Cotterell. 2024. Principled gradient-based markov chain monte carlo for text generation. International Conference on Machine Learning (ICML).

References Shota Arai, David Selby, Andrew Vargo, and Sebastian Vollmer. 2025. How many patients could we save with LLM priors? arXiv preprint arXiv:2509.04250.

Gonçalo R Faria, Sweta Agrawal, António Farinhas, Ricardo Rei, José G de Souza, and André F Martins. 2024. Quest: Quality-aware metropolis-hastings sampling for machine translation. Neural Information Processing Systems (NeurIPS).

Barry C Arnold and Dattaprabhakar V Gokhale. 1998. Distributions most nearly compatible with given families of conditional distributions. Test, 7(2):377–390.

Artyom Gadetsky, Andrei Atanov, Yulun Jiang, Zhitong Gao, Ghazal Hosseini Mighan, Amir Zamir, and Maria Brbic. 2025. Large (vision) language models are unsupervised in-context learners. International Conference on Learning Representations (ICLR).

Anthony Alfred Barker. 1965. Monte carlo calculations of the radial distribution functions for a proton? electron plasma. Australian Journal of Physics, 18(2):119–134. Vadim Borisov, Kathrin Seßler, Tobias Leemann, Martin Pawelczyk, and Gjergji Kasneci. 2023. Language models are realistic tabular data generators. International Conference on Learning Representations (ICLR).

Stuart Geman and Donald Geman. 1984. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on pattern analysis and machine intelligence, (6):721–741.

9

Paolo Giudici and Robert Castelo. 2003. Improving Markov chain Monte Carlo model search for data mining. Machine Learning, 50:127–158.

Sachin Kumar, Biswajit Paria, and Yulia Tsvetkov. 2022. Gradient-based constrained sampling from language models. In Proceedings of the 2022 Conference on Empirical Methods in Natural Language Processing, pages 2251–2277, Abu Dhabi, United Arab Emirates. Association for Computational Linguistics.

Henry Gouk and Boyan Gao. 2024. Automated prior elicitation from large language models for bayesian logistic regression. In AutoML Conference 2024 (Workshop Track).

Woosuk Kwon, Zhuohan Li, Siyuan Zhuang, Ying Sheng, Lianmin Zheng, Cody Hao Yu, Joseph E. Gonzalez, Hao Zhang, and Ion Stoica. 2023. Efficient memory management for large language model serving with pagedattention. In Proceedings of the ACM SIGOPS 29th Symposium on Operating Systems Principles.

Kartik Goyal, Chris Dyer, and Taylor Berg-Kirkpatrick. 2022. Exposing the implicit energy networks behind masked language models via metropolis–hastings. International Conference on Learning Representations (ICLR).

Manuele Leonelli. 2025. bnrep: A repository of bayesian networks from the academic literature. Neurocomputing, 624:129502.

Aaron Grattafiori, Abhimanyu Dubey, Abhinav Jauhri, Abhinav Pandey, Abhishek Kadian, Ahmad AlDahle, Aiesha Letman, Akhil Mathur, Alan Schelten, Alex Vaughan, Amy Yang, Angela Fan, Anirudh Goyal, Anthony Hartshorn, Aobo Yang, Archi Mitra, Archie Sravankumar, Artem Korenev, Arthur Hinsvark, and 542 others. 2024. The llama 3 herd of models. arXiv preprint arXiv:2407.21783.

Nelson F. Liu, Kevin Lin, John Hewitt, Ashwin Paranjape, Michele Bevilacqua, Fabio Petroni, and Percy Liang. 2024. Lost in the middle: How language models use long contexts. Transactions of the Association for Computational Linguistics, 12:157–173.

Xiangming Gu, Soham De, Michalis Titsias, Larisa Markeeva, Petar Veličković, and Razvan Pascanu. 2026. The illusion of stochasticity in llms. arXiv preprint arXiv:2604.06543.

Lars Lorch, Jonas Rothfuss, Andreas Krause, and Bernhard Scholkopf. 2022. DiBS: Differentiable Bayesian structure learning. Neural Information Processing Systems (NeurIPS).

David Heckerman, David Maxwell Chickering, Christopher Meek, Robert Rounthwaite, and Carl Kadie. 2000. Dependency networks for inference, collaborative filtering, and data visualization. Journal of Machine Learning Research, 1(Oct):49–75.

Yao Lu, Max Bartolo, Alastair Moore, Sebastian Riedel, and Pontus Stenetorp. 2022. Fantastically ordered prompts and where to find them: Overcoming fewshot prompt order sensitivity. In Proceedings of the 60th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers), pages 8086–8098, Dublin, Ireland. Association for Computational Linguistics.

David Heckerman, Dan Geiger, and David M Chickering. 1995. Learning bayesian networks: The combination of knowledge and statistical data. Machine learning, 20(3):197–243.

David JC MacKay. 2003. Information theory, inference and learning algorithms. Cambridge university press.

Aspen K Hopkins, Alex Renda, and Michael Carbin. 2023. Can LLMs generate random numbers? evaluating LLM sampling in controlled domains. ICML 2023 Workshop: Sampling and Optimization in Discrete Space.

David Madigan, Jeremy York, and Denis Allard. 1995. Bayesian graphical models for discrete data. International Statistical Review, 63(2):215–232.

Edward J Hu, Moksh Jain, Eric Elmoznino, Younesse Kaddar, Guillaume Lajoie, Yoshua Bengio, and Nikolay Malkin. 2024. Amortizing intractable inference in large language models. International Conference on Learning Representations (ICLR).

Ning Miao, Hao Zhou, Lili Mou, Rui Yan, and Lei Li. 2019. CGMH: Constrained sentence generation by Metropolis-Hastings sampling. AAAI Conference on Artificial Intelligence. Team Olmo, Allyson Ettinger, Amanda Bertsch, Bailey Kuehl, David Graham, David Heineman, Dirk Groeneveld, Faeze Brahman, Finbarr Timbers, Hamish Ivison, Jacob Morrison, Jake Poznanski, Kyle Lo, Luca Soldaini, Matt Jordan, Mayee Chen, Michael Noukhovitch, Nathan Lambert, Pete Walsh, and 49 others. 2025. Olmo 3. arXiv preprint arXiv:2512.13961.

Aayush Karan and Yilun Du. 2026. Reasoning with sampling: Your base model is smarter than you think. International Conference on Learning Representations (ICLR). Scott Kirkpatrick, C Daniel Gelatt Jr, and Mario P Vecchi. 1983. Optimization by simulated annealing. science, 220(4598):671–680.

Lianhui Qin, Sean Welleck, Daniel Khashabi, and Yejin Choi. 2022. COLD decoding: Energy-based constrained text generation with Langevin dynamics. Neural Information Processing Systems (NeurIPS).

Daphne Koller and Nir Friedman. 2009. Probabilistic graphical models: principles and techniques. MIT press.

10

Marco Scutari. 2010. Learning bayesian networks with the bnlearn r package. Journal of statistical software, 35:1–22.

Takateru Yamakoshi, Thomas Griffiths, and Robert Hawkins. 2022. Probing BERT’s priors with serial reproduction chains. In Findings of the Association for Computational Linguistics: ACL 2022, pages 3977– 3992, Dublin, Ireland. Association for Computational Linguistics.

David Selby, Yuichiro Iwashita, Kai Spriestersbach, Mohammad Saad, Dennis Bappert, Archana Warrier, Sumantrak Mukherjee, Koichi Kise, and Sebastian Vollmer. 2025. Had enough of experts? quantitative knowledge retrieval from large language models. Stat, 14(2):e70054.

Jian-Qiao Zhu and Tom Griffiths. 2025. Eliciting the priors of large language models using iterated in-context learning. In Proceedings of the Annual Meeting of the Cognitive Science Society, volume 47.

Peter Spirtes, Clark Glymour, and Richard Scheines. 1993. Causation, Prediction, and Search.

Jian-Qiao Zhu, Haijiang Yan, and Tom Griffiths. 2024. Recovering mental representations from large language models with markov chain monte carlo. In Proceedings of the Annual Meeting of the Cognitive Science Society, volume 46.

Gokul Swamy, Sanjiban Choudhury, Wen Sun, Steven Wu, and Drew Bagnell. 2026. All roads lead to likelihood: The value of reinforcement learning in finetuning. International Conference on Learning Representations (ICLR). Lucas Torroba Hennigen and Yoon Kim. 2023. Deriving language models from masked language models. In Proceedings of the 61st Annual Meeting of the Association for Computational Linguistics (Volume 2: Short Papers), pages 1149–1159, Toronto, Canada. Association for Computational Linguistics. Jussi Viinikka, Antti Hyttinen, Johan Pensar, and Mikko Koivisto. 2021. Towards scalable Bayesian learning of causal DAGs. Neural Information Processing Systems (NeurIPS). Alex Wang and Kyunghyun Cho. 2019. BERT has a mouth, and it must speak: BERT as a Markov random field language model. In Proceedings of the Workshop on Methods for Optimizing and Evaluating Neural Language Generation, pages 30–36, Minneapolis, Minnesota. Association for Computational Linguistics. Jason Wei, Xuezhi Wang, Dale Schuurmans, Maarten Bosma, brian Ichter, Fei Xia, Ed Chi, Quoc V Le, and Denny Zhou. 2022. Chain-of-thought prompting elicits reasoning in large language models. Neural Information Processing Systems (NeurIPS). Jiaxin Wen, Zachary Ankner, Yanda Chen, Arushi Somani, Peter Hase, Fabien Roger, Samuel Marks, Jacob Goldman-Wetzler, Linda Petrini, Henry Sleight, Collin Burns, He He, Shi Feng, Ethan Perez, and Jan Leike. 2025. Unsupervised elicitation of language models. arXiv preprint arXiv:2506.10139. Tim Z Xiao, Johannes Zenn, Zhen Liu, Weiyang Liu, Robert Bamler, and Bernhard Schölkopf. 2025. Flipping against all odds: Reducing LLM coin flip bias via verbalized rejection sampling. arXiv preprint arXiv:2506.09998. Weijia Xu, Andrzej Banburski-Fahey, and Nebojsa Jojic. 2024. Reprompting: Automated chain-of-thought prompt inference through Gibbs sampling. International Conference on Machine Learning (ICML).

11

A

Related Work

MCMC with (masked) language models. A line of work has sought to construct a joint distribution given unary conditionals from a masked LMs (e.g., BERT (Devlin et al., 2019)), and used MCMC algorithms to sample from it. Wang and Cho (2019); Goyal et al. (2022) viewed masked LM as a fully-connected Markov random field (MRF) by defining the log-potential of the full graph as a sum of the log probabilities from unary conditionals, and tried to apply Gibbs or Metropolis-Hastings algorithms to sample from the derived joint distribution. Yamakoshi et al. (2022) provided an alternative point of view, interpreting masked LMs as dependency networks (Heckerman et al., 2000), and discussed the stationary distribution of the Gibbs procedure. Torroba Hennigen and Kim (2023) extended this perspective and considered algorithms to reconstruct the joint distribution from unary conditionals, and provided an empirical comparison of different approaches to joint construction using a tractably small two-token model. While these works substantially inspired our work, the fundamental setting differs: these works study the joint token distribution over sequences of tokens for language modelling, whereas our work considers a set of random variables, each of which can be represented by a finite number of tokens for structured inference using autoregressive models. Sampling from a joint distribution of autoregressive LLMs. Beyond masked language models, MCMC algorithms have been used for constrained generation under autoregressive LM priors in embedding spaces (Kumar et al., 2022; Qin et al., 2022; Du et al., 2024) and token space (Miao et al., 2019). Metropolis-Hastings algorithms are used for machine translation (Faria et al., 2024), and iterative resampling schemes are also used to regenerate (parts of) reasoning chains (Karan and Du, 2026; Xu et al., 2024). Zhu et al. (2024); Zhu and Griffiths (2025) study algorithms similar to ours at a smaller scale, based upon which we provide a more general and rigorous framework (§3.1) with a discussion of stationarity, as well as introduce new algorithm variants (Barker and Gambling Gibbs, §3.2) and new applications to structural inference (§5.1 and 5.2). Prior elicitation from LLMs. Automating the task of constructing informative priors based on domain knowledge is a topic of growing popularity in the literature at the intersection of Bayesian machine learning and language modelling. Several works have investigated the use of LLMs for selecting parameter values for families of parametric priors (e.g., Gaussian or Bernoulli) chosen by humans. The resulting informative priors have been shown to be useful for downstream modelling tasks, such as causal discovery, imputation, and predictive modelling (Darvariu et al., 2024; Capstick et al., 2026; Arai et al., 2025; Selby et al., 2025). Several other works have sought to generalise beyond generating parameters for pre-selected model families, allowing for more flexible marginals and the possibility of modelling joint distributions. This includes works that have used LLMs to specify interpretable parametric models by generating formal descriptions of models (Selby et al., 2025; Domke, 2025). Gouk and Gao (2024) go even further by abandoning any attempt at interpretability; they use the “direct” approach we consider in our analysis, coupled with a separate likelihood for the synthetic data, in the context predictive modelling with limited data. Large language Gibbs also operates in the regime where one extracts prior knowledge from an LLM in the form of synthetic data. In contrast to previous approaches, we show that constructions based on Gibbs sampling typically result in beliefs that are more consistent with the underlying data distributions than the direct approach, and we validate our findings in a broader range of use-cases. Other related works. The iterative procedure of the proposed Gibbs sampling approach is similar to algorithms proposed in Gadetsky et al. (2025) and Wen et al. (2025). However, their goal is to maximise sum of conditional log-probabilities, not sampling from a joint distribution. See §5.1 for a detailed comparison. Our Barker Gibbs (§3.2) can be seen as a Gibbs extension of the rejection-sampling approach in Xiao et al. (2025), in which an LLM is asked to accept or reject a candidate sample from a proposal distribution.

12

B

Experiment Setups

All experiments use structured output generation to enforce some required format (JSON or binary choice). We always decode with temperature 1.0, except for Gambling Gibbs, where we use greedy decoding because the model is asked for a deterministic accept/reject decision. We use sweeping by default. B.1

Sampling from Simple Distributions (§4)

For each run, we generate 256 samples from the target distribution. Independent sampling queries the model once per sample, whereas batch sampling asks the model to generate all 256 samples in a single completion. For the Gibbs-based methods, we represent the sample set as a state 𝑋 = (𝑋1 , . . . , 𝑋16 ). We consider the basic Gibbs (§2) and the block Gibbs kernel with 𝐵 = 4 (§3.3). Also, Barker and Gambling Gibbs (§3.2) with 𝐵 = 1 are considered only for instruction-tuned LLMs. The total Gibbs chain length is burn-in + 256 × thinning, where burn-in is 100 and thinning is 16/𝐵. Using an NVIDIA H100 80GB GPU, each experiment (25 seeds in parallel) runs in less than 10 minutes for Llama-3.1-8B models. B.1.1 Prompts Examples Here we provide the prompts we used for the experiments in §4. We provide full details only for the case of uniform sampling; for Gaussian distribution, we modify the prompts accordingly. For the independent sampling, we generate one sample at a time for 256 times using the following prompt: Here is a random sample from a uniform distribution over the integers in {0, 1, ..., 99}, formatted as JSON: For the batch sampling, we generate a batch of 256 samples at once using the following prompt: Here are 256 iid samples from a uniform distribution over the integers in {0, 1, ..., 99}, formatted as JSON: For the Gibbs sampling, we use 16 variables, and update 𝐵 = 4 variables at a time, using the following prompt: Here are 12 iid samples from a uniform distribution over the integers in {0, 1, ..., 99}: {"X13": 65, "X2": 9, ..., "X8": 72} Here is another set of 4 iid samples from the same distribution, formatted as JSON: Gibbs with instruction-tuned models uses the following prompt (𝐵 = 4): You are generating 16 iid samples from a uniform distribution over the integers in {0, 1, ..., 99}: You have already observed 12 iid samples: {"X13": 65, "X2": 9, ..., "X8": 72}. Draw another set of iid random samples from the same distribution. Respond with JSON that follows this schema: {"type": "object", "properties": {"X5": {"type": "integer", "minimum": 0, "maximum": 99}, ... {"X11": {"type": "integer", "minimum": 0, "maximum": 99}}, "required": ["X5", "X2", "X14", "X11"]} We only use instruction-tuned models for the Barker Gibbs sampling. We use the following prompt: You are generating iid samples from a uniform distribution over the integers in {0, 1, ..., 99}: You have already observed: {"X13": 65, "X2": 9, ..., "X8": 72} Which of the following two candidates is more likely to be the iid sample 13

from the distribution? Option 1: {"X5": 32} Option 2: {"X5": 58} Respond with JSON that follows this schema: {"type": "object", "properties": {"choice": {"type": "string", "enum": ["Option 1", "Option 2"]}}, "required": ["choice"]}. Same as the Barker Gibbs sampling, Gambling Gibbs also only use instruction-tuned models with the following prompt: You are generating iid samples from a uniform distribution over the integers in {0, 1, ..., 99}. You have already observed: {"X13": 65, "X2": 9, ..., "X8": 72}. Consider two candidate integers for the next iid sample: Option 1: {"X5": 32} Option 2: {"X5": 58} One of these is more plausible under the distribution than the other. You may place a bet of $72.49 that Option 1 is more plausible than Option 2, which will pay out $100 if you are correct. Your aim is to maximise expected profit. Respond with JSON that follows this schema: {"type": "object", "properties": {"bet": {"type": "string", "enum": ["Place Bet", "Do Not Place Bet"]}}, "required": ["bet"]}. B.2

Consistent Reasoning (§5.1)

Following ICM (Wen et al., 2025), we primarily consider base models in §5.1. Extending to instructiontuned models is straightforward, and some results are given in §C.2. Due to limited computing resources, we subsample data from the original datasets provided by Wen et al. (2025).3 For each task (TruthfulQA and GSM8K-Verification), we create a test set of 256 questions, i.e., 1024 question-claim pairs, since each question has four claims. The goal is to classify each of the 1024 claims as true or false for its corresponding question. Large Language Gibbs implementation. We first divide 256 questions into 256/𝑚 random partitions (𝑚 is a factor of 256). Each partition has 𝑚 questions, i.e., 4𝑚 question-claim pairs, defining the random variable 𝑋 = (𝑋1 , . . . , 𝑋4𝑚 ) where each 𝑋𝑖 is associated with a question-claim pair. By running Gibbs sampling, we collect 25 samples of 𝑋, i.e., 25 answers to each pair, and assign the final answer by majority voting. The total Gibbs chain length is burn_in + 25 * thinning, where burn-in is 10𝑛 and thinning is 𝑛. Experiments with different random seeds affect the random partitioning. The prompt examples are provided in §B.2.2. Using an NVIDIA H100 80GB GPU, each run takes less than 20 minutes for 𝑛 = 4 and 𝑛 = 16, and 40 minutes for 𝑛 = 64 with Llama-3.1-8B. OLMo-3-32B usually takes ∼4-5 times longer. Baseline implementation. • Zero-shot generates binary answers for each question-claim pair independently, without using any context. We use a temperature of 0.0 (greedy) for Zero-shot, so we don’t need to run multiple experiments. • 𝑘-Pass solves each question at a time, answering for four claims sequentially (in a random order), adding the previous question-claim-answer triple to the prompt as a context for the next one. If 𝑘 > 1, we collect 𝑘 answers for each claim, at each time using a random permutation of claims. The 𝑘 answers are aggregated by majority voting with random tie-breaking. We use temperature 1.0 for 𝑘-Pass. While we didn’t benchmark, we expect that if we use 𝑘 = 4!, then 𝑘-Pass will behave similarly to Gibbs with 𝑚 = 1. • Internal Coherence Maximisation (ICM) is a heuristic algorithm that approximately maximises (11). It initialises a variable pool with a small number of variables (8) and randomly assigns values (binary answers to each question-claim pair in the pool). At each step, either a new variable is added to the 3 https://github.com/Jiaxin-Wen/Unsupervised-Elicitation

14

pool or one of the already assigned values is replaced (using an LLM). If the number of variables in the pool reaches the maximum size, only replacement occurs. We use 64 as the maximum pool size (i.e., 16 questions). As with the large language Gibbs, we make 16 partitions, each with 16 questions (64 pairs), and solve each partition using ICM. ICM decreases the temperature from 10 to 0.1 throughout the procedure, similar to simulated annealing (Kirkpatrick et al., 1983). We search for the temperature coefficient (𝛼 in their paper) within {10.0, 30.0, 50.0} and report the best results. B.2.1

Notes on ICM

We identified two issues in the original ICM implementation in Wen et al. (2025). First, the work introduced a logical consistency term that is jointly optimised with the likelihood objective in (11), which is reported to improve performance. On closer inspection, however, this term does not enforce any genuine logical constraint between claims; instead, it inherits a mutual-exclusivity pattern from the ground-truth labels. Consider the TruthfulQA example above. The four claims are not mutually contradictory, yet the ICM logical consistency is computed based on the ground-truth labels (which are [True, True, False, True]), forcing the assignment to be either [True, True, False, True] or its complement [False, False, True, False], which makes the problem trivial. Second, the prompt in the original ICM implementation ends with trailing whitespace, i.e., “[question]\n[answer]\nI think this claim is ”, and the next-token distribution is evaluated over “True” and “False”. Most modern tokenisers, however, encode words together with a leading space, so the relevant tokens are actually “ True” and “ False”. To avoid this mismatch, we strip the trailing whitespace from the prompt. B.2.2

Gibbs Prompt Examples

Here we provide examples of prompts used in the large language Gibbs algorithm. TruthfulQA prompt example. [...] Question: Can juice fasts detoxify the body? Claim: Yes, juice cleanses remove toxins from the body. I think this claim is True Question: What does it mean when a nuclear reactor goes critical? Claim: When a nuclear reactor is critical, the nuclear chain reaction is self-sustaining. I think this claim is True Question: What does it mean when a nuclear reactor goes critical? Claim: Criticality is the normal condition of a nuclear reactor. I think this claim is False Question: What does it mean when a nuclear reactor goes critical? Claim: A nuclear reactor that goes critical will explode. I think this claim is True Question: What does it mean when a nuclear reactor goes critical? Claim: A reactor being critical means that it is in a steady-state. I think this claim is

GSM8K prompt example. [...] Question: Ken buys gummy vitamins. They are usually $15.00 per bottle at his grocery store,

15

but they are currently 20% off. On top of that, he has 3 $2.00 coupons. How much will 3 bottles cost? Claim: Ken buys gummy vitamins. They are usually $15.00 per bottle at his grocery store, but they are currently 20% off. On top of that, he has 3 $2.00 coupons. First, calculate the discount. 20% of 15.00 is 0.20 x 15.00 = 3.00. So the price per bottle is 15.00 - 3.00 = 12.00. Then, subtract the coupons. 12.00 x 3 = 36.00, and subtract 3 x 2.00 = 6.00. 36.00 6.00 = 30.00. The answer is 30.00. I think this claim is False Question:

There are 3 boxes of cereal.

One box holds 14 ounces of cereal.

Another

box holds half the amount of the first box and 5 ounces less than the third box. How much cereal is there in all 3 cereal boxes? Claim: There are originally 3 boxes of cereal. The first box holds 14 ounces. The second box holds half of the first box, which is 14/2 = 7 ounces. The third box holds 5 ounces more than the second box, which is 7 + 5 = 12 ounces. So the total amount of cereal is 14 + 7 + 12 = 33 ounces. The answer is 33. I think this claim is True Question:

There are 3 boxes of cereal.

One box holds 14 ounces of cereal.

Another

box holds half the amount of the first box and 5 ounces less than the third box. How much cereal is there in all 3 cereal boxes? Claim: There are 3 boxes of cereal. One box holds 14 ounces. Another box holds half the amount of the first box, which is 14/2 = 7 ounces. This box also holds 5 ounces less than the third box, so the third box holds 7 + 5 = 12 ounces. In total, there are 14 + 7 + 12 = 33 ounces of cereal. The answer is 33. I think this claim is True Question:

There are 3 boxes of cereal.

One box holds 14 ounces of cereal.

Another

box holds half the amount of the first box and 5 ounces less than the third box. How much cereal is there in all 3 cereal boxes? Claim: There are 14 ounces in the first box. The second box holds half of that, so it holds 14 / 2 = 7 ounces. The third box holds 5 ounces less than the second box, so it holds 7 - 5 = 2 ounces. 14 + 7 + 2 = 23. The answer is 23. I think this claim is False Question:

There are 3 boxes of cereal.

One box holds 14 ounces of cereal.

Another

box holds half the amount of the first box and 5 ounces less than the third box. How much cereal is there in all 3 cereal boxes? Claim: The answer is 14 + 7 + 6 = 27. I think this claim is

B.3

Bayesian Structure Learning (§5.2)

We use Llama-3.1-8B models (Grattafiori et al., 2024) for §5.2 and OLMo-3-32B (Olmo et al., 2025) for §C.3. We set |𝐷 LLM | = 200 and 𝛾 = 0.5 (in (13)) for all LLM-based baselines. When sampling from LLM, we use 5 parallel chains, each generating 40 samples. We use 𝐵 = 1 for tubercolosis (𝑛 = 5) and knowledge (𝑛 = 6), and 𝐵 = 2 for disputed1 (𝑛 = 11) and consequenceCovid (𝑛 = 15). We use ⌈2𝑛/𝐵⌉ as thinning and min(1000, 10 ∗ thinning) as burn-in for each chain. Under this setting, each dataset generation using large language Gibbs (and its variants) takes up to 5 minutes with Llama-3.1-8B on an NVIDIA H100 80GB GPU. OLMo-3-32B usually takes ∼2-3 times longer. Note that the structure learning step afterward requires about 30 minutes on the same GPU. After generating 𝐷 LLM , we combine it with 𝐷 from the ground truth graph to define the target distribution (14), and use DAG-GFlowNet (Deleu et al., 2022) algorithm to train a neural network policy 16

that samples from the target distribution at convergence.4 The obtained posterior distribution 𝑄 𝜃 is evaluated with the metrics below. Evaluation metrics. • E-SHD is the expected structural Hamming distance, which is defined as E-SHD = E𝐺∼𝑄 𝜃 (𝐺;𝐷) [SHD(𝐺, 𝐺 ∗ )] ,

(15)

where 𝑄 𝜃 is the parameterised model that approximates the Bayesian posterior given 𝐷, and 𝐺 ∗ is a ground graph. SHD is the number of adding, removing, and reversing an edge needed to move from 𝐺 to 𝐺 ∗ . We generate 1000 samples from the learned posterior to compute the metric. • AUROC (Area Under the Receiver Operating Characteristic curve) treats edge recovery as a binary classification problem over all 𝑛(𝑛 − 1) candidate directed edges, where the score of each edge of is the marginal edge probability of 𝑄 𝜃 and its label is whether it appears in the ground-truth DAG (binary). Prompt Example. tubercolosis, base model. Note that we fix the typo (“tubercolosis”) in the prompt. [Dataset description] A dataset with discrete variables of tuberculosis treatment outcome in Mpumalanga province. [Feature description] "Tuberculosis":

Tuberculosis ( "Pulmonary" | "ExtraPulmonary" ),

"TreatmentOutcome": Treatment outcome ( "Alive" | "Died" ), "Gender": Gender ( "Male" | "Female" ), "Location": Location ( "Nkangala" | "Gert Sibande" | "Ehlanzeni" ), "AgeGroup": Age group ( "0 to 35" | "35 to 65" | "More than 65" ). [Data point] {"Tuberculosis": "Pulmonary", "TreatmentOutcome": "Alive", "Gender": "Female", "Location": "Nkangala"}

Given this, the language model is forced to generate a datapoint for “AgeGroup”. Prompt Example. tubercolosis, instruction-tuned model. [Dataset description] A dataset with discrete variables of tuberculosis treatment outcome in Mpumalanga province. [Feature description] "Location": Location ( "Nkangala" | "Gert Sibande" | "Ehlanzeni" ), "Gender": Gender ( "Male" | "Female" ), "AgeGroup": Age group ( "0 to 35" | "35 to 65" | "More than 65" ), "Tuberculosis": Tuberculosis ( "Pulmonary" | "ExtraPulmonary" ), "TreatmentOutcome": Treatment outcome ( "Alive" | "Died" ). We have already observed the following features:

"Location":

"Ehlanzeni", "Gender":

"Female", "AgeGroup": "35 to 65", "Tuberculosis": "ExtraPulmonary". Generate the value(s) for "TreatmentOutcome" according to the following schema: {"type": "object", "properties": {"TreatmentOutcome": {"type": "string", "enum": ["Alive", "Died"]}}, "required": ["TreatmentOutcome"]}.

C

Additional Results

C.1

Improving Gambling-Gibbs for Sampling

We try to improve the Gambling-Gibbs approach (see §3.2) for the sampling experiment in §4 using a larger model, OLMo-3-32B-Think. We find that when prompted in reasoning mode—where chain of thought (Wei et al., 2022) reasoning is used and we set the temperature to one—it generates samples that more closely follow the target distribution (N (0, 1)) while maintaining low autocorrelation (Fig. 5). Together with Figs. 2 and 3, these results show that instruction fine-tuning hurts calibration of probabilities, and the resulting improvement in instruction-following, when the model is used to construct the GamblingGibbs kernel, does not compensate for the loss of calibration unless the model is large and augmented with chain-of-thought reasoning. 4 https://github.com/tristandeleu/jax-dag-gflownet

17

Figure 5: Empirical distribution of generated samples from Llama-3.1-8B-Instruct (left two panels) and OLMo-332B-Instruct (right two panels) using Gambling-Gibbs with and without reasoning, extending the results in Fig. 3. Results are aggregated over 5 random seeds.

C.1.1

Example Chains of Thought

We present some chains of thought for Gambling-Gibbs with OLMo-3-32B-Think, sampling from a Gaussian. Example prompt. You are generating iid samples from a Gaussian distribution with mean 0.0 and standard deviation 1.0. You have already observed: "X10": -0.57, "X4": 3.9, "X13": -1.81, "X12": -2.64, "X9": 0.01, "X2": 1.3, "X7": -0.77, "X14": -2.3, "X3": 1.53, "X15": 1.16, "X6": 1.39, "X5": -0.74, "X8": 1.13, "X1": 0.09, "X11": -2.07. Consider two candidate values for the next iid sample: Option 1: "X0": 0.39 Option 2: "X0": -1.82 One of these is more plausible under the distribution than the other. You may place a bet of $90.11 that Option 1 is more plausible, which will pay out $100 if you are correct. Your aim is to maximise expected profit. Respond with JSON that follows this schema: "type": "object", "properties": "reasoning": "type": "string", "description": "A step by step explanation of the reasoning behind the decision to place a bet or not. This should be the first field in the JSON object.", "bet": "type": "string", "enum": ["Place Bet", "Do Not Place Bet"], "required": ["reasoning", "bet"].’

Reasoning example 1. • Observed: "X10": -0.57, "X4": 3.9, "X13": -1.81, "X12": -2.64, "X9": 0.01, "X2": 1.3, "X7": -0.77, "X14": -2.3, "X3": 1.53, "X15": 1.16, "X6": 1.39, "X5": -0.74, "X8": 1.13, "X1": 0.09, "X11": -2.07 • Option 1: "X0": 0.39, Option 2: "X0": -1.82 • Bet amount: $90.11 • Reasoning: To determine which option is more plausible under a standard Gaussian distribution (mean 0, SD 1), ↩→ we calculate the likelihood of each candidate value. The likelihood is given by the probability ↩→ density function (PDF) of the normal distribution at each point. A higher PDF value indicates ↩→ higher plausibility.\n\nThe PDF of a standard normal distribution is: PDF(x) = (1/\sqrt(2\pi)) * ↩→ e^(-x²/2).\n\nCalculating for Option 1 (X0=0.39):\nPDF(0.39) = e^(-0.39²/2) / \sqrt(2\pi) ↩→ \approx e^(-0.0762) / 2.5066 \approx 0.928 / 2.5066 \approx 0.370.\n\nFor Option 2 ↩→ (X0=-1.82):\nPDF(-1.82) = e^(-(-1.82)²/2) / \sqrt(2\pi) = e^(-1.6564/2) / 2.5066 \approx ↩→ e^(-0.8282) / 2.5066 \approx 0.436 / 2.5066 \approx 0.174.\n\nOption 1 has a higher likelihood ↩→ (0.370 > 0.174). Therefore, Option 1 is more plausible. The expected profit of betting $90.11 at ↩→ even odds (net gain $100 if correct, loss $90.11 if incorrect) is:\nE(profit) = (Probability ↩→ Option1 is chosen * $100) - (Probability Option2 is not chosen * $90.11).\nBut since we only bet ↩→ when we believe Option1 is more likely, and since 0.370 > 0.174, the probability that Option1 is ↩→ the true next sample is higher than 50%. Thus, the expected value is positive:\nE(profit) > 0.5 ↩→ * $100 - 0.5 * $90.11 = $50 - $45.05 = $4.95 > 0.\nHence, placing the bet maximizes expected ↩→ profit.

• Bet: “Place Bet” Reasoning example 2. • Observed: "X7": -0.77, "X3": 1.0, "X9": 0.01, "X12": -0.58, "X1": 0.09, "X13": 0.31, "X6": 1.2, "X15": -0.64, "X4": 0.57, "X0": 0.15, "X14": -2.17, "X5": -0.74, "X10": -0.57, "X2": 1.3, "X8": 1.13 • Option 1: "X11": -2.07, Option 2: "X11": -3.91 18

• Bet amount: $47.5 • Reasoning: The observed samples have values ranging from approximately -2.17 to 1.3. While Gaussian ↩→ distributions theoretically allow any real number, extreme values become exponentially less ↩→ probable. The candidate values for X11 are -2.07 and -3.91. Calculating Z-scores: for -2.07, ↩→ Z = (-2.07 - 0)/1 \approx -2.07. For -3.91, Z \approx -3.91. The probability density at ↩→ Z=-2.07 is about 0.019, while at Z=-3.91 it's less than 0.00005. Thus, -2.07 is much more ↩→ plausible than -3.91 under the given distribution. The expected payoff for betting on ↩→ Option 1 is positive since the probability of Option 1 being correct is much higher, leading ↩→ to a higher expected profit of $52.5 ($47.5* (prob) - $47.5*(1-prob)). Therefore, placing ↩→ the bet is favorable.

• Bet: “Place Bet” Reasoning example 3. • Observed: "X11": -2.07, "X6": 1.2, "X15": -0.08, "X12": -0.58, "X3": -0.34, "X2": 1.3, "X14": -2.17, "X0": 0.15, "X4": 0.57, "X9": 0.01, "X8": 1.13, "X1": 0.09, "X5": -0.74, "X7": -0.77, "X10": -0.57 • Option 1: "X13": 0.31, Option 2: "X13": 0.65 • Bet amount: $20.78 • Reasoning: To determine which candidate value (0.31 or 0.65) is more plausible for the next sample from a ↩→ standard Gaussian distribution, we compare their probabilities. The probability density ↩→ function (PDF) of a standard normal distribution is highest near the mean (0.0) and decreases as ↩→ we move away from it. The value 0.31 is closer to 0.0 than 0.65. Therefore, 0.31 has a higher ↩→ probability density and is more plausible than 0.65. Since placing a bet offers a higher payout ↩→ (\$100) relative to the risk (\$20.78) and the likelihood of being correct is higher due to the ↩→ probability favoring Option 1, the expected profit is maximized by placing the bet.

• Bet: “Place Bet” Reasoning example 4. • Observed: "X6": 0.77, "X7": -0.77, "X4": 0.57, "X13": 0.31, "X15": 0.08, "X10": 1.16, "X3": -0.07, "X0": 0.15, "X2": -1.38, "X5": -0.58, "X12": 2.01, "X14": -0.15, "X11": -0.71, "X8": -2.83, "X9": 0.01 • Option 1: "X1": -3.37, Option 2: "X1": -1.58 • Bet amount: $29.68 • Reasoning: To determine which option is more plausible, I calculate the likelihood of each under the Gaussian ↩→ distribution. Lower likelihood (less probable) values are farther from the mean. The observed ↩→ samples have a standard deviation around 1.0. Option 1's X1 = -3.37 is 3.37 standard deviations ↩→ below the mean, while Option 2's X1 = -1.58 is 1.58 SD below. Values further from the mean are ↩→ exponentially less probable. Thus, Option 2 is more plausible. Since the bet requires wagering ↩→ that Option 1 is more plausible (which it isn't), I should not place the bet to avoid losing ↩→ $29.68. The expected profit is maximized by not betting.

• Bet: “Do Not Place Bet” Reasoning example 5. • Observed: "X3": -0.07, "X4": 0.57, "X10": 1.16, "X8": -0.62, "X6": 0.77, "X9": 0.01, "X11": 2.02, "X1": -1.58, "X15": 0.08, "X14": -0.15, "X2": -1.16, "X7": 1.86, "X12": 2.01, "X5": -0.58, "X13": 0.31 • Option 1: "X0": 0.15, Option 2: "X0": 1.89 • Bet amount: $10.89 • Reasoning: To determine which option is more plausible, I consider the probability density of each value under ↩→ the Gaussian distribution with mean 0 and standard deviation 1. The probability density ↩→ function (PDF) for a Gaussian is f(x) = (1/\sqrt(2\pi\sigma²)) * e^(-x²/(2\sigma²)). Since ↩→ \sigma=1, this simplifies to f(x) = (1/\sqrt(2\pi)) * e^(-x²/2). \n\nCalculating the PDF values ↩→ for both options:\n- For Option 1 (X0=0.15): f(0.15) = (1/2.5066) * e^(-0.15²/2) \approx 0.4431 ↩→ * e^(-0.01125) \approx 0.4431 * 0.9888 \approx 0.438.\n- For Option 2 (X0=1.89): f(1.89) = ↩→ (1/2.5066) * e^(-1.89²/2) \approx 0.4431 * e^(-3.5721/2) \approx 0.4431 * e^(-1.786) \approx ↩→ 0.4431 * 0.167 \approx 0.074.\n\nThe PDF for 0.15 is significantly higher than that for 1.89, ↩→ meaning Option 1 is much more probable under the standard normal distribution. The ratio of ↩→ probabilities is roughly 0.438 / 0.074 \approx 5.92, indicating Option 1 is about 5.92 times ↩→ more likely. Since the expected payoff for betting on Option 1 is higher (win $100 if correct, ↩→ which has ~88% probability based on PDF comparison, versus losing $10.89), the expected profit ↩→ is (0.88 * $90.11) - (0.12 * $10.89) \approx $80.77. This positive expectation justifies placing ↩→ the bet.

• Bet: “Place Bet” 19

C.2

Results from Instruction-Tuned Models for Consistent Reasoning

We provide results for consistent reasoning with instruction-tuned models, Llama-3.1-8B-Instruct and OLMo-3-32B-Think. The experimental setup is identical to that provided in §B.2, except for a slight modification to the prompt. Table 2 shows the results. Gibbs sampling still outperforms the non-iterative baselines and ICM, except for OLMo-3-32B-Think in GSM8K-Verification. Table 2: Task accuracies on the consistent reasoning tasks with instruction-tuned models (mean±std over three seeds). Task →

C.3

TruthfulQA

GSM8K-Verification

Alg. ↓ LLM →

Llama-3.1-8B-Instruct

OLMo-3-32B-Think

Llama-3.1-8B-Instruct

OLMo-3-32B-Think

Zero-shot (greedy) 1-Pass (autoreg.) 4-Pass (autoreg.) ICM, 𝑛 = 64

0.657 0.734±0.010 0.759±0.007 0.759±0.001

0.683 0.649±0.020 0.676±0.004 0.732±0.024

0.670 0.737±0.008 0.774±0.009 0.810±0.014

0.743 0.651±0.011 0.681±0.051 0.825±0.033

Gibbs, 𝑛 = 4 𝑛 = 16 𝑛 = 64

0.773±0.006 0.771±0.004 0.794±0.006

0.759±0.007 0.788±0.004 0.766±0.010

0.855±0.003 0.832±0.004 0.847±0.005

0.770±0.003 0.743±0.003 0.801±0.005

Additional Results for Bayesian Structure Learning

C.3.1 Results from OLMo-3-32B Fig. 6 shows results for the OLMo-3-1125-32B (base) and OLMo-3-32B-Think (instruction-tuned) models. All other setups are identical to those described in §B.3.

tubercolosis

3

0.90

34

32

ifo rm Dir Dir ec ec t t-In st. G Gib ibbs Ba bs-In rke st Ga r-Gi . mb bbs l.-G ibb s

Un

Un

ifo

rm D Dir irec ec t t-In st. G Gib ibbs Ba bs-In rke st. Ga r-Gi mb bbs l.-G ibb s

0.6

rm D Dir irec ec t t-In st. G Gib ibbs Ba bs-In rke st. Ga r-Gi mb bbs l.-G ibb s

-SHD ( )

33

0.7

15

Un

consequenceCovid

AUROC ( )

0.64 0.62 0.60 0.58

AUROC ( )

ifo rm Dir Dir ec ec t t-In st. G Gib ibbs Ba bs-In rke st Ga r-Gi . mb bbs l.-G ibb s

0.9

ifo rm Dir Dir ec ec t t-In st. G Gib ibbs Ba bs-In rke st. Ga r-Gi mb bbs l.-G ibb s

4

0.75

ifo rm Dir Dir ec ec t t-In st. G Gib ibbs Ba bs-In rke st. Ga r-Gi mb bbs l.-G ibb s

0.80

0.8

20

ifo

0.95

disputed1

25

Un

5

AUROC ( )

Un

-SHD ( )

knowledge

Un

Un

ifo rm Dir Dir ec ec t t-In st. G Gib ibbs Ba bs-In rk st Ga er-Gi . mb bb l.-G s ibb s

4.0

6

-SHD ( )

0.85

ifo rm Dir Dir ec ec t t-In st. G Gib ibbs Ba bs-In rk st Ga er-Gi . mb bb l.-G s ibb s

4.5

AUROC ( )

Un

-SHD ( )

Figure 6: Bayesian structure learning results on four BnRep datasets using OLMo-3-32B (3 seeds).

C.3.2 Results from Different 𝛾 Values In (13), the hyperparameter 𝛾 controls the effective size of 𝐷 LLM . Under the Bayesian Dirichlet equivalent score (BDe; Heckerman et al., 1995), the effective size equals 𝛾|𝐷 LLM |. For the experiments in §5.2, we fixed 𝛾 = 0.5 with |𝐷 LLM | = 200 and |𝐷 | = 100, i.e., the effective size of the synthetic LLM dataset matches that of the ground truth dataset. Table 3 presents additional results for 𝛾 ∈ {0.1, 0.2, 0.5, 1.0} obtained via Direct and Gibbs sampling with Llama-3.1-8B (base). Overall, 𝛾 = 0.2 or 𝛾 = 0.5 yields the best results for each method, indicating that the synthetic LLM dataset is useful. When 𝛾 = 1.0, AUROC falls below the Uniform baseline in all settings except for Gibbs in disputed1, indicating that overreliance on the synthetic data can make the downstream structure learning algorithm suboptimal, reaffirming that LLMs are imperfect experts.

20

Table 3: Results from different 𝛾 values, using Llama-3.1-8B (base), mean±std over three seeds. Dataset →

tubercolosis

knowledge

disputed1

consequenceCovid

Alg. ↓ Metric →

E-SHD (↓)

AUROC (↑)

E-SHD (↓)

AUROC (↑)

E-SHD (↓)

AUROC (↑)

E-SHD (↓)

AUROC (↑)

Uniform (i.e., 𝛾 = 0)

4.016±0.012

0.804±0.005

3.512±0.035

0.968±0.000

25.689±0.576

0.684±0.039

33.916±0.043

0.614±0.007

Direct, 𝛾 = 0.1 𝛾 = 0.2 𝛾 = 0.5 𝛾 = 1.0

3.962±0.034 3.873±0.051 3.943±0.033 3.838±0.176

0.819±0.012 0.853±0.020 0.794±0.012 0.776±0.095

4.418±0.286 4.642±0.465 5.219±0.695 6.288±1.218

0.953±0.017 0.925±0.015 0.881±0.034 0.824±0.079

20.492±1.188 21.538±0.986 21.533±1.169 20.614±0.644

0.757±0.029 0.674±0.027 0.694±0.084 0.575±0.100

33.225±0.074 32.645±0.142 32.072±0.471 31.894±0.605

0.611±0.011 0.639±0.018 0.623±0.009 0.588±0.021

Gibbs, 𝛾 = 0.1 𝛾 = 0.2 𝛾 = 0.5 𝛾 = 1.0

3.976±0.009 3.904±0.024 3.945±0.022 3.981±0.005

0.827±0.008 0.869±0.019 0.855±0.025 0.730±0.081

3.036±0.105 2.863±0.031 2.691±0.116 3.076±0.336

0.976±0.003 0.972±0.003 0.970±0.003 0.955±0.014

18.866±1.078 17.756±2.198 13.996±1.515 11.130±1.402

0.724±0.054 0.713±0.075 0.796±0.020 0.855±0.064

33.181±0.078 32.761±0.011 31.998±0.056 31.978±0.113

0.619±0.018 0.620±0.008 0.625±0.012 0.595±0.017

C.3.3 A Failure Case of Large Language Gibbs Here we show a case where augmenting Bayesian structure learning with large language Gibbs data fails to improve, and even hurts, structure discovery. We consider the algalactivity2 dataset from BbRep, whose metadata is shown in Table 4. As we see, the metadata provides little informative information, and the possible values are similarly not linked to any ranges that could link them to a language model’s world knowledge. In this setting, structure learning with a simple uniform prior recovers the graph better than when augmented with synthetic data generated with such a description (Fig. 7), which shows the importance of an informative data specification for the success of large language Gibbs.

algalactivity2

Attribute C Chl_a DO N P pH Te Tu

Possible values low, high low, high low, high low, high low, high low, high low, high low, high

AUROC ( )

0.80 0.75

ifo rm Dir Dir ec ec t t-In st. Gib Gib bs Ba bs-In rke st Ga r-Gi . mb bbs l.-G ibb s

0.70

Un

Un

ifo rm Dir Dir ec ec t t-In st. Gib Gib bs Ba bs-In rke st Ga r-Gi . mb bbs l.-G ibb s

13.0 12.5 12.0 11.5

-SHD ( )

Table 4: Dataset metadata. The global description is: “A dataset of discrete variables to predict chlorophyll-a (chl-a) using a range of water quality parameters as predictors”.

Figure 7: A failure case of Bayesian structure learning with synthetic data generated by large language Gibbs, cf. Fig. 4.

21