Large Language Bayes
Authors: Justin Domke
Venue: NeurIPS 2025 poster ()
Keywords: bayesian inference
variational inference
Primary area: probabilistic_methods

Abstract: Many domain experts do not have the time or expertise to write formal Bayesian models. This paper takes an informal problem description as input, and combines a large language model and a probabilistic programming language to define a joint distribution over formal models, latent variables, and data. A posterior over latent variables follows by conditioning on observed data and integrating over formal models. This presents a challenging inference problem. We suggest an inference recipe that amounts to generating many formal models from the large language model, performing approximate inference on each, and then doing a weighted aver- age. This is justified and analyzed as a combination of self-normalized importance sampling, MCMC, and importance-weighted variational inference. Experimentally, this produces sensible predictions from only data and an informal problem description, without the need to specify a formal model.

Large Language Bayes
Justin Domke University of Massachusetts Amherst
Abstract
Many domain experts do not have the time expertise to write formal Bayesian models. This paper takes an informal problem description as input, and combines large language model and probabilistic programming language to define joint distribution over formal models, latent variables, and data. A posterior over latent variables follows by conditioning models. This presents a challenging inference problem. We suggest an inference recipe that amounts to generating many formal models from the large language model, performing approximate inference on each, and then doing a weighted aver- age. This is justified and analyzed combination of self-normalized importance sampling, MCMC, and importance-weite tally, this produces sensible predictions from only data and informal problem description, without the need to specify a formal model.
observed data and integrating over formal
variational inference. Experimen-
Introduction
Why isn't Bayesian inference more popular? Arguably, a major reason is simply that creating probabilistic models hard. Most people interested in analyzing data are neither programmers nor experts in statistics. Yet, writing a probabilistic model requires learning a probabilistic programming language (PPI.,, fluency with a range of statistical distributions, and experience formalizing problem intuitions. Even for experts, this is difficult and error-prone.
natural idea is ask the Language Model (LLM) to generate formal probabilistic model. While LLMs can write serviceable models, practice they struggle in the same way as sesrans--ommimmee the models they create are good and sometimes they aren't tt [27].
describe their problem in plain language and then use a Large
Yet L.MS have major advantage different candidate models for the same problem.
humans: they don't mind being asked create many
The central idea of this paper is to mathematically "glue" an LLM to PPL. Given informal description problem, a joint distribution is defined over (1) formal probabilistic models, (2) observed data, and (3) unobserved target variables. We then condition on data and marginalize out the space formal models to get a final posterior over target variables.
The main contributions of this paper are:
A new problem definition, in which an informal description and dataset define a posterior via an LLM and a PPL. (Sec. 2)
A broad algorithmic recipe for solving the resulting inference problem. (Sec. 3)
Experiments illustrating that the final approximated posterior captures user intent and typically better than taking naive average of formal models. (Sec. 4)
Theory analyzing the expected accuracy of the suggested inference recipe. (Sec. 5)
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Informal text
Data
Marginal likelihoods
I' ve recorded the prscipttation in Nascachusetts each
"precipitacton" 16.
01,
3:36.
benesber 2824. tnerels some aind of patters?
6:01.
9.45.
Weighted
Hake sure to reflect that
Predict
2025.
0.01. 0.07, 6.821
0.50 5 0.04. a.01,
2-02.
Average
precipltation to never regative.
Posteriors
Final posterior
DATA
tst run dayd
days] praciptiatiio
results Becember
Inference
GOA..
tst // oxteree for nect
LLM
Formal models
motel
gannats, 9-11; for 11 Traum days
#odet
#odel {
srecis expenentiaa((u));
prec cauctyle, 1881;
nomalle. 18);
stodey
for 1 in 1:mun L:nua_ preciptiation
sigra - ceuchylE, 2.51:
pracall
prec inwcipitaiion]1 exporettial(eppte))
exproential rngleuty
pric_stodee)
precipitattonit] exporerttate
expinu + #ht-Iprecipptaationn--11 2011))
mdal for It diyun
real mixt pormal_nnjjrec, proc. proc_st00le)
generated quantities {
exponential precisttationi
nuxt feacte, next]:
real expetentiial.rrgt phillpreciptaaatinnnnus
3
peserated quantities
rual exdonent La Lambda);
Figure The basic idea. Given informal user text, an LLM generates a set of candidate formal models. Inference is performed on each and the posteriors are combined with weight proportional to the marginal likelihood. Here four (real) LLM-generated formal models the Stan language are shown in different colors, with corresponding colors for marginal likelihoods and posteriors.
2 The basic idea
The user provides plain- languag description of the problem, including anything they know about the phenomena under study, the types and shapes of the input data, what assumptions should be used, and what target variables should be predicted. For example, user might write the informal text shown at the top left of Fig. 1. Given such a string, and asked to create models in the Stan PPL [5], an LLM could output any of the formal models shown at the bottom of Fig. 1.
Denote the informal user text and formal model as m. Since LLMs are stochastic, we can think of one (along with system prompt, hyperparameters, etc.) as defining a distribution
P m
formal
textual
model
description
Meanwhile, think of a PPL as defining a distribution target variables and data x, conditional on the formal model. That is, we may think of a PPL as defining a distribution
p
target
observed data
formal
variables
model
The core idea of this paper is to glue the above two distributions together to define the joint
p(z,x, = =p(lt) p(2,x/m)
LLM
PPL
Our hypothesis is that the distribution defined in Eq. (3) is good one. We are interested in the posterior over z, conditioning on and I, but integrating out m. It is not hard to show that this
p(z|z, t)
Σ
p(m|x,t) t) p(z(r,m)
(4)
final posterior
posterior posterior of model weight model (PPL)
Here, p(2|a, m) is the posterior for model (as defined by the PPL in Eq. (2)) and
p(m|r; 1) & p(ml()
p(alm)
(5)
posterior model weight
prior model marginal likelihood weight (LLM) of model (PPL)
Note that the final posterior in Eq. (4) can be
instance of Bayesian model averaging
[16, 32], just with a prior over models that's 't's defined using an LLM.

---

Algorithm Theoretical exact LLB algorithm (intractable)
1. Input textual description and data x.
For all possible model strings m:
(a) Compute model probability p(m]t).
I/ using LLM // under PPL.
(b) Compute posterior p(= jr, m) and marginal likelihood p(a) m)).
3. Set w(m) & p(m]))p(==))) m)), where Em ((()) = 1.
4. Return final posterior p(z|ar,)) = =mm polm)p(222z.)
2.1 Varying latent spaces
Note that the user is not expected to specify all latent variables, only the target variables to predict. This is crucial since, as illustrated in Fig. I, different formal models typically have different latent spaces. This is fine. All that's needed is that each model produced by the LLM contains the target variables a specified in the user prompt t. Formally, suppose that model m defines some distribution
p(z, v. u ((m)mm)) mm),
(6)
where (m) varies in meaning (and dimensionality) for different models m. Then Eqs. (1) to (5) are all still correct. The final posterior remains as in Eq. (4), just with p(2, x/m) interpreted as Eq. (6) after marginalizing out u( (m"), The posterior weights in Eq. (5) make no reference to the latent space, and so need no modification.
2.2 Comparison to flat averaging
is useful to contrast the final posterior p(z|0, t) to the result of taking a "flat" average of of models sampled from the LLM, i.e.
Paa Ir, ()
Σ
p(m I)
p(=|ar, m))
model weight (LL.M)
posterior of model (PPL)
Clearly (2a(2 1) not equal to (2 1) as defined by Eq. (4), but can be seen as an "ensemble of posteriors". The difference is that Pfat gives model m weight equal to the prior p(m]t) whereas p(2|2, t) gives weight that is also proportional to the marginal likelihood p(x|m) (See Eq. (5)). Thus, the fundamental difference between the true posterior and "flat" posterior is that the former gives more influence to models that are more consistent with the observed data.
3 Inference
While mathematically simple, the posterior p(z(a,t) in Eq. (4) is computationally difficult. In principle, one might imagine computing in Alg. but this will rarely be practical. One familiar issue is that for given model m, it's usually difficult to compute the posterior p(2|ær, m) or the marginal likelihood p(z|m). In addition, the space of models m very large or possibly infinite.
A more subtle issue is that while we assume one can sample from p(m](), the probability p(m(() (t) is often unavailable. Many commercial LLM providers decline to share this information. Further, it's often beneficial to ask LLMS to "think out loud" before producing a final answer. This too makes p(m|t) intractable, since the same model could appear after many different "thinking" passages.
So, while the posterior in Eq. (4) can be seen as an instance of Bayesian model averaging, existing algorithms for Bayesian model averaging (e.g. reversible jump MCMC [8, $11J) seem difficult to apply, as they involve explicitly iterating over all models and/or using evaluations of p(m|tt) For these reasons, p(2[1r t) appears to represent novel inference challenge.
To approximate the final posterior, this paper will use the recipe Alg. 2. The simplest interpretation of this recipe is as heuristic approximation of Alg. I where the sum over all models is replaced by random sampling, and p(z[r, m) and f(z(m) denote approximations of the posterior and marginal likelihood for given model Sec. will give more formal justification and analysis.
3

---

Algorithm 2 Suggested generic approximate LLB recipe.
1. Input textual description and data x.
For = 1,2,-..,,N
(a) Sample model (x) p(m)t).
I/ uuiin LLM II under PPL.
(b) Approximate posterior p(z|ar, ml()) and marginal likelihood p(x:) (M))).
Set w(()) x p(xlmm!))), where En=1 w(n) 1.
Return final approximate posterior p p(z|z,t) x Em= p222,
3.1 Preliminary observations
We experimented with using various LLMs to generate formal Bayesian models in various PPLs, including Stan [5], NumPyro [26] and PyMC [1]. LLMs seemed better at generating Stan code, perhaps since more Stan code available and included LLM datasets. We also found it was helpful to give LLMs system prompt that instructed them to describe their modeling strategy in words before producing a formal model. Finally, to aid in-context learning, we found was helpful to provide examples of user prompts along with "good" responses.
3.2 Generating models
These experiments used system prompt (Appendix F.1) that told the LLM first generate "thoughts" block that would describe the modeling strategy in words, and then create the formal Stan model. A few additional instructions were used avoid common mistakes and encourage the LLM to create well-normalized distributions.
For in-context learning, provided the LLM with six examples of ostensible inputs, along with high quality outputs (Appendix G). Each input followed the format illustrated in Fig. with PROBLEM block, describing the problem, the target variable(s). Each output followed the instructions in the system prompt, with THOUGHTS block describing the modeling strategy and MODEL block with formal Stan code.
block describing the data, and GOAL block describing
For the problems below, 1024 models were generated using Llama-3.3-700 [14, 21] with 4-bit AWQ quantization [19]. Since many models were generated from the same prompt, continuous batching (parallel inference) greatly increased inference speed. Using a single A100 GPU, generating 1024 models took 10-15 minutes, depending on the problem.
3.3 Inference
Outputs from the LLM were rejected if they did not contain GOAL and MODEL blocks or if the code in the model block did not compile. For models that compiled, 2 chains of Stan's 'un's implementation of NUTS [17] were run for 10,000 iterations each. These samples were then considered as representing the approximate posterior p(2) 22, m²))).
It remains only to estimate the marginal likelihood p(x)m(())) While methods exist to estimate the marginal likelihood from set of posterior samples [23], these are often considered unreliable [22]. Instead, elected to variational bounds. The typical approach is create some variational
distribution q(2) (e.g. a Gaussian) and optimize to maximize the lower-bound
((z))
ELBO(g) = E log
q(z) < log (())m())).
(8)
zr-q
This bound will be loose the true posterior p(2)a, mn)) is not the variational family. For a bighter-bound, pursued importance-weighter bounds [3] of the form
K
p(z ((mm)))
L(g)
E
log
5 log (ogl(m))()))
(9)
K
g(zb)
where z),... z are independent samples from q.

---

Algorithm 3 The variant of Alg. 2 used in the experiments of this paper.
1. Input textual description and data x.
For = 1,2,,-.,,N...
(a) Sample model (x) ~ p(m]t).
I/ using LLM // under PPL.
(b) Draw samples =(,,1), ((n.)) p(2|2,, ml))) using MCMC.
(c) Compute a lower bound L (n) on the marginal likelihood log p(z(m(() using Eq. (9).
3. Set w(()) exp (m), where E w(x) K-
4. Return the samples (z(n.k) where ((,,) is given weight w(x)
remains to choose the distribution q. The typical approach with sepportance-weightee bounds [3, 10, 11] is to explicitly optimize to maximize Eq. (9). However, is difficult make stochastic optimization fully automatic and reliable across thousands heterogenous models. Instead, we exploit the fact that we have already (approximately) sampled from p(= [a, m(*) using MCMC and simply set q to be Gaussian distribution matching the empirical mean and variance of the MCMC samples.
This seemingly-heuristic heuristic choice is justified by the fact that Eq. (9) is asymptotically tight (as K - oo) with an asymptotic rate of 1 K and constant determined by the x2 divergence between q and [20, 11]. From the perspective of alpha-divergences, maximizing ELBO(g) is equivalent to minimizing ki(glpp) p)) Pp(p| q), an "exclusive" divergence. Eq. (9) would be (asymptoticallly tightest if we could minimize the "inclusive" divergence x" (pllg) D2(pl (), but this is difficult to optimize [13]. Matching the mean and covariance of the posterior is equivalent to minimizing kL(glPp) D2 (pllg), and so is pragmatic compromise. Experimentally, this performed at least as well as explicitly optimizing Eq. (9), but was faster and more reliable.
The full inference recipe is summarized as Alg. 3. More details are in Appendix F.
4 Experiments
It is likely that all standard models and datasets included in LLM training data. To avoid the risk that the LLM would simply "remember" human- human-written written models, these experiments use all-new problems and datasets.
Since this paper will presumably also be included in future LLM datasets, any future work this direction should not use these problems for evaluation with any LLM with knowledge cutoff after the date this paper was first made public, namely April 21, 2025.
PROBLEM
("num_ "'num_days":2 d2,, "rain": (1, 0,000,000000000,
I've been recording if it rains each day, with for rain and for no rain. Haybe there's some kind of pattern? Predict if it will next day.
1.0
0.9
DATA
num_ days
int rain // results
so far, in order
02-
GOAL
00- 80
40
40
pizix Puulzix, O
int next
outcome for next day
log (x (l))
Figure 2: The rain problem. Left: Informal user text t. Top right: The given data æ. Bottom center: Estimated marginal likelihoods (m)) and posterior means E|z|a, m(()) for each generated model (x) Markers for the four models Fig. are colored. Bottom right: The final posterior mean, compared to flat average.
5

---

got the coin from the mint, ((4) I'm almost completely sure that it's standard penny.
pCrjx elit
likelinood
glance, it appears (02) standard US examined it closely.
although I haven't
can see that the penny quite bent, (()) so I'm pretty sure the bias is dirferent from though can't tell in what direction.
platac 120g
plate e
0.4
0.s
0.6
0.7
0.8
0.9
bias
Figure 3: The coin problem. Left: Snippets from three different user prompts. Right: Resulting final posteriors, which appear to reflect user intent. Predictions from Pitat are shown as faint dotted lines.
4.1 Rain
In this problem (Fig. 2), a user recorded rained on sequence of days and wishes to predict if will rain on the following day. The data provides two contradictory signals. One the one hand, most days did not have rain. However, adjacent days are correlated and rained on the last three days.
Out of 1024 generated models, 960 (93.8%) compiled and allowed for inference. Four example models are shown Fig. 7 (supplement). Many models treat each day as independent and give predictions near the base rate of 8/22 = 0.364. Others give predictions near the between -ay consistency rate 14/21 0.667. (Both these base rates are shown dotted lines in Fig. 2.) A flat average includes many models of both types and predicts 0.573. However, the final posterior is dominated by single model that models multi-day dependencies and predicts an even-higher value of 0.793. Detailed results are in Appendix B.1.
4.2 Coin
This problem (Fig. 3) tests the impact of different assumptions stated the user text. Three prompts given describing flipping a coin but with different assumptions about the true bias. In all cases, 20 flips were observed, out which 14 were heads.
Out of 1024 generated models, between 848 (82.8%) and 989 (96.6%) compiled and allowed for inference, depending on the prompt variant. For all three variants, many generated models had marginal likelihoods high enough to contribute the final estimated posterior. The final estimated posteriors are compared in Fig. 3, where the different assumptions seem reflected in the final posteriors. Detailed results are in Appendix B.2, Appendix B.3, and Appendix B.4.
4.3 Polling
In this problem (Fig. 4), user describes candidate with fluctuating popularity who is polled by three different pollsters at various times throughout year, and asks to predict the true popularity on each day. See Appendix B.5 for the full prompt.
100
popularity
80
popularity
80
60
40
pizix
pollster in
20
Prac 0)
20
pollster #3
pollster #2
0
100
200 day
true
300
00
100
200
300
Figure The polling problem. Left: Observed data æ. Right: Final estimated posterior, compared to a flat average.

---

Out of 1024 generated models, 568 (55.5%) compiled and allowed for inference. For this problem. the estimated posterior shows little benefit over flat average. single model captured essentially all posterior weight, which may indicate inaccurate inference. Detailed results are in Appendix B.5.
4.4 City Temperature
100
In this problem (Fig. 5), the temperature is given for set random days in a handful of random cities around the world, along with the temperature on subsequent days. The user gives set of test days and asks to predict the temperature on the following days.
50
476
479
+-
F-+
50
100 3-* +--
1-2 477
6-
6- 3-4
50
100 "F
j-+ 4-2
4-4
f-e
50
100
Out of 1024 generated models, 736 (71.9%) compiled and allowed for inference. In this problem, again the final estimated posterior isn't obviously better than flat average. Detailed results are in Appendix B.6.
+-+
#--
8-9
two
t-
f-+ 4-+
50
Test days
pizjx.
truez
Figure 5: Medians and 90% credible intervals for each city and test day for the the city temperature problem.
4.5 Gold
In this problem (Fig. 6). the user describes rod with a varying density of gold. Binary tests have been done at different positions and the goal is to infer the true density. Inference is done with two different datasets, one with 30 observations, and one with 150.
This problem proved quite challenging, with only 1-2% of generated models being syntactically valid and allowing inference to proceed. Thus, 16,384 models were generated, out of which only 192 (1.2%) compiled and allowed for inference with 30 data and 71 (0.4%) with 150 data. The marginal likelihoods result in almost all weight concentrating in just three models in the smaller dataset and just one model with the larger dataset. Using these gives sensible results, and clearly better than "flat" average, at least with 150 data. Detailed results are in Appendix B.7 and Appendix B.8.
Finally, note that some of the above results may be impacted by the fact that truncated distributions are created through the use of constraints in Stan [5], densities are not re-normalized to reflect the truncation. So, for example, if a model includes standard half-normal, computed log-densities for that model will be lower they should be by a factor of log s 0.693 nats. Unfortunately, this seems difficult to address automatically. This "safe" in the sense that it can decrease but not increase
1.00
gold density
0.75
0.50
true density
0.25
gold detacted
t2xx under pizík, O
not detected
0.00
n
mi
+zix under Puulz|x,
0.0
0.
10
2.0
2.5 0.0
0.5
1.0
1.5
2.
25
position
position
gold density
0.75
0.50
truc density gold detected
0.25
not detected
ELpx under PiZTA. O
0.00
WiN HB*))
12ix. under PuulzZj,
0.0
05
1.
15
2,0
2.5 0.0
0.5
1.0
1.5
2.0
25
position
position
Figure 6: The Gold problem. Left column: Observed data. Right column: Final estimated posterior, compared to a flat average. Top row: Results with 30 observations. Bottom row: 150 observations.
7

---

estimated marginal likelihoods, but it means that some "good" models may have lower weights then they deserve all the above experiments.
Theory
This section seeks to understand how accurately Alg. 2 will approximate the true posterior, in particular how the accuracy depends on the number of samples N, and the quality of the posterior and marginal likelihood approximations, and the LLM-defined prior over models p(m|t).
5.1 Importance sampling in model space
Recall the final posterior p(z|,r t) = Em p(m]z, t))(2 m) from Eq. (4). To cope with the large model space, one idea would be to perform importance sampling to draw m p(m]t) and reweight p(z|æ, m) by the ratio p(m]a jr, t)/p(m/t). Unfortunately, this is infeasible, since the normalizing constant for ((m] in Eq. (5) is intractable and the probabilities p(m(() are unavailable (Sec. 3).
An alternative is self-normalized importance sampling (SNIS)-to draw set of models from p(m]() and give each a weight proportional to p(m|a, /mmmm)), but normalized to sum to one. It can be shown (Appendix C.I) that this estimator simplifies into
p(xlm((m)
p(z|r, 2
Σ (m) p(z|a, (*),, (")
(10)
Ew p(xlm((x))
where m(1), (N) ptm)t). This is equivalent to Alg. 4 (Appendix D).
This will rarely be practical, since the posteriors p(z) x. m()) and marginal likelihoods (()
must be approximated. But is useful to analyze it to isolate error due purely to SNIS. Suppose we are interested in the expectation of some function f' with respect to posterior, i.e. lr = R2 f()). p(z(z, f) is the approximation in (Eq. (10)), one can estimate μ with
N
AN
E f(z)
f(=)=
= Σ (n)
E
[plz|a, ml]))]
(1I)
p(ale,))
p(z)z.m())
n=1
Standard results (Appendix C.2) show that juN is asymptotically unbiased with asymptotic variance
VjūN VN
1
E p(m)r)
p(mmlz t)
(g(m) - u1² p(mll)
N
Further, if (g(m) - 111 S 8 (as would be true, for example, if [f(z) - 14| S 8), then,
8² (1 (1+x² (p(mlz, l)l((m)))))
Vx S N
(12)
where x" (pllg) again denotes the chi chi-squared divergence. Eq. (12) suggests that for Alg. 4 to achieve a given accuracy, N must proportional to x² (p(m) ())m)( This is an "inclusive" divergence, meaning that if p(m)t) is small but p(m| t) is large for some model this greatly increases the divergence, while the reverse situation only causes a modest increase.
This suggests that it is important that the LLM-defined prior p(m[t) be relatively broad: The accuracy of the estimated posterior harmed more by the prior being "overconfident" than being "underconfidentt. In particular, since p(m]2,, () only varies from p(m[() by factor of p(z) (and normalization), is important that p(m(() should not give vanishingly small probability to models m where the marginal likelihood p(æ: m)) is large.
5.2
Approximate inference
This section considers the impact of approximating the marginal likelihood and the posterior for given model. Suppose temporarily that for each m, some approximation q(2|2, m) = p(2|2, m) is known. Define the joint approximation
g(2, m//x = g(m)))g(21z,,))
(13)

---

where g(m[2) remains to be chosen. In principle, we would like to optimize ((m(a) to make the marginal ((2(2) close to p(z|æ, t). but this seems difficult, since even evaluating p(z|z, t) is intractable (Eq. (4)). Instead we propose to choose ((m) a:) to minimize the joint divergence
kL(g(z, m m(r) p(z. m)z,t))
(14)
which (by the chain rule KL-divergence [9, Thm 2.5.3)) an upper bound on the divergence from (<(22) to p(2|2, t). The following result gives the the optimal g(m]2z). (Proof in Appendix E.1)
Theorem 1. Suppose p(z, M2.,2,ml0 and m)) are fixed. Then kL(g(z, m|xx)]pp(z,n) m m|x, ()) is
minimized by
q(ma)) p(ml()((<(m) (exxxxxxxxxxxxxexe (op(-Ku(az)z; m m)llp(((,,)) m)))
(15)
with resulting joint divergence of
KL(g(z, m/) m|x, t)) = log
exp (xp_-L(((z)2,, m)) m)))
(16)
p(m]z,t)
Informally, the optimal ((m)a)) in Eq. (15) can be seen taking p(mlar, t) x p(m|()p(x))) from Eq. (5) and down- weighting models m where g(>) m) is worse approximation of p(z[r, m). Eq. (16) is less intuitive. As KL-divergence is non-negative, the inner expectation is between 0 and 1, with larger values giving a smaller joint divergence. Thus, inaccuracy in ((aaæ, (,, m) "hurts" more when p(m/ar, i) is larger. This same intuition follows from noting that Eq. (16) can be upper-bounded by Ep(m)=,)) (L(g(z)z,, m)llp(z)ar, m)) (See Thm. in Appendix C.5).
Thm. 1 is not immediately useful, since the marginal likelihood p(a[mm) and KI. divergence for given m are typically intractable. However, recall the ELBO decomposition
()//()) = (L + (g(2ar, m) [[](z] m)),
(17)
where the (tractable) evidence bound (ELBO) (now written as function of the model m rather than the distribution q as in Eq. (8)) is
ELBO(m) = E log p(z, x/²) (z|r,, m)
(18)
(z)ar,m)
Variational inference (VT) algorithms maximize the ELBO which is equivalent (since log p(a (m) fixed) to minimizing the diverregee Using this decomposition, Thm. I can be given in the following alternate form. (Proof in Appendix F.2.)
Corollary 2. Under the same assumptions as Thm. 1, the optimal q(m]ax)) can also be written as
((m/x) p(m)t) (ELBO(m)).
(19)
with a resulting joint divergence of kL(g(z, m|x))lp(z,m m|r, 1)) = logp(z|() log exp (ELBO(m)). m(m(t)
(20)
Note that Eq. (19) was previously shown by [18, Sec. 2.2] and Ohn & Lin [24, Thm 2.1].
Taken literally, Thm. 2 suggests that to find the joint approximation q(2, m|x) closest to p(2, m|zz t), one should loop over all models m, independently do variational inference on each, and then average the variational distributions using the weights from Eq. (19). Such algorithm is shown as Alg. 6 (Appendix D), but is not practical since again the space of models m large and the LLM probabilities p(m[t) are unavailable.
However, one can also incorporate the SNIS ideas from Sec. 5.1. Suppose we would like to estimate (<2[2)) = {m))) (2)22, m), where q(2|2, m) is approximation of p(z(x, m) and q(m[x) are the weights from Eq. (19). If we use proposal distribution p(m](), then it's easy to show (Appendix C.3) that the self-normalized weights are proportional to (x)(F..BO(m)). This results in Alg. This an instance of the recipe in Alg. 2, since ELBO(m) is bower-buund on log p(z

---

5.3 Algorithmic variants
It can be beneficial to use approximating distributions (22)2r, m.) where only lower-bound on the ELBO is available. One example of this would be to take advantage of Monte Carlo VI methods like importance-weighted VI [3. 12, 2. 10, 71. Another example would be when one believes MCMC can efficiently sample from p( [z, m.) If one thinks of those samples as representing an approximating distribution {2]2 m), one might believe it is essentially exact. But even if this is true, estimating the marginal likelihood from samples is famously difficult [23, 22], so one might bound the marginal likelihood using a different method (e.g. variational inference).
To justify this, imagine that one first finds a variational approximation g(z[z, m) for each model m, and sets the model weights ((m[2z) using Eq. (19) Then, imagine replacing each distribution p(z|ar, m) with one with a higher BLBO (lower KI--divergence) while leaving q ((m[z) fixed. By the chain rule of KL-divergence, kL(g(z, m| ll)l m) x t)) KL(q(m) (x)) [|]|[mm Jr, t)) + (L(g(z)m, (2)m()) x)). Improving (=||m, x) will decrease the second KL-divergence on the right while leaving the first divergence unchanged, meaning the joint divergence can only decrease.
Appendix C.4 gives a generalized version of the above results when the model weights g (m))t) are com- puted based on inexact ELBO values. A particularly interesting special case is when the distributions p(z|2, m) can be computed exactly (e.g. using MCMC) but the weights ( m) are computed based some bound on the marginal likelihood. (See Alg. in Appendix D.) Then the joint divergence result from Appendix C.4 reduces to (L(g(z, m m|x))]] |p))z m| m|ar, t)) log Ep(m|z.() exp(-8(m) + 8), where 8(m) is the gap between log p(=|m)) (m) and the lower bound used when computing q(m|ax) and ö (m)x)) d(m), This result has somewhat similar intuition as Eq. (16): Constant errors have no effect, and errors on models where p(m| ar, 1) very is small have little effect. What really matters if
g(m) varies among models where p(m|ar, f) is large.
Discussion and limitations
This paper proposed the L.L. scheme for defining a posterior from natural language (Sec. 2), sug- gested broad approximate inference strategy for it (Sec. 3), tested it experimentally (Sec. 4) and analyzed the error of the inference strategy theoretically (Sec. 5).
One direction for future work would be to investigate better ways of using LLMs to define distributions better distributions p(ml(t. There are many obvious options, such as
over models different system prompts, more/better examples, different formatting of user inputs, or training (or fine-tuning) an LLM specifically for this task.
Another direction is better inference methods. The suggested recipe involves doing inference independently each sampled model. While fully parallelizable, this is expensive, and practical only on relatively small problems. We intend no claim of optimality. When doing inference on many models for the same task, many models are often quite similar, suggesting it might be possible to share work between models. It might also be possible to use constrained generation to guarantee that only valid models were created, eliminating the overhead of generating and then rejecting some syntactically invalid models.
6.1 Related work
In terms of usage of LLMs for Bayesian modeling, Wong et al. [33] suggest system where an LLM translates natural language text into a PPL. This is related to our approach, but assumes that "facts" as provided the LLM as text, rather than being described in text and then provided to inference as numerical data. Capstick et al. [4] consider using an LLM to take natural language description and predict the parameters of Gaussian priors (to be used in concert with known likelihood models). Selby et al. [29] also consider a similar idea, and compare the results from several LLMs to priors elicited from human experts Stefan et al. [30]. Requeima et al. [28] suggest a scheme to make probabilistic prediction from natural language plus data. Other work considers the use of L.MMs for time-series predictions [15, 34] or as regressors [31]. Choi et al. [6] suggest a general framework for using natural language information to inform learning procedures via an LLM.
10

---

Acknowledgement
This material is based upon work supported in part by the National Science Foundation under Grant
No. 2045900.
References
[1] Abril-Pla, O., Andreani, V., Carroll, C., Dong, L.,, Fonnesbeck, C. J., Kochurov, M., Kumar, R., Lao, J., Luhmann, C. C., Martin, O. A., Osthege, M., Vieira, R., Wiecki, T., and Zinkov, R. PyMC: modern, and comprehensive probabilistic programming framework in Python. PeerJ
Computer Science, 9:e1516, September 2023.
[2] Bachman, P. and Precup, D. Training Deep Generative Models: Variations Theme. In NIPS Workshop: Advances in Approximate Bayesian Inference, 2015.
[3] Burda, Y., Grosse, R., and Salakhutdinov, R. Importance Weighted Autoencoders. In ICLR, 2015.
[4] Capstick, A., Krishnan, R. G., and Barnaghi, P. Using Large Language Models for Expert Prior Elicitation in Predictive Modelling, December 2024.
[5] Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. Stan: A Probabilistic Programming Language. Journal of Statistical Software, 76(1), 2017.
[6] Choi, K., Cundy, C., Srivastava, S., and Ermon, S. LMPriors: Pre- PraiTre Language Models as Task-Specific Priors. In NeurlPS 2022 Foundation Models for Decision Making Workshop, November 2022.
[7] Christian Andersson Naesseth. Machine Learning Using Approximate Inference: Variational and Sequential Monte Carlo Methods. PhD thesis, Linköping University, 2018
[8] Christian Robert. Monte Carlo Statistical Methods. 2004.
[9] Cover, T. M. and Thomas, J. A. Elements of Information Theory. Wiley-Interscience, Hoboken, N.J, 2nd ed edition, 2006. ISBN
[10] Cremer, C., Morris, Q., and Duvenaud, D. Reinterpreting Importance-Weigee Autoencoders. arXiv:1704.0291 [stat], April 2017.
[11] Domke, J. and Sheldon, D. Importance Weighting and Variational Inference. In NewrIPS, 2018.
[12] Domke, and Sheldon, D. Divide and Couple: Using Monte Carlo Variational Objectives for Posterior Approximation. In NeurlPS, pp- 11, 2019.
[13] Geffner, T. and Domke, J. On the difficulty of unbiased alpha divergence minimization. In ICMI., pp. 3650-3659, July 2021.
[14] Grattafiori, A., Dubey, A., Jauhri, A., Pandey, A., Kadian, A., Al-Dahle, A., Letman, A., Mathur, A., Schelten, A.. Vaughan, A., Fan, A., Goyal, A., Hartshorn, A., Yang. A., Mitra, A., Sravankumar, A., Korenev, A., Hinsvark, A., Rao. A., Zhang, A., Rodriguez, A., Gregerson, A., Spataru, A., Roziere, B., Biron, B., Tang, B., Chern, B., Caucheteux, C., Nayak, C., Bi, C., Marra, C., McConnell, C., Keller, C., Touret, C., Wu, C., Wong, C., Ferrer, C. C., Nikolaidis, C., Allonsius. D., Song, D., Pintz, D., Livshits, D., Wyatt, D., Esiobu, D., Choudhary, D., Mahajan, D., Garcia-Olano, D., Perino, D.. Hupkes, D.,, Lakomkin, E., AlBadawy, E., Lobanova, E., Dinan, E., Smith, E. M., Radenovic, F., Guzmán, F, Zhang, F., Synnaeve, G., Lee, G., Anderson. G. Thattai, G., Nail, G., Mialon, G., Pang. G., Cucurell, G., Nguyen, H., Korevaar, H., Xu, H., Touvron, H., Zarov, I., Ibarra, I. A., Kloumann, Misra, I., Evtimov, Zhang J., Copet, J.. Lee, J., Geffert. J., Vranes. J., Park, J., Mahadeokar, J., Shah, J., van der Linde, J., Billock, J., Hong, J., Lee, J., Fu, J., Chi, J., Huang, J., Liu, J., Wang, J., Yu, J., Bitton, J., Spisak, J., Park, J., Rocca, J., Johnstun, J., Saxe, J., Jia, J., Alwala, K. V., Prasad, K., Upasani, K., Plawiak, K., Li,
K., Heafield, K., Stone, K., El-Arini, K., Iyer, K., Malik, K., Chiu, K., Bhalla, K., Lakhotia,
K., Rantala- Yeary, L., van der Maaten, L.. Chen, L., Tan, L.,, Jenkins, L.. Martin, L., Madaan,
11

---

L.. Malo. L., Blecher, L., Landzaat, L... de Oliveira, L., Muzzi, M., Pasupuleti, M., Singh, M., Paluri, M., Kardas, M., Tsimpoukelli, M.. Oldham, M., Rita, M., Pavlova, M., Kambadur, M.. Lewis, M., Si, M., Singh, M. K., Hassan, M., Goyal, N., Torabi, N., Bashlykov, N., Bogoychev, N., Chatterji, N., Zhang. N., Duchenne, O., Çelebi, O., Alrassy, P., Zhang, P., Li, P., Vasic, P., Weng. P. Bhargava, P., Dubal, P., Krishnan, P., Koura, P. S., Xu, P., He, Q.. Dong. Q. Srinivasan, R., Ganapathy, R., Calderer, R., Cabral, R. S., Stojnic, R., Raileanu, R., Maheswari, R., Girdhar, R., Patel, R., Sauvestre, R., Polidoro, R., Sumbaly, R., Taylor, R., Silva, R., Hou, R., Wang. R., Hosseini, S., Chennabasappa. S., Singh, S., Bell, S., Kim, S. S., Edunov, S., Nie, S., Narang, S., Raparthy, S., Shen, S., Wan, S.. Bhosale, S., Zhang. S., Vandenhende, S., Batra, S., Whitman, S., Sootla, S., Collot, S., Gururangan, S., Borodinsky, S., Herman, T., Fowler, T., Sheasha, T., Georgiou, T., Scialom, T., Speckbacher, T., Mihaylov, T., Xiao, T., Karn. U.. Goswami, V., Gupta, V., Ramanathan, V., Kerkez, V., Gonguet, V., Do, V., Vogeti, V., Albiero, V., Petrovic, V., Chu, W., Xiong, W., Fu, W., Meers, W., Martinet, X., Wang, X., Wang. X., Tan, X. E.,, Xia, X., Xie, Xi, X., Wang. X., Goldschlag, Y., Gaur, Y., Babaei, Y., Wen, Y., Song, Y., Zhang, Y., Li, Y., Mao, Y. Coudert, Z. D., Yan, Z., Chen, Z.,, Papakipos, Z.., Singh, A., Srivastava, A., Jain, A., Kelsey, A., Shajnfeld, A., Gangidi, A., Victoria, A., Goldstand, A., Menon, A., Sharma, A., Boesenberg, A., Baevski, A., Feinstein, A., Kallet, A., Sangani, A., Teo, A., Yunus, A., Lupu, A., Alvarado, A., Caples, A., Gu, A., Ho, A., Poulton, A., Ryan, A., Ramchandani, A., Dong, A., Franco, A., Goyal, A., Saraf, A., Chowdhury, A., Gabriel, A., Bharambe, A., Eisenman, A., Yazdan, A., James, B., Maurer, B., Leonhardi, B., Huang, B., Loyd, B., Paola, B. D., Paranjape, B., Liu, B., Wu, B.. Ni, B., Hancock, B., Wasti, B., Spence. B., Stojkovic. B., Gamido, B., Montalvo, B., Parker, C., Burton, C., Mejia, C.,, Liu, C., Wang, C., Kim, C., Zhou, C., Hu, C., Chu, C.-H., Cai, C., Tindal, C., Feichtenhofer, C., Gao, C., Civin, D., Beaty, D. Kreymer, D., Li, D., Adkins, D., Xu, D., Testuggine, D., David, D., Parikh, D. Liskovich, D., Foss, D., Wang. D., Le. D., Holland, D., Dowling, E., Jamil, E., Montgomery, E.,, Presani, E. Hahn, E., Wood, E., Le. E., Brinkman, E., Arcaute, E., Dunbar, E., Smothers, E., Sun, F. Kreuk, F. Tian, F., Kokkinos, F,, Ozgenel, F., Caggioni, F., Kanayet, F., Seide, F., Florez, G. Schwarz. G., Badeer, G., Swee. G., Halpemn G., Herman, G., Sizov, G., Guangyi, Zhang, Lakshminarayanan, G., Inan, H., Shojanazeri. H., Zou, H., Wang, H., Zha, H., Habeeb, H., Rudolph, H., Suk, H., Aspegren. H., Goldman, H., Zhan, H., Damlaj, I., Molybog. I., Tufanov, I. Leontiadis, I., Veliche, I.-E., Gat, I., Weissman, J.. Geboski, J., Kohli, J., Lam. J., Asher, J., Gaya, J., Marcus, J., Tang. J., Chan, J., Zhen, J., Reizenstein, J., Teboul, J., Zhong, J., Jin, J. Yang. J. Cummings, J., Carvill, Shepard, J., McPhie, J., Torres, J., Ginsburg, J., Wang. J., Wu, K., U, K. H., Saxena, K., Khandelwal, K., Zand, K., Matosich, K., Veeraraghavan, K., Michelena, K., Li, K., Jagadeesh, K., Huang. K., Chawla, K., Huang. K., Chen, L.,, Garg, L., A, L.,.Silva L.,, Bell, L., Zhang, L.. Guo, L., Yu, L., Moshkovich, L., Wehrstedt, L., Khabsa, Avalani, M., Bhatt, M., Mankus, M., Hasson, M., Lennie, M., Reso, M., Groshev, M., Naumov, M., Lathi, M., Keneally, M., Liu, M., Seltzer, M. L., Valko, M., Restrepo, M., Patel, M., Vyatskov, M., Samvelyan. M. Clark, M., Macey, M., Wang, M.. Hermoso, M. J., Metanat, M., Rastegari, M., Bansal. M.. Santhanam, N., Parks, N., White, N., Bawa, N., Singhal. N.. Egebo. N., Usunier, N., Mchta, N., Laptev, N. P., Dong, N., Cheng, N., Chernoguz, O.,, Hart, O., Salpekar, O., Kalinli, O., Kent, P., Parekh, P., Saab, P., Balaji, P., Rittner, P., Bontrager, P., Roux, P., Dollar, P., Zvyagina, P., Ratanchandani, P., Yuvraj, P., Liang, Q., Alao, R., Rodriguez, R., Ayub, R. Murthy, R., Nayani, R. Mitra, R., Parthasarathy, R., Li, R., Hogan, R., Battey, R., Wang. R., Howes, R., Rinott, R., Mehta, S., Siby, S., Bondu, S. J., Datta, S., Chugh, S., Hunt, S., Dhillon, S., Sidorov, S., Pan, S., Mahajan, S., Verma, S., Yamamoto, S., Ramaswamy, S., Lindsay, S., Lindsay, S., Feng. S., Lin, S., Zha, S. C., Patil, S., Shankar, S.. Zhang, S., Zhang, S., Wang, S., Agarwal, S., Sajuyiggee, S., Chintala, S., Max, S., Chen, S., Kehoe, S., Satterfield, s., Govindaprasad. S., Gupta, S., Deng. S., Cho. S., Virk. S., Subramanian, S., Choudhury, S., Goldman. S., Remez, T., Glaser, T., Best, T., Koehler, T., Robinson, T., Li, T., Zhang, T., Matthews, T., Chou, T., Shaked, T.. Vontimitta, V., Ajayi. V., Montanez, V. Mohan, V., Kumar, V. S., Mangla, V., Ionescu, V., Poenaru, V., Mihailescu, V. T., Ivanov, V., Li, W., Wang, W., Jiang, W., Bouaziz, W., Constable, W., Tang, X., Wu, X., Wang, X., Wu, X., Gao, X., Kleinman, Y., Chen, Y., Hu, Y., Jia, Y., Qi, Y., Li, Y., Zhang, Y., Zhang, Y., Adi, Y., Nam, Y., Yu, Wang, Zhao, Y., Hao, Y., Qian, Y., Li, Y., He, Y., Rait, Z., Vitoo Z., Rosnbrick, Z., Wen, Z., Yang, Z., Zhao, Z., Ma, Z., and and others. The Llama 3 Herd of Models, November 2024.
[15) Gruver, N., Finzi, M., Qiu, S., and Wilson, A. G. Large Language Models Are Zero-Shot Time Series Forecasters. In NeurlPS, 2023.
12

---

[16] Hoeting, J. A., Madigan, D., Raftery, A. E., and Volinsky, C. T. Bayesian Model Averaging: A Tutorial. Statistical Science, 14(4):382-401, 1999.
[17] Hoffman, M. D. and Gelman, A. The No-U-turn sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15():11993-163, 2014.
[18] Kejzlar, V., Bhattacharya, S., Son, M., and Maiti, T. Black Box Variational Bayesian Model Averaging. The American Statistician, 77(1):85-995 Jamuary 2023.
[19] Lin, J., Tang, J., Tang, H., Yang. S., Chen, W.-M., Wang, W.-C., Xiao, G., Dang, X., Gan, C., and Han, S. AWQ: Activation-awan weight quantization for LLM compression and acceleration. In MLSys, 2024.
[20] Maddison, C. J., Lawson, J., Tucker, G., Heess, N., Norouzi, M., Mnih, A., Doucet, A., and Teh, Y. Filtering Variational Objectives. In Guyon, I., Luxburg. U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan. S., and Garnett, R. (eds.). NeurlPS, 2017.
[21] Meta. Llama 3.3 Model Cards and Prompt formats. http://www.llaaa..ccc/o////meees poddproomp--frmats/llie
_3, 2024.
[22] Neal, The Harmonic Mean of the Likelihood: Worst Monte Carlo Method Ever. https://radfordneall wwwp://fooonaaaa.ooooooo..... worst-monte-carlos mattiie---coo-mmmod--eeee method-ever,, ever/, August 2008.
[23] Newton, M. A. and Raftery, A. E. Approximate Bayesian Inference with the Weighted Like- lihood Bootstrap. Journal of the Royal Statistical Society: Series (Methodological), 56(1):
3-26, January 1994.
[24] Ohn, I. and Lin, L. Adaptive variational Bayes: Optimality, computation and applications. The Annals of Statisties, 52(1):335-363, February 2024.
[25] Owen, A. Monte Carlo Theory, Methods and Examples. 2013.
[26] Phan, D., Pradhan, N., and Jankowiak, M. Composable Effects for Flexible and Accelerated Probabilistic Programming in NumPyro. arXiv: [cs, stat], December 2019.
[27] Price.
code so I don't have to. htt.cc//mmii..i.///////////aaaae have-to/, April 2023.
P
ChatGPT4
writes Stan
[28] Requeima, J., Bronskill, J. Choi, D., Turner, R. E., and Duvenaud, D. LLM Processes:
Numerical Predictive Distributions Conditioned on Natural Language. In NeurlPS, November 2024.
[29] Selby, D. A.. Spriestersbach, K., Iwashita, Y., Bappert, D., Warrier, A., Mukherjee, S., Asim, M. N., Kise, K., and Vollmer, J. Had enough of experts? Elicitation and evaluation of Bayesian priors from large language models. In NeurlPS 2024 Workshop on Bayesian Decision-making and Uncertainty, October 2024.
[30] Stefan, A. M., Katsimpokis, D., Gronau, Q. F. and Wagenmakers, E.-J. Expert agreement in prior elicitation and its effects on Bayesian inference. Psychonomic Bulletin & Review, 29(5): 1776-1794, October 2022.
[31] Vacareanu, R., Negru, V. A., Suciu, V., and Surdeanu, M. From Words to Numbers: Your Large Language Model Is Secretly A Capable Regressor When Given In-Context Examples. In First Conference on Language Modeling, August 2024.
[32] Wasserman, L. Bayesian Model Selection and Model Averaging. Journal of Matheratical Psychology, 44(1): 44(1):92-107, March 2000.
[33] Wong, L. Grand, G. Lew, K., Goodman, N. D., Mansinghka, V. K., Andreas, J., and Tenenbaum, J. B. From Word Models to World Models: Translating from Natural Language to the Probabilistic Language of Thought, June 2023.
[34] Xue, H. and Salim, F. PromptCast: A New Prompt-Base Learning Paradigm for Time Series Forecasting. IEEE Transactions on Knowledge and Data Engineering, 36(11):6855 -6664, November 2024.
13

---

A
Acronyms and notation
Aeronyms
LLM: Large Language Model
PPL: Prohabilistic Programming Language
BMA: Bayesian Model Averaging
SNIS: Self-normalized importance sampling
VI: Variational inference
Notation
t: Plain-language description of inference problem.
Formal model definition (in PPL)
Input data
Query variables latent variables
n: Index for sampled model
m (n) rith sampled model
w()); weight given to nth sampled model
Terminology
Single-model posterior: p(2|2, m)
Posterior model weight: p(mlar t)
p(m 4)p(arlm))
p(mlt)p(armm
P(m (Dp(zlm))
Final posterior: t)
14

---

B Experimental details
This section gives detailed results for the experiments. For each problem, describe the prompt and data, show example generated models, attempt to visualize the posterior of each model, show marginal likelihoods and weights, and the final estimated posterior, compared against "flat" average as in Eq. (7).
B.1 Rain
In this problem, the user describes recording if it rained not on a series days, and wishes to predict if will rain on the next day. Mcanwhile, the given data æ provides two contradictory signals: On the one hand, most days did not have rain. However, adjacent days are correlated and it rained on the final day.
PROBLEM
=22,
day, rein end 0 for no rein, Maybe there's some kind of pattern Predict rain day.
recording
DATA
int days
arrey days] days]
results far,
int next // outcore for next day
Out of 1024 models generated, 742 (72.5%) compiled and allowed for inference. Fig. shows four example models. Fig. shows the estimated marginal likelihoods and posteriors for each of the valid models. Finally, Fig. shows the estimated marginal likelihoods and corresponding weights for the top-scoring 50 valid models, as well the final estimated posterior p(=lar, t) compared to "flat" (=|ar,,
average Piat (2(ar,
In this case, there are common "clusters" of models. One consists of models similar m (247) (Fig. 7c) that treat each day as independent and give predictions near the base rate of & 22 ~ .364. A
second consists of models similar to (463) (Fig. 7b) that model the probability that rain is followed by no-rain or vice-versa and give predictions near the base rate 21 s 0.667 for between-day consistency. Models from the first cluster tend to have lower-marginal likelihoods than models from the second cluster. However, the top-scoring model m (742) (Fig. 7a) has still-higher marginal likelihood, and so receives essentially all weight the estimated final posterior p( p(2(1, ().
15

---

int nun.doys days;
int nun.days;
array (runy days] Intelower=o, upper=ls raln;
array days] Intelowerso, upper=l> rali;
paraneters {
paraneters {
realstover-.19. 18, upper=18> Logst, initial;
realstower=e, upper=l>
real outoregrassive_cmee confficients array(nun.days-
real-touerso, upper=lx niven
noise;
model
trans forsed paraneters(
purain given beta(),
arrey(nun.dayy real logit_prob;
p_rain. given..ennnnnniinn beta(1,1);
logit. prob]]] logit initial,
rain[]] bernoull1 // unifors prior for first day [ays)
tor(1 d
1ogit_probff1 proit_probff-11 autoregressive coe moise(i-1];
If (rainfi-11=111
rein(i] bernouttiipp revnn.
else { rain[i] bernoutli rwi.. given...eeenneeeee.....
nodel
logit instial uniform(-1e, 18);
autoregresstvve.oo coneffictttttttfiiiienntt 18,11
moise cauchy(0,101;
redd quantities {
1:mun daysi(
rainfi] bernoulli, toggt
if [rainfmun days]mm1I
next bernoulli.rng gew.:eeerreinn....eeeeeeeee
gonerated quantities [
next sernoulti_rng(p.. ratiannnnnnnn.eeeeeee.ee ravingiven, given, no. rain);
real logit_neet probit.pooffumm.c days] sutoregressive.oom coef ficieit cauchy..rofe 10) Logit. matlogit next);
(b) m (463) typical of one "cluster" of models.
(a) (742), top-scoring model.
data
data
nun_deys;
nunudeys:
array[nun. days] int=looerso, uppersle rain;
ar anyay days] int-looerss) uppersle rain; // resultss so
paraneters { reatstover=o, uppersi> R1
paraneters { reatstower=0, upper=1> real<loverr0, upper=ls pll;
sodel uniform(0.1) for (1 in 1:num daysi( rain[i] bernoulli(p);
model beta(29, 291:
pll beta(28, 281;
rain]]] bernoutli // no information about the tirst
generated quantities int next rnggoulli.nggp))
for (1 in denux
rain(i] bernolliitrein(i1] pil po13;
(c) m (247) typical of a second cluster.
generated quantities { int next: iffrainfrun days]1f next bernoull..ng(p11);
next bernoulls.rng(p1])
(d) m (6) low scoring model.
Figure 7: Example models for the rain problem, sorted by estimated marginal likelihoods.
16

---

1.0
0.8
E 0.6
EIZ
0.4
0.2
0.0-
3
100
300
300
count
count
30
10
3
1-
-80
60
-40
-20
20
log pm (n))
Figure 8: Top left: Estimated log marginal likelihoods log logp(zlm(()) and the estimated probability of rain on the next day E[2|c, mee],,,("))) (n]]] for each of the valid models for the rain problem. Markers colored for four example models from Fig. Top right: Histogram for estimated probability of 8 rain. Bottom left: Histogram for estimated marginal likelihoods. The dotted lines at 22 = 0.364 and
14 = .667 show the base rates for rain and between-day consistency in rain, respectively.
17

---

1.00
well Meignt M (m (m)
log 1 (m)) (n)
20
0.75
0.50
0
4009
0.25
0.00
700
710
720
730
740
model number n
1.0
0.8
X
0.6
0.4
0.2
0.0
plz|x, t) Pnat(z|x, t)
Figure 9: Left: Estimated marginall likelihoods log p(z ("+" symbols, left axis) and weights (m) (bars, right axis) for the top 50 models n for the rain problem. Right: The final estimated posterior, compared "flat" average. The dotted lines 22 A 0.364 and rates for rain and between-day consistency in rain, respectively.
2T R .667 show the base
18

---

B.2 Coin (standard)
In this problem, the user has flipped a coin several times and would like to predict the true bias. In this first version, the user strongly implies that it is "standard" coin with bias near 0.5, The observed data are 20 flips, out of which 14 were heads.
PROBLEM
Plips" ; 20, "nun.heeds"
I have coin. I've Nlipped It bunch of Lines. wondering what the true bias of the coin
got the fron corpletely sure that it's standard US penny-
alnost
DATA
int flips: nunber of times coin flipped int nun, heads; // nurber of times heads was observed in flips
real bias; true bies of the coin
Out of 1024 models generated, 989 (96.6%) compiled and allowed for inference. Fig. 10 shows four models. Fig. 11 shows the estimated marginal likelihoods for each of the valid models. Fig. 12 compares the estimated posteriors for each valid models. Finally, Fig. marginal likelihoods and corresponding weights for the top-scoring 50 valid models, as well as the final estimated posterior p(2|2: t) compared to "flat" average Paat (2| t).
shows the estimated
In this case, while there two high-scoring models similar m (968) (Fig. 10a), they essentially "outvoted" in the posterior by lower-scoring but more numerous models similar to m (938) (Fig. 10b). Thus, the final posterior in Fig. 13 is similar to the posterior from m (938).
data
data
int nun. flips: hoads:
int nun. flips:
paraneters real Logit. bias;
paraneters
realstoverso, upper=l> bias;
transforred paraneters {
sodel
real bies inv. wtw.looii(looitt..o bias);
bias nomal(0.5, 0.01):
binosial (nun.
model {
logit.bias nornal(8, 8.2);
nur. heads oinomtalinun. rlips, blas];
(b) (938) model from the most common "cluster".
generated quentities {
// Nottinng to generate here, just inference on the bias
paraseter
(a) m (968) the highest-scoring model.
data
data int fuip:
int nun. Flips: int nun. heads;
int nun. heads)
parareters { realstover-e, uppar-l> bias;
paraneters { realxtover-o, upper-l> bias;
sodel
sodel
biss beta(509, 4981;
bias beta(19999, 10000);
nun. heads binomial tnun. tlips, bias]:
num. heads binomia ionssaatttuu trun. www.saaattunnttaaaaa rlips, biasl;
(c) m (385)
low-scoring model
(d) m (29)
very low-scoring model
Figure 10: Example models for the coin (standard) problem, sorted by estimated marginal likeli- hoods.
19

---

300-
300
100
count
100-
count
30
30-
10
3-
3-
1-
1-
600000
400000
200000
30.0
27.5
25.0
22.5
20.0
17.5
15.0
log (m)
log (xg((((m(n)
Figure 11: Estimated marginal likelihoods p(z[m for each of the valid models n for the coin (standard) problem. Because of the many order of magnitude, two different ranges are shown. Bars are colored for four example models from Fig. 10. multiple models map to the same bin, all colors are shown stacked.
30
[biaas]
20
10
0.0
0.2
0.8
1.0
bias
Figure 12: The estimated posteriors for each of the valid models for the coin (standard) problem. Posteriors are colored for the four example models from Fig. 10.
-15-
0.05
pizix, 1)
log mm (x)mln) (n))
-16-
-0.4
Prxx
weeon in w w(m
20
max likelihood
-17
-0.33
plblas
-18-
-0002
10-
-19-
-0.01
-20
0..00
model number
bias
("+" symbols, left axis) and weights (ml) (bars, right axis) for the top 50 models n for the coin (standard) problem. There are two models similar to m (968), which each receive a weight of approximately 0.036, while there are several hundred models similar to m (938), which each receive weight of approximately 0.0032. Right: The final estimated posterior, compared to "flat" average. The maximum of 14/20 0.7 is also shown for reference.
Figure 13: Left: Estimated log-marginal likelihoods log
maximum-likelihor likelihood estimator
20

---

B.3 Coin (looks)
In this second version of the coin problem, the user states that the coin "looks" like a standard coin, but implies less confidence. As before, the observed data are 20 flips, out of which 14 were heads.
PROBLEM
("non. Clpps 20, "punsheeds": 14)
I have a coin, I've tlipped it bunch of tines. 1'm blas
wandering
clance, appeers
penny,
though haven't It examined it closely.
coin flipped int wun.heads; // number of tines heads was observed in flips
int
GOAL
real blas; // true blas of the coin
Out 1024 models generated, 983 (96.0%) compiled and allowed for inference. Fig. 14 shows four models, while the results of inference shown in Fig. 15, Fig. 16, and Fig. 17.
In this case, there are seven models similar to the top-scoring model m (960) (Fig. 14a) that each have weight around 0.1, and seven models similar to m (949) (Fig. 14b) with weight of around 0.03. All of these have priors with reasonable amount of uncertainty. There many models similar to m (661) that have prior with low variance, but these have low marginal-likelihon and so get very low weight. The final estimated posterior in Fig. 17 thus resembles mixture of the posteriors of m (960) and m (949).
data
data nun.rLips;
nun. tlpps; int nun. heads;
int nun. heads;
parareters
realstover=0, upper=1> bias;
real realstover=0, upper=1> blas;
model
model normal(B.5, 8.13;
uniform(9, 11:
nux. heads binosial foun. rlips, blas];
nun. heads binosial (nun. flips, blasl;
m (960), the highest-scoring model.
(b) m (949)
a slightly lower-scoring model.
data
data
1nt nun. tips; nun. heads;
int nun. nup.r nun. heads;
paraneters { realstover=o, upper=l> bias;
paraneters { realstover=0, upper=l>
sodel
bias beta(60,600 681:
bias betar638, 668);
heads ninovial (nun. flips, bias]:
binosi (mun. ips, biasl:
(c) m (661) from the most common "cluster" of mod- cls.
(d) m (1) the lowest-scoring model
Figure 14: Example models for the coin (looks) problem, sorted by estimated marginal likelihoods.
21

---

300
300-
100
100-
count
count
30
10
10
3-
1-
800
600
400
200
log (m)
log (x (n)
Figure 15: Estimated marginal likelihoods (1) for each of the valid models n for the coin (looks) problem. Because of the many order of magnitude, two different ranges are shown. Bars are colored for four example models from Fig. 14, If multiple models map to the same bin, all colors are shown stacked.
25-
20-
plbias
10
5-
0.0
0.2
0.4
0.8
1.0
bias z
Figure 16: The estimated posteriors for each of the valid models for the coin (looks) problem. Posteriors are colored for the four example models from Fig. 14.
0.15
-15-
direee fnt w(m (n))
pizix. 1)
log ppx(m (m))
Puxlz|x max likelhhod
0004
0.10
-20
p[bias]
0.05
-30
940
950
960
970
980
0.00
0.0
0.2
0.4
0.6
0.8
1.0
model number
bias
Estimated log-marginal likelihoods 1(z m (*)) ("+" symbols. left axis) and weights (m) (bars, right axis) for the top 50 models n for the coin (looks) problem. Right: The final estimated posterior, compared a "flat" average. The maximum likelioooikeooood estimator of 14/20 0.7 is also shown for reference.
Figure 17: Left:
22

---

B.4 Coin (bent)
In this third version the coin problem, the user states that the coin looks "bent", implying great uncertainty about the true bias. As before, the observed data are 20 flips, out of which 14 were heads.
PROBLEM
Lipss : 28, "nun. heads":
tines. wondering what the true blas of the coin 1s.
coin. tlipped
I that the penny quite bent, preity the hiss i5 different from 0.3, though can't tell in what direction.
int nun. flips: nurber Lines coin wes Rlipped int nun. heads; number of tines heads observed in the tlips
GOAL
real blas; true blas of the coin
Out of 1024 models generated, 848 (82.8%) compiled and allowed for inference. Fig. 18 shows four models, while the results of inference are shown Fig. 19, Fig. 20, and Fig. 21.
In this case, many models contribute substantially to the posterior. While m (811) has slightly lower estimated marginal likelihood than the highest scoring model (m (m2(0))), models similar to it were generated many times, and so wields the most influence in the posterior.
data int nun. flips; nun.heads;
int nun. flips: nun.heads;
paraneters {
paraneters { realclovern0, upper=l> blas;
real blas_Togit:
transformed paranaters
realslover=0, upperml> bies = inv. inv. logit);;
bias beta(0.5 0.51; // broad prior, but with more moss
model
num. heads binomial (num, rlips, blasl;
blas_ uniforet-10, 10); diffuse prior that bistributes fairly evenly over [0.1] after logistic transfor
(b) m (811), typical model
binosialinun. flips,
cuantities
// nothing
(a) m (830) the highest-scoring model
int nun. flips; nun.heads;
int fun_flips; nun.heads]
parareters {
peraneters {
realslover=0, uppersis bias;
real realslover=0, upper=i> bias;
bias beta(2,2); // symnetric around 8.5, Low density
bias binanial riips, bias!:
beta(29,39))
nux, heads (1nn blas];
(d) m (1) the lowest-scoring model
(c) (531) from the most common "cluster" of mod- els
Figure 18: Example models for the coin (bent) problem, sorted by estimated marginal likelihoods.
23

---

300-
100
100
count
count
30
30
10
10
3-
3-
1-
1-
50
40
30
20
20
18
16
14
12
log ((x
log (((mmm))
Figure 19: Estimated marginal likelihoods (()[m(()) for each of the valid models n for the coin (bent) problem. Because of the many order of magnitude, two different ranges are shown. Bars are colored for four example models from Fig. 18, multiple models map to the same bin, all colors are shown stacked.
10
8-
pibiast
6-
0.0
0.2
0.8
1.0
bias
Figure 20: The estimated posteriors for each of the valid models for the coin (looks) problem. Posteriors are colored for the four example models from Fig. 18.
0003
pizix, 1)
-12-
weight w(min) Ini)
log pmxmm (m)))
Prai
(1xmml))
-0.02
max likelihood
-13-
blbia
2-
-0.01
-14
-15
0.00
model number
bias
Figure 21: Left: Estimated log-marginal likelihoods logp(s//m ('+" symbols, left axis) and weights (m) (bars, right axis) the top models n for the coin (bent) problem. Right: The final estimated posteriors, compared to a "flat" average. The maximum-likelihon estimator of 14/20 = 0.7 also shown reference.
24

---

B.5 Polling
In this problem, the describes a politician with fluctuating levels of popularity and different pollsters who take polls at different times. The goal is to predict the true popularity on each day.
PRODLEM
rating lbetween 1881 that changes pollsters that measure the poputarity of the politician. Each pollster does palls on sone different subset of days, Fach polister does polls
polstician.
There
ditrering nosse.
have no idea what their popularity night the beginning of year. Make sure model the variability in the troe popularity, and nake D4TA
nodel distribution the noise levels of different pollsters.
int nuv days; // mamber of days over which poputarity is andeted
nur-pollsters:; nurber pollsters
int polls; // maximur nurber of palls performed by pollster
array (nray ponsssnslls int mon.polls; mumber of polts done by each pollster
array [nun pollsters, palls] int num. palls[t])
what each poll done? [day[i,3] only meaningful
array wonyynunnnuottteeoo pollp
actual polting data (polis]f.]] only meaningful j num_polls[11)
GBAL
array days) popularity; // true popularity of the politician on each day
sample of the observed data is shown below, while the full data plotted in Fig. 22.
(
days": 365,
*num pollsters":
*nex. polls"; 100,
polls": [5, 20, 198).
*day": [1291, 26, .... 1I.
[365, ....
[314, 66, .... 5711.
"polts": [[51, 1158388552505786244, 48.741758128772864, .... 8.91.
248.4419596 144896. 48. 09.82574746044996.
0.01.
(46. 236522244660244 34.54022446996244
.... 22.90157200195961
100
popularity
80
60
40
pollster #1
20
pollster #2
pollster #3
0
100
200
300
day
Figure 22: Observed data for the polling problem.
Out of 1024 models generated, 568 (55.5%) compiled and allowed for inference. Fig. 23 shows three models. Fig. 24 shows the estimated marginal likelihoods for each of the valid models. Fig. 25 compares the estimated posteriors for each valid models. Finally, Fig. 26 shows the estimated marginal likelihoods and corresponding weights for the top-scoring 50 valid models, as well as the final estimated posterior p(z|x, t) compared to "flatt" average (2),(2)
In this case, the highest-scoring model m (568) (Fig. 23a) very "simple" yet has high enough marginal likelihood to take almost all weight in the final posterior in Fig. 26. The LLM did attempt to create many more sophisticated models, e.g. hierarchical models. But these often failed to compile or, as with m (89), simply had lower marginal-likelihood. likelihood.
25

---

data
data (..
int nun, pun.polisters
arraytnun days] real poputarity;
pol.p;
prrey(run__oolstee pollsters] realelovernd pollster.noise:
array[nun pollstars] int nun polls;
array fnun. pollsters, sax polts] int day: arrey(nun, pollsters, rex.polls] real polls;
sodel
// popularity changes tine according simple randon walk
transtormed
popularity[1] norral
idea what the initial popularity is deys
vectorfnun_dayst zeros rep_vector18, nurudays);
paraneters {
popularity[1] normal (pouuaarttyl1- 53: not 100 such change from one day to the next
vector vector[nun.coe days] popularity;
array(mun.pollstere realslowersde
1/ observed polls are noisy versions of the true popularity for (1 in 1:nux_ polsstotsss
sodel
//-olts are modeled by a distribution with some underlying and noise sigre
1:nun. polls(i))
Wwhiich is potlster-dependee popularity(1] uniforn(0, 100); 2
poltsfi.j1 11
noraal hormll (poutarityfoay(i..]l
pollster.noise(Al):
poputarsty(s] normal ( pprralppopuaarrtyttt. 18);
I/ prior noise tevels poltsters for (1 in 1:nux, pullsterrs
for (i in 1:men pollteslls
signa(1] inv garra 20); for Ienux, polls(i]M
pollster. noise(i] cauchy(a, 5);
polls(i,1] norrelppoooaariyylyyyl.... sioms(i));
generated ouantities
array(nun [ays] real poputarity.sarpte poputarity;
(a) m (566)
the highest-scoring model.
(b) (553) lower-scoring model.
data [....
paraneters {
real poputatt.
real-couerdds step-size;
arrey errey[nun, [run.deys1]
steps;
vectoneloweride Inux, pollsters] pollster.
realstover-de pollster, neise_nean;
realstower-de collster_nolse_signn
transforned
array prrayfrun.day) real poputarity; popularity(1) initiel. popelartty for (1 in 2:num days]t
poputarityE1] popularsty(s-1]] stops[t-1];
initial..opularity (50:at(50 // very vegue prior over initist, popularity step size normal(8,5);
for (1 in 1:(num, days-13it
stops[t] nomall1..tep_sst2)) randon walk poputarity
caltster.noise.mae sean (S, 181: // very vague prior over the nean of the notse
pollster.noise sigsa cauchyto,21: very vague
signa
pollster.noise
hwwpppoolsteeeneiee.eeeeoee 1snum pollsters)(
noiseerr.niss.siiee // hierarchical prior for the noise of each pollster
for in 1:nin, pan__palls(t]( pollsfi.j] noreal (popstarityyddyylii JI1,, poltster_noise(ilb)
(c) m (89) a lower-scoring (hierarchical) model.
Figure 23: Example models for the polling problem, sorted by estimated marginal likelihoods. The data block was always identical, so is truncated for space for some models.
26

---

300
100
count
30
10
3-
400
200
200
400
600
800
log log (()(m
Figure 24: Estimated marginal likelihoods p(a ()) for each of the valid models for the polling problem. Bars are colored for four example models from Fig. 23.
100
popularity
80
60
40
20
0 0
100
200 day
300
Figure 25: The estimated posteriors for each of the valid models for the polling problem. Posteriors are colored for the four example models from Fig. 23.
log pgppp// [m|m)) (m)
600
1.00
weight w ((min))
popularity
80
0.75
60
400
-0.50
40
200
0.25
piz|x t)
latent 300
520
5š0
540
sšo
560
0.00
0
100
200
model number
day
Estimated log-marginal likelihoods log p(a m ()) ("+" symbols, left axis) and weights (() (bars, right axis) for the top 50 models n for the polling problem. Right: The final estimated posteriors, compared to "flat" average.
Figure 26: Left:
27

---

B.6 City temperature
In this problem the user describes observing temperatures on neighboring pairs days for a set of different cities. Then, after observing temperatures on set of "test" days for each city, they wish to predict the temperature on the following days.
PROBLEN
There is training dala for e bunch of different citlies on rendon days, along with the Lemgerature for the following data. there
for bunch different days sane Predict tesperature the following days.
Assume different eifies So when guessing the weather the next day, both the typical temperature in the city and the weather the pravinus day witt natter.
different weather. assune weather tends stailar neighboring
All terperatures are in l'hhrenheit..
Be sure to use informative priors for all quantities.
DATA
nor.cities;
included
int trr.days.trsin; // number of training days in each city tays
troining days prraytnun_cition nire dadddaysssannn real thyp_terp_tree train; // temperature on the first day, training set array [nun cities, mun_days, train] real terp..er.ttati // tesperature on following day, training set array cities, test] resl teyp.terp. terp.terettt test; // terperature the first day, test set
city
GOAL
array Inwniccttioos, nun. dww.ddy..tee.... test] toy2.torp.test. // terperature on fallowing day, test set
In the given data, there are 5 cities, with 10 training day pairs, and 10 test days.
"non.cities": "nun. trantit 18,
days. tost": 18,
Mayl_tenp_trafn
[161.0, 74.9, 61.5, 73.2,
78.1
69.8,
55.2,
74.2,
53.9. 75.07.
[73.0, 78.3, 89.2.
84.2.
79.0
81.3.
77.3.
76.0).
179.3, 82.1, 82 76.0.
81.2,
84.9.
79.0,
79.0,
81.1. 85.71.
[01.0, 76.5, 82.0.
70.0.
77.6,
77.1
76.9,
82.2.
03.7.
76.5).
179.8, 78.8, 77.8,
79.6,
76.8,
75.4,
88.4,
75.6,
81.0,
76.2]1.
"dayl_ tomp.
1169.3, 77.5. 74.8.
65.3.
58.5.
87.0.
54.9.
48.5.
59.4,
72.91.
176.5, 4, 84
6.
85.1
81
01.0.
77.4,
71.9.
77.7).
82.2.
79.7].
102.7, 74.7, 8.
77.9
8,
76.2.
81.8,
76.5,
79.1.
02.11.
[78.4, 75.1, 91.9.
75.0.
74.5.
74.2,
00.5.
76.5.
79.0.
77.8]).
tay2_ttep.ttramm
[161.5, 69.9, 62.8,
73.4,
78.8,
68.8,
53.1
74.5.
49.4,
76.21.
[76.4, 78.0. 82.4,
76.9.
80.
82
79.4,
79.1.
78.3.
00.01,
179.5, 83.8,
7..9.
7..0
78.5,
8,
79.0,
81
89.0.
05.5).
180.5,
888.8, 7,
75.3.
1.
74.0,
38
81.6,
77.8,
79.5,
72.811
Out of 1024 models generated, 736 (71.9%) compiled and allowed for inference. Four example models shown in (Fig. 27) and (Fig. 28). A histogram of estimated marginal likelihoods is in (Fig. 29), while Fig. 30 shows the estimated marginal likelihoods and corresponding weights for the top-scoring 50 valid models. Finally, Fig. 31 compares the final estimated posterior p(2|ær, t) to "flat" average Pmat (2(
In this case, the top scoring model m (736) Fig. 27a models per-city mean temperatures and shared regression coefficients to predict test temperatures. Many generated models used more shared parameters. For example, m (688) (Fig. 28a) has global parameters for the distribution of temperatures and m (487) (Fig. 28b) also has shared autocorrelation parameter. However, these more complex models had lower marginal likelihoods. Thus the final posterior is essentially the same that of m (736)
28

---

data
citicit
int nun days, train;
nun.doys_rest:
array[nun. cities, nue. day..ttrosnnn real day2.. tew.-tenp.ttttt. array cities, nun. days. train] doy2-ttemp. train:
array(nun. cities, nue. days. test] real dayl. Temp.tests test:
paraneters
| real city_mean;
brrayfnun, cities, nun_days.train) teyp_temp.ad]: real beta;
real alpha;
rodel
beta nommal1..8 8.21:
alpha nomel(a, 51;
i-nuncitties)(
city_neon[s] (58,al 20);
(J T:mm days. train)
deyl. ter lij no 181;
day2 temp trainji, il nomalicity rean[i] alpha beta (dayl tewp.teerpttee.... train[i, 1] (ity..een(14] + dayl terp.adiie adi
generated quentities { array cities, tesst real tey2.terp.testt:
for (1 in 1:num, cities]t
heyp_tenp_testfa f1 norna_. rng.oltyynncitt..eenee.. alpha beta (dayl_terp.esttte f1 dity.neanfil), $1;
(a) m (736), the highest-scoring model.
data
nun. cities;
nun.days.train,
int nun.days... test.
array[nun cities, day.days. train] daya. temp
array crrayfnuu_citiin cities, nun_days, tmn_days_trainl real day2_tenp, trapin arrey (non, cities, nun. den.deys.test] dayl. tests
paraneters
(nrny real city.nean;
arrayinun. cities, man.days..trai real day.effect; realstower-.1, upper=ly change:
realclovernts tesp_sd;
rodel
for(1 1snun [
city.neen[1] normal(50,200
tarlj in daw..ddyy.tttann
effectif1) nornal (0,10);
terp. change uniform(-1,, 332
tesp_sd cauchy18, 51:
tarti 1snum. cities]l
tar(j in 1:nus days.trainj(
day2. trainti.]] tornaliiity reanfil day.effectfi..]1 terp.change tewpyltttt.triifi.. rety.re day dfy.effect(i..)) temp.sdi;
generated quentities { array frun, cities, nun days.test] test] real day2. tosp: forti leran citiest( for() in 1:num days.teysttett day2, temp.ttttettteeeee testis,11 noral, rngicity, [ea] crrnne (dayl, temp, test[i,1] cdty_nean[s)). tosp_sd);
(b) (785) aciond-highest scoring model
Figure 27: Example models for the city temperature problem.
29

---

data
non.cities;
int mun. wun.deys.traim,
days. test;
arraytnun, cities, nue mun_days_tral real dayl_teng, temp.trattnn
array[nun, (nunccctiionn cities, mun_days_
arrey cities, nun. days. test] real Jeyl_terp..ee
paraneters
real
realslover-d- sigeD:
array[nun_ cities] real city_mu:
arrey (nunc cities) resteloverm-l, uppersle city.ary
sodel
normal(59,
sigra nomalç 28]:
1snun_cities][
city_mufe] norsal sign));
city.ar[e] uniform(-1, 1);
for tyssnt
tepp.tenp.trainfcc d] nommlicill.roftc) city.ar[c] (enyLlenn..teinnfe, d) city.nufe)), 53s
generated quentities {
arrey cities, nun. days. real day2. terp.tertttt.....ee test: 1snum cities)
for (d in 1:nux,days dans_days..tst)/
day2. tez-ttepp.testic, romattylglcttt..ulee city.aric] (dsyL.terp.testie, dj city-rutc]), 5)2
(a) Y71 (688) lower-scoring model.
data
int nun.cilies;
traint
int munudeys den.ddys..ess
array anray[nun. cities, nun. dedddyyytrrrii dayd. temp. arraytnun. cities, nue. day.yyyyttte train] real toy2..tenp.tree train: array(nun,citionn nue, day__ttsssss tost] hoyl_tonp_tost
paraneters
vectorinun, cities] nean;
vectarclovers0, upper=1=[nues coclll(nnnnciiieee cities] ettyyeutocorrllation;
Lrean
realxtouer=0, upper=1x buttcorratttc
real-towersds sigsa;
rodel
global.nesm n6emal
autellautocoorelatiin - unifom(9, 11;
cinuncc
city neynne (groal 181: city_ rly.auuocorrelationn for (j 1:nux days. trysntrain)
unifora(0, 11:
day2. temp, trainji, 11 noymaldayl terp. trainii, 11 city. cttyssutocorrelation1 + city rean[i] + (1.
attpyautccorrllatinn signo);
sigea cauchyla, 51:
generated quantities {
arrayinun, eities, nun days, real toyp,torp.tte test: cininec
for (j in 1:nux dew.ddyyyttteeee
Jay2.tenp. www.ttttttttpppp. tastit, 31 noratt Dommt_.ragddop/ ragiday1. tomp.tost/1. 11 city (ty.auuttcorrrlatiinnt)) city.reants) (1
city autocorrelation signa);
(b) (487) typical but lower-scoring model.
Figure 28: Example models for the city temperature problem.
30

---

31

---

300-
100
count
10
3-
1-
600
400
200
log log (ogpmmmm)))
Figure 29: Estimated marginal likelihoods p(z (m) for each of the valid models for the city
lemperature problem. Bars are colored for four example models from (Fig. 27) and (Fig. 28).
150
1.00
M (
log p (m (n))
0075
50
-0.50
weghtt
-0.5
700
710
720
730
0.00
model number
Figure Estimated Nog-marginal likelihoods p(z/mm()) (n)) and weights w(m(1)) for the top 50 models for the city temperature problem. markers are colored for four example models from Fig. 27 and Fig. 28.
100
1-+ f-+
47+
50 OF
100 OF
50
100 oF
Cities
50 OF
100 oF
+-#
50 oF
100 oF
--+ 4-+
4-4
50 OF
Test days
p(z|x,0) t)
(nat t)
true
Figure 31: Medians and 90% credible intervals computed for each city each test day for the city temperature problem. In this case, both LLBB and Pfst have coverage of 44/50=88% of the true values, indicating reasonable calibration, though the posteriors under LLB slightly more narrow.
32

---

B.7 Gold (small)
In this problem, the user describes a "rod" of some length, with a varying density of gold atoms. Samples were taken at various positions and tested to be gold or not. Finally, a set of future test locations are given. The user "forgot" to give types for some variables.
PROBLEM
have metal rod. I've vartous positions along the rod and campledd single atons and tested if they were gold. or no.) I'd Like to predict gold would detected tried ressurenents at various test locations. essure the density of gold changes smoothly, but the actual censity could be quite complicated.
DATA
real Longth; Length rod
int train; // number tines rod was sarpled
array [nun. train] train. locs: where samples taken in training data train] train;
detected each location (0 if int nur test; // hypothatical test locations
array [nun_ test_locs;
future locations sarpling
GBAL
array [nun_ (est] test; // would gold be detected at future lecations
Synthetic data was generated by taking training locations uniformly from 2.5. The true density of gold was a (1 + sin (9x) x=²)) where (s) 1/(1 + e²*) denotes sigmoid function. The test locations were 100 points equally spaced from 0 2.5. A sample of the data is shown below, while the full data is plotted in Fig. 32.
'rod_length" 2.5,
"nun.train": 30,
*train, 28.8647744441726, 8.916637936887462, 188597733118777252,, trainn 10, 1, .... 11,
5.34499994338811
test":
"test. locs': 10:0, 0.9222222222222222
0.85993000000000000 2.5),
1.00
gold density
0.75
true density
0.25
gold detected
not detected
0.00
0.0
0.5
1.0
1.5
2.0
2.5
position
Figure 32: Observed data for the gold (small) problem.
In preliminary testing, it was difficult for the LLM generate valid models for this problem. Thus, instead of 1000, 24 16384 candidate models were generated. Of these, 192 (1.2%) compiled and allowed for inference. The top three models and low scoring model are shown in Fig. 33, Fig. 34, Fig. 35, and Fig. 36. Histograms of estimated marginal likelihoods are in Fig. 37, while Fig. 38 shows the estimated posterior for each model, and Fig. 39 shows the weights for the top 50 models and the final estimated posterior.
In this case, essentially all posterior weight is given to the top-three scoring models, m (192) Fig. 33, m (191) Fig. 34, and ml (190) Fig. 35. The final posterior appears bit better than flat average.
33

---

data {
rod, Length;
1nt nun_traio; errey(num
loci:
arrayInun, train] int-louerr0, upper=l> gold. train;
int run, test:
arrey(nun. test] real lost.
transforeed
int-looer=l> nun_train nus,test;
errey(N]
append treayy locs, Lescl)
paraneters
vector[N] realetower=s> rha;
trans formed paraneters vector[N]
INtrix(N,
for
KES, 1] 1 8.61; // diagonal etement
for 1:f1-133
real. dist sert Sort(aaareeall K[A, (Brt(B///orttBrrhool Bathho/opooi-aaatttdda////t/t
all tocsssijj
K[j, Kli, 21: // fitl in lower triangle
K etaz
nomel10, 11:
rha noraal8, Length)) // use rod length to specify the input scale the function
for in 1:nm. train) gold.train(i] Pemooutiitogitiiiii
generated quantities
arrayinun.tese tost] intslover-o, upper-1> gest.
for test) gold. testlil tegit. train ill:
Figure 33: (192), the highest-scoring model for the gold(small) problem.
34

---

data {
rod, Length;
1nt nun_traio;
errey(nun,
loci:
array[nun, traann) int-louer=0, upper=l> gold. train;
int run, test:
arrey(nun. test] real lost.
transforeed
// creating ing a grid of locations to corpute covariance satrix testt
array[nun. tatent] real alt.locs:
for(1 traint
[lL] (4)):
fort: Iinin_testje all_tocsfienun, train| locsfli
paraneters {
vectorinun. Latent]
transfored
vector(nun. latent] f;
reurix(nun latent, mun.latent] K; for (1 in latent)
latent)
real dist (art( aot_tostf])>21:
sert(3 distil exp(-sqrt(3
}
180 latent) deoleskyydeconppseelle
sodel [
-stddnomell);;
for(i in 1:nux trainj(
gold troins] bernotlis(aw logitif(1)));
gonerated
arreyfnun. testl intelovernO, upper=ls gold.teste test; for(L testif
real true. qold, density inv [f train]]: rngitrue. gold_test[s] tost[s] bernouls,, tertuuts, wnggtrue.godddeeee gald_denstty):
Figure 34: (191), the second-hiiheest-scoring model for the gold() problem.
35

---

data {
rod, Length:
1nt trantrain;
arrey
loci:
arrayInun, tranttn int-louer=0, upper=l> gold. train;
int run, test:
arrey(nun. test] real lost.
transfored
int nun, nun_Tocs nue_ train nunttsst; arrey(num. locs] alLss
Tor(1 1:nux train] atL.tocsf]] (i))n
tor(1 1:num, testil all_tocs(Sinun_tee test_lacs[1];
parareters resteloversde length;
vectorinun. locs]
transforeet
vec vector(naa cetttorrnun..ooo locs]
retrix (nun. locs, nun.locs] for (1 in locst
locs)(
KES,J1 - (ortt(( (arrsts)-ll_eeecttte /onghh))/
K
rodel
(9,nat(9,11)
length normel(0.1.0.1))
for (1 in 1inux. train](
gold troin]] bemoottss2/(Heeetttee
gonerated
arrey[nun. test] inteloverno, upper=ls test:
test;
vector[num. locs] smooth;;
pisnnoth a inv. logitip);
for
test][
gold_ test/]] be bernoulli__ temm
Figure 35: m (190), the third-tighest-scoring model for the gold(small) problem.
36

---

data {
rod, Length;
1nt nun_train;
arrey (nun,
loci:
anrayInun. trantnn int-lluer=0, upper=l> gold. train;
int run, test:
arrey(nun. test]] real lost.
real-d real-towersde rho: // length scale for the kernel restslover=te alphe; scale kemne vector[nun. train] eta; // Latent variables
transforned paraneters vectorinun, troin] nu:
votrix(nun train nunu train] K;
train)
KЕЕ, 1 alpha: for train)
alpha exo(-((train.locak) locs(s) tortnn)ocs))))"))) moo2))3
K els:
sodel {
couchyl8, alpho cauchyl8, 51: normel(0,
gold.trein - bernoulli. toritttl
generated quentities pectorinanttes]] test.nu;
vatrixfnum test, nun. nun.tran train] traast lests
for in 1inun. trein)
traann, 11 alpha expt-((test,locsss)) tocs:]] / mor211
test.nu (K.test. tratn
arrey(nun.ti test] intsloversO, uppersis test; far 1inux. test) gold_tost tost[s] tog.ttts_ooott.nngggooooooo nu[s]);
Figure 36: m (171), low-scoring model for the gold(small) problem.
37

---

300-
300-
count
count
10
10
3-
3-
1-
1-
200
100
100
80
85
100
105
log ((x
log p(x|min)
Figure 37: Estimated marginal likelihoods p(=|mm (m) for each of the valid models for the gold (small) problem. Because of the many order of magnitude, two different ranges are shown. Bars are colored for four example models from Fig. 33, Fig. 34, Fig. 35, and Fig. 36. multiple models map to the same bin, all colors are shown stacked.
1.00
gold density
0.75
0.50
0.25
0.00
0.0
0.5
1.0
1.5
20
25
position
Figure 38: The estimated posteriors [[) (m)
x for each the valid models n for the gold (small) problem. Lines colored for four example models from Fig. 33, Fig. 34, Fig. 35, and Fig. 36.
6* 0.6
1.00
100
weeign ht w(m inly
log pX m())
gold density
0.4
50
0.50
25
0.5
E[Z|x, t] under pizíx, t)
Fizix, t] under Pralz((x t)
0.0
0.00
150
160
170
180
190
0.0
0.5
2.5
model number
position
Figure 39: Left: The log-marginal likelihoods and weights for the top 50 models for the gold (small) problem. Right: The final estimated posteriors, compared to a "flat" average. The true density is shown for reference.
38

---

B.8
Gold (large)
This problem concerns the same gold problem, except with larger dataset with more observations. Synthetic data was generated by taking 150 training locations 2 uniformly from 0 to 2.5. The true density of gold was (1 sin( f xx²) where ((s) 11((1 + e e*) denotes sigmoid function. The test locations were 100 points equally spaced from 0 to 2.5. A sample of the data is shown below, while the full data plotted in Fig. 40.
*rod. Length": 2.3,
train":
"train_locs" (B.666344472396544, 0.449875883313362,
1.8815478996618811
trainnt ....
"nun.tast": 186,
Tocs": 1e.0, 0.85525295955525555. .... 2.51,
1.00
density
0.75
gold golddensit
0.50
density gold detected
not detected
0.00
I1
l H Nn 1 0
I I
0,0
0.5
2.0
25
position
Figure 40: Observed data for the gold (large) problem.
The same 214 = 16384 candidate models were considered as in the previous problem. Of these, only 71 (0.4%) compiled and allowed for inference with this larger dataset. The top three models and a low scoring model are shown in Fig. 41, Fig. 42, Fig. 43, and Fig. 44. Histograms of estimated marginal likelihoods are Fig. 45, while Fig. 46 shows the estimated posterior for each model, and Fig. 47 shows the weights for the top 50 models and the final estimated posterior.
In this case, essentially all posterior weight given to the hop-three scoring models, m (71) Fig. 41. This produces a final posterior quite close to the true density, and much better than a flat average.
39

---

data {
rod, Length;
1nt nun_traio; errey(num_
loci:
arrayInun, traint int-louerr0, upper=l> gold. train;
int run, test:
arrey(nun. test] real lost.
transforeed nun_train nus,test;
int-looer=l>
errey(N]
append treayy locs, Loscl)
paraneters
vector[N] realetower=s> rha;
trans formed paraneters vector[N]
INtrix(N,
for
KES, 1] 1 8.61; // diagonal etement
for 1:f1-133
real. dist sert Sort(aaareeall K[A, (Brt(B///orttBrrhoo Bathho/opooi-aartttdda////ttt
all tocsssijj
K[j, Kli, 21: // fitl in lower triangle
K etaz
nomel10, 11:
rha noraal8, Length)) // use rod length to specify the input scale the function
for in 1:nm. train) gold.train(i] Pemnoutiitogitiiiin
generated quantities
arrayinun.tess tost] intslover-o, upper-1> geld.
for lsnux gold. testlil togit. train ill:
Figure 41: m (71), the highest-scoring model for the gold(large) problem.
40

---

data {
rod, Length;
1nt nun_traio;
errey(nun,
loci:
array[nun, traaan) int-louer=0, upper=l> gold. train;
int run, test:
arrey(nun. test] real lost.
transforeed
int-looer=l> N = nun_train nus,test;
arrey(N]
append trray locs, Locs))
paraneters
realsloversdo rho: Length scale Metern covariance realetower=s= alpha; // scale for Natern covariance vector[N] latent value
location
rho cauchy10, 21: prior for length scale
alphe cauchy18, 21: // prior for scale multinnormal (roporettroeeeetoorreeee
gp.nstern32, covialt, locs, sipho, rhol): // sample fran GP prior for 1snun_train](
told_trein(j] bemoulli
logitleta(i)); -emmottiitoiiitteeeiiill
generated quantities
tast] intelover-o, upper=l> gest. for 1:nun_test){ gold_test(i] termotttittgiii.nnnieeeaa nantrainl1:
Figure 42: m (70), the second-higheste sigiigt model for the gold(large) problem.
41

---

data {
rod, Length;
1nt tunntrain;
arrey train] loci:
arrayInun.tren int-louer=0, upper=l> gold. train;
int run, test:
arrey(nun. test]] real lost.
pectorfnun_trainn latent;
realslover=de sigra; real-douersss rho:
rodel {
sigro couchyl8, 21:
rha cauchyl8,
{
ratrix(nun train, train)
for (1 train)
K[a, = sigre"2;
for
KE3, 11 = sigrar2 expt-([tratn tocsta] train, ttratt_lesssppp. locs(51)"2) / rhor2))s
31:
latent multi (rep.
vector (rep.
train),
gold.troin Logpntttis.lootttaaaat
generated quentities
vector[nun. test] test. latent;
retrix(mun.ti test, nun.train] Kitest: retrax(nun.train, nunutrain] K.train; linum.train K.trainfi, 11 1 signar;;
K. K.trainj,, 1] siona"2 expt-((train locs[A] train locs(1)) / (2 rhor2)):
K.troinij, Kutrainjs, 11:
1snux
for (J. in Isnun. train)
Kitestfi, signe"2 locs))) (ocs) rho*2));
test_1atent (K.test InverselK.trainn) latent;
arrayfnun. test] intcloverso, upper=ls gold.test: for (1 in lesur.test) {
test[i] bernuuli, Logit. (tgttest latent[i]); (
Figure 43: m (69), the third- hhrrd-highest-scoring model for the gold(large) problem.
42

---

data {
rod, Length;
1nt trantrain;
errey
toai:
anrayInun, [rurtann int-cluer=0, upper=l> gold. train;
int run, test:
arrey(num. tesit] real lost.
real-towerrde gigsigr:;
reetslover=do gp.rhe;
vector[nun. train] gp-2:
transformed paraneters vectorinun,trattt eta;
gPu7:
for (1 in traun_train)
zun
for
Tsnun Irnux.train)
sun == gp.sigre exp(-((train tocin.locs(i]) op.rho) gp.2(i):
eto[i] sum:
gp.sigre inv www.gereefee (1;
gpurto (sw.gasraaa1111
nomaa[8,2]] gold_train = bernaulli_ (onnoltt_oogtt(eta))
generated quenticies { vectorinun test] test; for lesus
real sun 8: linun.tratn)
sun am 3p sigra exp(-[(test tocsti] lorsff] / rhoi gp->f31:
cta. test[i]
arrayInun_ [ect]] int-looer=o, upper=1s tost; for (1 in t:nux_test) gold tost]]s togoettt..ogttt.nniete. test(1]);
Figure 44; m (31), a low-scoring model for the gold(large) problem.
43

---

1.00
gold density
0.75
0.50
0.25
0.00
0.0
0.5
1.
1.5
20
25
position
for each of the valid models for the gold (small) problem. Lines are colored for four example models from Fig. 41, Fig, 42, Fig. 43, and Fig. 44.
Figure 46: The estimated posteriors [[a|m (m),
1.00
log pix|m mn))
100
1.00
Vagg ni win (u)
density
0.75
0.75
50
0.50
1000 (d
0-
025
025
Elz|x, t] under piz[k, t)
Elz|x, t] under Praileix,))
-50
60
70
0.00
0.00
0.0
0.5
1.0
1.5
2.0
2.5
model number ber
position
Figure 47: Left: The log- marginal likelihoods and weights for the top 50 models for the gold (small) problem. Right: The final estimated posteriors, compared to a "flat" average. The true density is shown for reference.
300-
300-
100
100
count
count
10
10
3-
3-
1-
1-
-400
-300
-200
-100
100
40
60
100 log (x(mm()) in))
120
140
log (n)
Figure 45: Estimated marginal likelihoods p(z|m(() (n)) for each of the valid models for the gold (large) problem. Because of the many order of magnitude, two different ranges are shown. Bars are colored for four example models from Fig. 41, Fig. 42, Fig. 43, and Fig. 44. If multiple models map to the same bin, all colors are shown stacked.
44

---

C
Additional analysis
C.1 Form of the self-normalized importance sampling weights
Consider self-normalized importance sampling with target distribution of p(m]ar; and proposal distribution of p(m](). The naive importance weights would thus be
p(m (e..) p(mmmije)
(n(n)
Σ plm(n
p(m((a))
However, we know from Eq. (5) that p(m|ar pr 1) x p(m()()(())m)) or, more precisely
p(m toolttp((mm)) (m)
p(m]ttp(alm)
t) =
Em p(mm()))(((mm)
p(zlt)
From this, obtain that
p(m|2z,) t) p(a|m)
p(m(() p(al()
Substituting this into the above form for w(n), obtain that
w (n)
p(2(()
N plalm(n))
=1=11
p(a(()
p(x)(m(n))
C.2 Variance of the seff-normalized importance sampling estimator Consider the importance sampling estimator from Eq. (11), i.e.
i w g(m(n))
Standard theory [e.g. 25, Eq. 9.8] says that the variance of ja is asymptotically
Ep(mle)
((m)² (g(m) - -n)
lim NVIA]
N00
Ep(mlt) [m](m]]
where w(m)= p(m] This simplifies into
p(m t)
p( p(m (g(m) - 14)
Ep(mle)
lim NV N000
p(mlz.e)]
Ep(ml()
p(m t)2
E
(a(m) 11)²
m(mle)
p(m[1)2
Now, if assume that [g(m) pl] < 8 then this be bounded by
pp(m |z, t)2
lim NV N
s 8² E
p(m
p(m]t)²
8² (1 x² (p(mlz,.))ll(m(tt)
45

---

C.3
Form of the self-normalized importance sampling weights with variational inference
Consider trying estimate
((222) = E g(2 ar m),
((m)2)
where q= m) is some distribution and ((mlar)) is as defined in Thm. 2. We are interested in estimating this expectation using SNIS with the proposal distribution p(m](t. It follows that we should use estimator
@(=|<>) Σ wh ( g(zlar (
where m( ,.., ,m())
, ~ p(ml(), and the self-normalized weights are
((m(()) (x)
(m)
p(mim) (t)
g(m(n [x] n==1 p(m(n")11)
Now, observe from Eq. (19) that
p(mlt) exp
ELBO m)) |p(z, m m2 t))
((m/a)= = Em p(m' i) exp
ELBO(g(Z)z, m') llpzz m m1 t))
where the denominator does not depend It follows that the weights
exp (LL
(x, t))
(m)
E exp
EIBO m(() m(()) [pp(( mn())) Ix. (t,tt))
C.4 Analysis with inexact ELBO estimates
Theorem 3. Suppose that p(2, mlt) and q(2/22, m) are fixed and choose
q(m/x) x p(m)z) exp
8(m)
where 8(m) 2 O represents the "slack" in between the true ELBO for model m and the bound used for computing (ml)x) Then the resulting joint divergence
KL(g(z, ml x [lp(z, ml m|x,t))
log
E
cxp
KL (m)
((m) +
(21)
p(m(ar,))
log p(z]t) - log
E pml())
exp
(ELBO((M) 5((m) +
(22)
where 5 = Eg 8(m), and (m) and ELBO(m) (m) represent the KL-divergence and ELBO from Eq. (17).
(Proof in Appendix E.3.)
To interpret this result, first note that if all ELBOs are exact, meaning g(m) O, then Eq. (21) is equivalent to Eq. (16) while Eq. (22) is equivalent to Eq. (20). This also shows that using "improved" distributions and tolerating "slack" in the ELBO can only help, compared to using distributions where the exact ELBO has the same value. That is, no matter the value of the slack variables 8(m)) Eq. (21) is never worse than Eq. (16) and Eq. (22) is never worse than Eq. (20).
The obvious implementation of Thm. 3 would be be algorithm that loops over all models. This is given explicitly as Alg. 7 in Appendix D. Again, such an algorithm is intractable, but can be made practical by combining it with SNIS. This is given Alg. 8.
interesting special case when MCMC is used to approximate each posterior p(2 m) and we assume that these samples in fact come from p(z|0z m). (See Alg. 3 in Appendix D.) If this case, the
46

---

best strategy is simply to use the best possible lower- bound on the marginal likelihood, computed by any method. It's see in this case that the joint divergence reduces to
log Σ p(m]x,t) exp (-8(m) + 8).
From this, we can see that if all errors 8(m) are equal, then 8(m) = and the joint KL-divergence remains zero. Also, if there are enrors on models where p(m æ. is negligible, these also have negligible impact on the joint divergence. However, if there are different levels of error models where [m]2 1) is large, this could impact the joint divergence.
C.5 Relaxed bound
Corollary 4. Under the assumptions of Thm. 1,
m)))
()) S
E p(mla,,)
Proof. We have that
log exp (-
E (1(g(2),,, (ill(al)m p(m|a.1)
E p(m|a.,)
kl(g(=/>, m)) m) (m))
log exp m))a(l p(m|2,))
kl(g(2, mla:) m| t))
The first line is obvious. The second line follows from using Jensen's 's inequality to see that
E (_p(-K(((2)z, m) m) m))) 2 exp(-
[l(g(z|,, m))lp(z(ar,))
p(m|z.()
p(m[æ,6)
The third line is the form of kL(g(z, (pl(,, m/z,, ()) from Thm. 1.
47

---

Algorithm 4 LLB with exact inference and SNIS.
1. Input textual description and data z.
For -,,,
(a) Sample mf p(m[t).
I/ using LLM I/ under PPL.
(b) Compute marginal evidence (()(m(()) and posterior p(z|22, m())).
Set w(()) x p(x|m!!)) where Σn w(() 1.
Use final posterior approximation
N Σ m=1
p(2 ~
w (n) p(z(z,
Algorithm 5 LLB with variational inference and SNIS.
Input textual description and data
2. For 1,,,
(a) Sample mn(") p(m[t).
II uuiigg LLM
(b) Maximize ELBOO( over some variational family ((2/0, m())). under PPL
3. Set (())) exp (mL where En w!n) 1.
4. Use final posterior approximation
p(2 (e,,,t) 1)
(n)) ((=/o,m!'))
D Other example algorithms
Algorithm 6 Variational LLB with exact enumeration of models (theoretical)
1. Input textual description and data x.
For all possible model strings m:
(a) Compute model probability p(m[t) (using LLM)
(b) Maximize ELBO(g(Z)z, "/0) over some variational family
3. Set (mm) p(m]t) exp
4. Use final posterior approximation
p(2lz,t) R Σ wo (m) (m) g(2|o, (())
48

---

Algorithm 7 LLB with variational inference with inexact ELBOs (theoretical)
1. Input textual description and data z.
For all possible model strings m:
(a) Compute model probability p(m[t) (using LLM)
(b) Compute some approximation to the posterior g(zlar, m) so p()
(c) Compute some bound on the ELBO
p(z, x/mm
(m)
E
log
((( m) logp(r)m).
((z/2,m)
Set w(m) x p(m[t) exp (L())).
4. Use final posterior approximation
t22 se Σ g(
Algorithm 8 LLB with ELBO bounds and SNIS.
1. Input textual description / and data
2. For n = 1, 1,2,
(a) Sample (n) ~ p(m[(t (using an LLM),
(b) Compute approximation the posterior ((2/,r, mlm)) p(
(c) Compute some bound on the ELBO
p(z, p(z,z/mm') xlm
log
(21|,, m(nl) I logp(z/m("))
3. Set (n)n) xx exp ((n)
Use final posterior approximation
N p(2|0, t) ~ Σ
wn m())
==1
49

---

E
Proofs
E.1 Proof of Thm. 1
Proof. First, notice that minimizing
m
m/,,
is equivalent to maximizing
p(z, x, m, (t) q(z ar, m)
log
We will therefore attempt to minimize the latter. Set up Lagrangian to find optimal ((mla)).
L Σ {(m])))(())m x log p(z,,, x, mll)
de A [m(ll)) - ((mlt)g(=lmmaa)
DC dq(m)
f (()lm)log 2l2,2,m/t dz / q(ml)))(2))m))) dg(ml(t
d
log q(m)t)dz X
(( x)
/ ((z)m) log 2(2, m1
de 1
((m]f) (2)
/ ((2lm) log
p(mlt)p(zl/mpp((p... m) de 1 -
q(mit)q(2|mmmz)
p(m[t)p(x/mm) kulg(lla, - 11 A
log
(m]])
Solving for Zg(m)
dC we obtain that
(m)(() x p(m|0p(x/mm) exp ( m) (a m))).
Now, for simplicity, define kL(g(2|2, m)]lp(2(2, m)) Am- Then, write
((m)a:) = p(m (m) exp (-4m)
Emo p(mn p(((()m exp -Amn))
p(mlx,.)
Ems p(mma())p(mmmm
Now, note that
= p(mlt)p(z((m) exp( Amm)
Emo p(m'1t))p(((m)
p(m|x,t)
Em p(m' (onn)(()((r//mm)) exp -mmr
p(m]t)p(z/m)
exp (-Am) Em p(m (mmmmm))))) (mmm(m)(((pmm))))
Emo
exp - -Am²) exp (Amm)
En p(m' exp (Amm)
Thus have that
50

---

E log g(m(k)) E log g(=lma
kilg(2, malllsmm)olloo
g(mlar))
p(m(x,
((z,m/)))
ple m,,))
exp
Am)
log
Am
a(maar)
Σ m² p(mm lor, () exp
-mmm)
a(m a))
log
Em p(m') expp -Amm)
= llog
E exp(-Am)
p(m(x.1)
Note that since Am 2 0, it follows that exp(-Am) < and thus Ep(m|z,) exp(-Am) S 1 and so log pp(lla.) exp( - Am) ) S 0. Thus kl(g(z, m|x))p((2,, m|x: 1)) is non-negative, as required.
E.2
Proof of Thm. 2
Proof. Recall that the Thm. 1 found that the optimal g(ml2z))
q(mlz) x p(m[()p(([mm) x exp (- m))l(tz))
(23)
And recall the "ELBO decomposition". i.e. the fact that
logp(z|m) ELBO(g(Z)2, m))ppz, xlm)) + kL(g(z/2, mm)l:// m)).
(24)
From which it follows that
p(x/m) exp mlll(ii,,, m))) (2)
Substituting into Eq. (23) gives that
(mlz) p(mlt) exp ELBO(g(ZIz, (p(z,, t))
Recall also that Thm. found the joint divergence resulting from using the optimal g(m/x) is
ml))llzz m|x, ())
log E exp --KK(g(tlosmmmlllzzizz m))
p(mlz.1)
log E exp (LLB)O((22
p(gp((m))
p(mlz.()
p(m|ær,)
log
exp exp expp|ELLB00(ggllssm)))tpttmm))))
p(z|m)
p(m[t) p(xlt)
log Σ
exp (
log p(x(()
log E exp (ELBO(((2),, p(m)e)
Where we have used that
p(m|x,1) p(m p(mlt)p((mm) p(x(()
and so
p(m]art) 1) p(m)t) p(x/m) p(x]t))
51

---

To understand this, suppose each was exact so that
(LBO(g(Z)x; m)) [l((, x//m) = log ((g
Then
kLlg(2, (ilol,,2m logp(<;(() log IE p(x/m) p(mit)
logp(x](t) log p(rlt)
0.
E.3
Proof with inexact ELBOS (Thm. 3)
Proof.
Theorem 5, Suppose that p(z, 2, ml() is fixed and q(z]a, = p(z|0, m)) is exact but ((m)z) oc
p(m]t)p(z(t) cxp(-8(m)) is not exact (e.g. because it is chosen based on ELBOs rather than true marginal likelihoods). Then,
<14g(=/=)llo(((oo00) s (I(((m))) [p(m|ar,)) 1))
log E exp( exp(-8(m)) p(m[z,t)
E
8(m).
Proof.
log g(m|z))
kulg(m/=))l(mmmz..)) (m)
a(m mn)
p(ml(,,(
p(m (p(((1(()) exp(- log En p(mllt)((x))))) (m) exp(-8(m))
En, p( ((mm)(((m)
E log
(m ar)
p(mlt)p(z((m)
Em p(mlt))((//m)
8((m) log En p(mlt)((r(m) expp
E
q(m ar)
Em ppm()p(mlm)
log
Em p(m)t)p(())) m)) exp(o - 8(m))
p(rlt)
log
Σn p(mm((()))p)) 1p exp( - 5(m)
log
E exp( ((m)
p(ma(e,)
log E exp( (8(m))) 8
p(m|x,()
Jensen's inequality (- log is convex) gives us the easy upper bound ku(g(mlz))lp(mmaa t))
llog E
cxp(ō 5((m)
p(m|x.))
<
E log exp ( 8((m)
p(m|a.1)
E g(mm)
E 8(m))
p(m]a.))
g(m x
To understand this, note that the difference ofp(m/a, t) and g(m|2) (z) is precisely that g(m (z) is smaller when (11) is larger.
52

---

F Experimental details
Models were generated using Llama-3.3-70B [14, 21] with 4-bit AWQ quantization [19]. In testing,
LLMs seemed much better at writing Stan code [5] than other PPLs like NumPyro [26] or PyMC [1], possibly due to more code being available. We thus used Stan, though the possibility of creating unnormalized distributions in Stan poses some difficulties (Appendix F.2).
The LLM system prompt (Appendix F.1) explains the PROBLEM/ DATA / GOAL format illustrated in the experiments and asks the LLM to first write a THOUGHTS block explaining how plans to model the problem, followed by MODEL block of Stan code. For in-context learning, provided six example inputs along with high quality outputs (Appendix G). We rejected any models that did not compile that used certain language constructs that might lead unnormalized models (Appendix F.2).
Models were generated using single A100. With continuous batching (parallel inference), model generation was reasonably fast-for example, it took around 14 minutes to generate 1000 models for the city temperature problem and around 11 minutes for the polling problem.
For inference, it necessary to reliably and automatically approximate posteriors and marginal likeli- hoods for thousands of models. To approximate the posterior, we simply used Stan's 's default NUTS [17] sampler with 10,000 iterations. To approximate the marginal likelihood, we used importance- weighted variational bounds [3], with Gaussian proposal distribution. In general, it is known that such bounds are improved when the proposal distribution would minimize the x2 divergence to the target. However, given the difficully of minimizing such divergences Geffner & Domke [13] elected instead to use a proposal distribution q , m) minimizing (L(p(z|ar, m) [] az, m)). which, given samples from p. amounts moment matching, i.e. matching the empirical mean and variance of the samples from MCMC. The full inference strategy is given explicitly as Alg. 9.
Algorithm 9 The hybrid MCMC/ / /I/ SNIS algorithm used in the experiments.
1. Input textual description and data a.
For 1,2,....44. .N:
(a) Sample mn() p(m(t) using LLM.
(b) Draw samples 2(n,1), (m, ~ p(z|æ, m()) using MCMC.
(c) Set q(2|.,, mn)) to Gaussian with mean and covariance matching the sample 2(.11) .... (n.K).
(d) Estimate the wmportance-weighhed ELBO (()) = É log pizla.m
g0 where the inner ex-
pectation is estimated using 25 samples q and the outer approximate expectation uses 10,000 repetitions.
3. Set (n.k) x exp ("),
Return the set samples {z(n.k) where (m.))) given weight (n).
Measuring the amount of compute used for inference is difficult since it was run on heterogeneous cluster. However, running MCMC inference on 1000 models course quite expensive and required tens of hours of CPU time for each model. Simply compiling 1000 Stan models (on CPU) is more time consuming than generating 1000 models using an LLM, though of course both of these steps are embarrassingly parallel.
53

---

F11 System prompt
The following system prompt was used when generating all models:
You are StanWriter. Given a description of problem, you write code.
You always first write THOUGHTS and then describe your model in words. Then write MODEL and write Stan code.
Be creative when coming up with your model!
When you declare arrays, ALWAYS
syntax "array[2,5] int NEVER syntax *int x[2,5]
In the generated quantities block ALWAYS the rng syntax like "x normal (mu, sigma)" NEVER
"x normal(mu, (miggm))
NEVER write 'stan before your code.
ALWAYS give priors for all variables. NEVER use implicit/improper priors.
NEVER "target syntax.
F.2
Validating models
We rejected any model that did not contain the string MODEL, where the text after model failed to compile.
Using Stan for the purpose of this algorithm, has one significant disadvantage: is easy in Stan to create models that are normalized. This of course poses issue when doing MCMC sampling, but is significant here since it can affect the marginal likelihood. Unnormalized models can be created several ways:
1. By default, Stan uses flat / unnormalized / improper priors.
2. One can use the same variable in sampling statement more than once, e.g. one may write "x normal 10) and then later write "x normal Both of these statement increment the density, meaning the result is no longer normalized.
One can directly manipulate the target using the "target += construct, essentially leading to arbitrary changes with probabilistic interpretation.
4. One can transform variables the "left-hand side" e.g. write "log(x) normal (0,1)," This again leads to an unnormalized density.
If variables are constrained in Stan, this constraint is not reflected in the normalization.
Except for the #5, none of these are too critical, and one can use system prompt (and in-context examples) that discourages use. This works quite well, but (probably due to frequent use example models) #3 remains common. We thus reject all models that contain the string "target ", as if they did not compile. Checking usage of #1, #2, and is not trivial to do programmatically. and thus we relied instead spot checks in models given high weight. Only a handful of cases were found.
This leaves issue #5, which unfortunately difficult avoid due to frequent usage of constraints in practice, c.g. to sample from truncated normal distributions. We simply accept this bias in practice, though note the bias is always downwards, e.g. towards lower marginal likelihood.
54

---

G Examples
This section provides the example inputs and outputs given the LLM aid in-context learning.
G.1
Weight
INPUT
PROBLER
Given the weight of bunch of students, infer true meem of all students. Use reasonably informative priors.
DATA
1nt students;
students
errey(nun, students] real. StudenL.weight // weight of each student in kitograms
GOAL
real weight_sean;; true mean weight for all students
OUTPUT
I will model this problen by creating latent variebles. "weight.rean reflect the true neen weight of all students, and "weight. signa'
reflect standard deviationn assune norsel prior for "weight. rean" with sean of 80kg and standard deviation of 20kg. For the 'weight.signa" I'll guess standard
instructed reasonably informative I'll
Lognornal prior mu=log(10) signeslog(20). around 19 ka, but could be significantly larger or snaller.
MODEL
int nue.s tudentss
array (trassnun.studeets] real seudent.weight;;
pa ronoters
real neight_neen:
realslower=t veight. sigrs:
nodel
weight. normali90,,20
weight signa (1ggoorol((oggg1))
(tudentss
wiight nornel(weightL..s weightsignal)
55

---

G.2 Future weight
INPUT
Güven the weight bunch of students, infer weight of new student drawn the same distribution. reasonably unforative priors.
students; number
array[nun. sTudents] real. student weight; // weight of each student kiTograns
Student, weight; weight student, sampled the distribution
OUTPUT
dato {
students;
array[nur_ [tur_stuentttt real Srudent_wwighh::
pa naranaters
real weight. neen:
realslower=s= veight.
rodel
weight_rean- normal normal|70..011 370;301;
weight signo 12255,0,5) // heavy
for (1 in 3onun students]{
weigtnt. nornel(weightL.een signs):
quantitles
resl. new. www.sttuden..eigghee normel. rng/eeihhgh veight.
56

---

G.3
Movies
INPUT
Güven bunch binary features novie, predict user will not.
DATA
1nt nun train; // number of movies to train on int test: number movies int nun features; number of binary features for each novie arrey(nun, train, [untfuat int train. festures: // binery Testures array[nun. train] train. Like; // arrayinun, train, resturrs int test, features; // binory teatures
on
like the
GOAL
anray Test] int test. like; // does the user like each novie the test set?
OUTPUT
I will model this problen by creating vector "beta" of size "nus. festures'. i Mhen the inner-product of this vector is taken like that novie. I'll sample eech of that
with the features vector trom a standard normal, which seens Like reasonable scale when being combined with binary teatures. I'll then sample 'train, 1ike" "tost_like" by taking inner-product "beta" respectively.
it logit for how Likely the
"train_ teatures' "tast, features
MODEL
int nur, train;
int nue. features
array [nur. train, features] Intxlower-8, upper-ls train. features;
array [nun.. train] IntelovernO, uppersl> train. array test, num. festures] intelloer=8, upper=1s fest.resture;;
pa ronotars
Tnun, features]
(0,
for (n in {:muntrrain)
bernoull1. hog.:tit..rr...ettttttttat..aee testurestn)
generated quantities
anrayInun, test) intelower--, upper-l> Like; for temun.testt( tiken] bernoull.. tngtt. wtw.://tt..mtttoooooooooooo wtw.:rrwww.eooorttttteetetee.eee beta];
57

---

G.4
Collaborative filtering
INPUT
PROBI
Glvenn bunch ratting of vartous
novies, predict users future movies. by
ferrring underlying features of and novies.
int users:
int nun.novies; // number of movies
irt nun ratings; number
ratings
arrey (aan(nnnn.atnnnng int // what user did rating
arrey(nun. ratinngs int novie: // what novie was rated arrayLnun, restures)
no)
Ant arrayl num. users, novies] all, ratings; //
OUTPUT
THOUGHTS
I will model this problen by assuming that both users and novies can be described in terms of "nun, features"-18 featuras. I'll create latent varsable teaturee that describes the reatures for each novie. will assune standard cauchy priors for all these features, indicating that most features are seall, but some wight be quite large, Finally, assure the probability given bernouti. logit distribution with create the desired output ratings',, sill create generating quantities block in which loop over all pairs of users
describes teatures for user
variable nevie.teatures
given user liking given sovie consisting the inmer-pruduct of the and rovie features, To
data
int
number of
int nur. rovies; // number of novies
ratings;; arrayTnux, ratinggt int user: // what user did rating
observed
array [nun. ratings] novie; //
array four. ratings] int like; // did user like the novie 0 if nol
formed datal
int nur, features 18:
remeters {
array (nun. users, features) user. feetures; array (nue, novies, features] real novie. features;
model
forlk in I
forti users))
user. festures(i,k) (0,chy(0,1);;
for(j in rovies) (etterfeat cauchy(0,1);
ratings)
array[nun. features] real = user festures(uzer(]]) arrayLnun teatures] real reaturess real pro..rrduct(x.)): Likein] bernoulli. toottssor))
generated quantities
anrayInun. users, nus movies] int all. ratings;
forl1 1inun usersit
for(i rovies)
arrayinun. Inun. testures] real user. fosturesstsss arrayinue festures] real
real score pro. alL.ratings(L..) bernoulli logit. rng(score);
58

---

G.5
Vaccines
INPUT
PROBLEM
There
There different types response rate for each trial type. Assune that different vaccines of the sane Type tand to have simitar response rates and different types vaccines usually response
different
1nt nun. types; // number of different types of vaccines
int trials;
arrayfnun trialsl int type; // what type of vaccine tested in each triat
errey(nun. triels] subjects; array[nun. trials] int responders; // number of subjects in each trial that responded
trial
array[nun.co types] response..
response
vaccine type
OUTPUT
THOUGHTS
I'll create hierarchicel model. The idea is thet each trial has distribution, Then the peraneters for the per-trial type distributions are all drawn from sone global. distribution. By sarpling the three levels information can be optinally shared.
that drawn fror per-triel-Lype
(individual,
per-trialitypa
global rates)
In nore detail, I'll färst create paraneters 'a' and "b" for each trial type. These determine the distribution over true response rates for different trials of each type. true response rate for trial 'I' "per. response ratefil', is drawn from beta distribution with paraneters "s[type[k]]" and "b[type[i]]', Then the lubserved) number of responders in each trial will be drown fron binerial cistriburion with 'nun sun_scojccss(() repeats and a "per.trial responssetss type
success norral
will distribution fon paraneterized by "au" and "L', This will be napped to "a" and "b" in the transforedd parareters block through an nultivariste
2-dinonsional
inverse logit trans formation. I'll sample 'au' that results fron sanpling froz broad inverse Wishart distribution and then computing the Cholesky factorization.
distribution
from distribution
MODEL
int nue, types;
triels;
anrayinun. trials] int type:
arrayInur, triols] subjects;
arraylnun.trialee respanders;
cholesky.fector.covv fector.cov[2] L;
vector[2] nu;
anraylnun, types] vector[2] 2: // paraneters for response rate dist for each type array (ocnynnnnrttttt triats] real<lower=0, uppor-l> per_trial_ responss___ss_______________eese rate;
trans formed parareterst
map
peraneters
dist type Interval
array Inun. types] real 8:
Inun. types]
anuayInur, types] real responso_ rate;
fortn [ypes)
www.looogtiitnn13222
binl inv. www.ttgiii>tn,211.
response.rateinl =fm3/(cd(]+bbml))
distribution a/feeb)
nodel.
nu not noLtinnora (15.00..11..00..0.555
wtshart cholasky(3.8, 18,181,18,111::
forin in Iseun. types) theel
forl1 trials]
pnr_trial_rrsponse rateft] betataltypelil1. bitypeiills:
responders[1] binonial (ubjesabje))s teratrial. rer.trral..espomeereeeeee resell:
3
59

---

G.6 Water
INPUT
I start kg water, Every day, some fraction of water (between 0 and evaporates, and meesure much is rensining. reesurenents have anything between and 1, but suspect it's something around 0.25, The tirst day with 1kg isn't measured.
of something around Predict how much water
day. It could be
DATA
Int nun days; array(nun, (nur, (eay days] reesurenents: // reosurements in ky
GOAL
real evaporation; // what traction evaporates each day?
OUTPUT
THOUGHTS
model problon by creating each day. I will sample this from beta distribution with alpha=1 and betae3 meaning mean of around 8.25. Gi the evaporation, I will use Toop to define an arrey corrent.vstee which the true (unobserved) arount of water resaining efter each Finally, will sarple 'nessurenents(1] for each day fron nornal distribution with of [ater[i] and a standard deviation of 8.81, chesen because 8.81 kg 18g.
variable 'evaporation" which measures traction
disappears
MODEL
data
int ney,
array[nun.ce days]
reslclowers0, upper=l> eveporation;
trans formed paraneters(
array (nys) realslower=0, upper=1> water;
real current,
forl1 in [
current_water (l-evaporation):
vater(i) wetrent.
nodel {
evaporation
forl1 in days)
nornel [i], 0.01);
60

---

NeurIPS Paper Checklist
Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The abstract and introduction reflect the contribution and scope of the paper, and do not include aspirational goals.
Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification:
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [NA]
Justification: Yes, typically in appendix.
4 Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: All details disclosed. 5. Open access to data and code
Question: Does the paper provide open access the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [No]
Justification: It is difficult to provide code runnable by third party given the usage of local cluster and many non-portable modifications. Every effort has been made to make the results reproducible from the given description.
Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, sypea--arameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification:
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification:
8. Experiments compute resources
Question: For each experiment, does the paper provide sufficient information the computer re- sources (type of compute workers, memory, time execution) needed to reproduce the experiments?
Answer: [Yes]
61

---

Justification:
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurIPS Code of Ethics?
Answer: [Yes]
Justification:
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification:
This paper does appear to present significant social impacts beyond the core goal (which is discussed) of making probabilistic models easier to write.
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification:
13. New assets
Question: Are new assets introduced in the paper well documented and the documentation provided alongside the assets?
Answer: [NA]
Justification:
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Does the paper deseribe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based the requirements your country or institution) were obtained?
Question:
Answer: [NA]
Justification:
16. Declaration of LLM usage
Question: Does the paper describe the usage of L.LMs it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing,
62

---

editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required
Answer: [Yes]
Justification: It was used only for generating probabilistic models described clearly in the text. LLMs were not used in any other fashion, e.g. for generating code or writing.
63


Reviews:
[{"review_id": "4nQOEJMdir", "reviewer": "Reviewer_eRLK", "summary": "This work deals with a key limitation of probabilistic programming: model specification. Domain experts, such as clinicians, are typically not familiar with the syntax and semantics of probabilistic models, posing a significant barrier to quantitative analysis and requiring a dedicated data scientist or engineer to encode their domain knowledge into the model. Knowledge elicitation is often error-prone and cumbersome. This work addresses this problem using LLMs, allowing experts to describe their problems in English and translating the problem description (t) into a probabilistic program (m) using a decoder-only LLM, such as GPT. The program m describes a model over target variables (z) given observed variables (x), that is, p(z | x, m). Since LLMs are stochastic, the program m can be interpreted as having been sampled from p(m | t). The authors exploit this fact to make the program more robust by defining it as an ensemble, a set of sampled programs, the final density over the target a weighted average of each model’s density. They consider two methods: flat averaging and Bayesian model averaging (BMA). While the former weighs each sampled model equally (that is, proportional to the prior distribution p(m | t)), the latter uses the posterior over the model $p(m | x, t) \\propto p(m | t) p(x | m)$, which has the additional p(x | m) marginal likelihood term. Exact inference is exceedingly difficult. First, the LLM’s distribution over the model m given the text description p(m | t) is usually inaccessible; only samples drawn from it are accessible. Second, the model’s posterior density over the target p(z | x, m) and its marginal likelihood p(x | m) are hard to compute exactly. The former is estimated via NUTS, while the latter is estimated via variational inference, setting the variational distribution to be a Gaussian matching the empirical mean and variance of the samples. The authors demonstrate the efficacy of their proposed method on five newly created domains: the fully Bayesian averaged model outperforms the flat average, especially in challenging domains.", "questions": "1. Since the models are LLM-generated and the experts are not familiar with probabilistic programming, there seems to be a risk of the model encoding undesirable biases from the LLM’s training data.\n\n2. Section 5.1 says that \"These results underline the importance that the LLM-defined prior be relatively broad.” Can you elaborate on this point? Would imposing syntactic validity constraints on the LLM output make it less likely to assign low p(m | t) to models having high marginal likelihood ?", "limitations": "yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths\n1. This framework exploits the complementarity between LLMs and probabilistic programs: LLMs excel at natural language processing but cannot reason; probabilistic programs excel at reasoning but cannot deal with unstructured information.\n2. Empirical evaluation on newly constructed domains\n3. Model-agnostic approximate inference procedure for Bayesian model averaging.\n\nWeaknesses\n1. Generated models failing to compile due to syntactic invalidity seems to be a major limitation. This seems to get worse in complex domains like Gold. This might be addressed by using constrained generation, e.g., Loula et al. (2023)\n2. While the authors construct new domains to avoid leakage with the LLM’s training data, the domains themselves are relatively less complex. It appears that the complexity of the domain, the size of the text description, and its conflicts with the LLM’s embedded information would affect the quality of models.\n\nReferences\n\nLoula, João, et al. \"Syntactic and Semantic Control of Large Language Models via Sequential Monte Carlo.\" ICLR, 2025.", "was_revised": true, "final_justification": "Based on the other reviews and the authors' response, I understand the manuscript better. I will keep my positive evaluation and increase my confidence to 3", "rebuttal": "Thanks very much for your review! We particularly appreciate the summary, which is, as far as we can see, accurate in every detail.\n\nGenerated models failing to compile is indeed a significant issue. Though, in practice it's not *as* severe as you might think. We take advantage of continuous batching (parallel inference) when generating models from the LLM which speeds things up enormously. The good news about invalid models is that they are quickly rejected during the compilation step, meaning no computational effort needs to be expended on them during inference. (See section F for some more discussion of this point.)\n\nWe have indeed considered using constrained generation. In fact, we initially were generating models using llama.cpp, which supports constrained generation (https://github.com/ggml-org/llama.cpp/blob/master/grammars/README.md), and we noticed that Stan has a published BNF language syntax (https://mc-stan.org/docs/reference-manual/syntax.html). The truth is that we ultimately did not to pursue this simply because llama.cpp did not suppose continuous batching and we could not find any LLM software that supported both continuous batching and constrained generation. As continuous batching made generation 1-2 orders of magnitude faster, we could not forego it.\n\nWe apologize for being somewhat unclear with the statement in section 5.1 that \"These results underline the importance that the LLM-defined prior be relatively broad.”\n\nWhat we were trying to say here is that the variance in Eq. 9 is determined by the chi-squared divergence between p(m|x,t) and p(m|t). To clarify notation, here p(m|t) is the conditional distribution over models from the LLM given only the model text, while p(m|x,t) would be the conditional distribution when you also condition on the observed data x.\n\nConsider, for example, a case where the LLM is extremely \"overconfident\" in that it assigns almost all probability mass to \"bad\" models m where the marginal likelihood p(x|m) is very low. In this case, it is conceivable that the posterior p(m|x,t) ∝ p(m|t) p(x|m) would be more influenced by the marginal likelihood p(x|m) then the LLM prior p(m|t) and so could concentrate in regions where p(m|t) is very low. If that were to happen, the chi-squared divergence between p(m|x,t) and p(m|t) would likely be very large, meaning that estimates from this method would have very high variance.\n\nFor this reason, we believe it is important that the LLM be instructed to be quite \"open minded\" about what a good model might be. This, of course, depends on the base LLM. But one can also try to encourage this by providing a good system prompt and a good set of examples, as we sought to do.\n\nImposing syntactic validity constraints would not have much impact on this issue. Since we reject all invalid models, the only impact of imposing those constraints would be to improve the fraction of valid models. Essentially, if you find 1000 syntactically valid models by generating and rejecting or imposing validity constraints, then the resulting accuracy would be the same.\n\nUltimately, we suspect that general progress in LLMs will ultimately improve the fraction of valid models, as it has done in the generation of other kinds of code (e.g. python code), enough that this would become a fairly minor issue. We also believe there is great opportunity for better inference methods that would make it practical to address more complex problems (see lines 304-307). We hope that this paper establishes the basic idea of defining posteriors directly from informal user text, and that with future progress in inference methods it will be practical to address large-scale problems."}, {"review_id": "4bzWfBu9L3", "reviewer": "Reviewer_u2vx", "summary": "This paper attempts to make Bayesian inference widely available by automatically\ndeducing the appropriate Bayesian model and performing inference under the model\nwith the user of the system only being required to provide some data and a text\ndescription of the latent variables that they are interested in.\n\nThe paper outlines a very straightforward recipe that can be very easily replicated.\nFirst, use an LLM to draw samples of valid PPLs and perform posterior inference over\nthe latent variable with each PPL. Then use these posterior samples to do variational\ninference of the marignal likelihood of the data given the PPL model. Finally, return \na likelihood-weighted average of the posterior from all the PPL model posteriors.\n\nThe paper also provides other methods where a pure variational inference approach can\nbe used.\n\nThere are some results provided on few examples with a few hundred data points\nincluding examples where the data type of the latent variable is not explicitly \nstated.", "questions": "Why is the proposed method in this paper not been compared to related work in BMA?\n\nWhy has the algorithmic variants in section 5.3 been provided when there are no\nexperiments to back up the claims here?\n\nWhat kinds of models would experts and non-experts have generated for these\nscenarios and how do the highest weight models compare to those?\n\nHow well calibrated are the posteriors?", "limitations": "yes", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 3, "originality": 2, "strengths_and_weaknesses": "## Strength : Novelty\n\nWhile there are other papers on PPL generation including generation of PPLs from LLMs\nthere are none that explicitly perform model averaging as described here in a simple\npractical recipe.\n\nThe authors have also provided a second inference method using a purely variational\napproach which might be a novel.\n\n## Strength : Broad Applicability\n\nThe ability to make Bayesian Inference broadly applicable with a simple text prompt\nhas the possibility to introduce an unprecedented level of rigor in scientific\ndiscourse where probabilistic inference is needed for an accurate interpretation of\ndata but where statistical expertise is generally lacking.\n\n## Weakness : Detached from BMA literature\n\nThis paper is a very natural continuation of the Bayesian Model Averaging literature\nbut it appears to attempt to negate that connection by a rather offhand remark,\n\"However, existing algorithms for Bayesian model averaging (e.g. reversible jump MCMC\n [8, chapter 11]) do not appear to be applicable here, since they involve explicitly \n iterating over all models (here very larger or possibly infinite) and/or using evaluations of p(m|t).\"\n\nThis is not a fair remark. In BMA literature model averaging does use the posterior \nof the model given the data as the weight which is very similar to the current paper.\nCurrent paper draws samples of the model using an LLM and then weights these samples\nusing the likelihood of the data given the model. Now, depending on how one chooses\nto interpret the samples from the LLM this can be considered technically identical.\n\nIn any case, there is a fair amount of literature in BMA including frequentist methods\nwhich could all be considered relevant. It is unclear if the two methods described in\nthis paper are indeed novel unless a thorough comparison is made to BMA.\n\n## Weakness : Scalability\n\nNow, PPL inference is already quite intractable and this paper adds another layer\nof intractability by generating a large number of PPLs and performing inference on \nall of them. This makes the approach appear to be impractical.\n\nThe paper demonstrates results on only very simple examples with a few hundred data\npoints at best. This doesn't provide any reassurance on the scalability of the approach.\n\n## Weakness : Baseline Comparisons\n\nThe experiments don't really compare the generated posteriors to those that would\nbe generated if say an expert was to write a PPLs. It would be fine to compare\nagains the performance of a PPL generated by a non-expert.\n\nThere is no discussion of whether the generated posteriors are in any way calibrated\nwhich makes it hard for the reader to convince themselves that the proposed method\nhere is sound.", "was_revised": true, "final_justification": "The authors have failed to place their paper in context with other BMA algorithms under a factually incorrect pretext. The authors claim that all other BMA algorithms explicitly iterate all models or require the ability to compute a model's prior probability. This claim is false. There are many decades of BMA algorithms that the authors are glossing over.\n\nWithout properly placing this paper in the context of other BMA algorithms, it is not possible to assess the novelty of this work. As such, we have to assume that this work is not novel.", "rebuttal": "Thanks very much for your review!\n\nRegarding this quote:\n\n> \"However, existing algorithms for Bayesian model averaging (e.g. reversible jump MCMC [8, chapter 11]) do not appear to be applicable here, since they involve explicitly iterating over all models (here very larger or possibly infinite) and/or using evaluations of p(m|t).\"\n\nWe believe this remark is correct as written, although we concede that we weren't as clear as we could have been.\n\nWe certainly agree that the final posterior can be seen as an instance of Bayesian model averaging, and we tried to be explicit about that (See lines 50-51).\n\nHowever, there is a crucial distinction. In the setup of this paper, we assume that it is not possible to evaluate p(m|t). We briefly allude to this on lines 75-79, but we should be more explicit as it is a subtle point.\n\nThere are two issues here:\n\nFirst, most commercial LLM providers (e.g. OpenAI) simply do not provide probabilities for their generations—probably to make it harder for other companies to distill their models.\n\nSecond, for open-source LLMs (like llama, which we use) generation probabilities are available. However, in order to encourage the LLM to generate better models, we ask it to generate a \"chain of thought\" where it informally describes the modeling strategy. (See the examples in Section G.) We have found this to be crucial for generating high quality models. But the impact is that the generation probabilities for the model (alone) are not available. Rather, the generation probabilities give you something like p(informal text, formal model). Finding the probabilities of the model (alone) is completely intractable as it would require marginalizing out all possible informal texts.\n\nRegarding scalability, we must be straightforward that we simply do not intend to claim great scalability. The current inference method is, frankly, extremely expensive, as it requires doing MCMC on ~1000 models. (See section F for some discussion of timing.) However, we believe that there is great opportunity for faster inference methods (see lines 303-307) that would greatly increase scalability. At the same time, we respectfully suggest that given that this is the first method to directly define this kind of posterior directly from user text, it is useful to establish the basic sanity of the framework, and that even a method that can only (currently) address modestly-sized problems constitutes a valid contribution.\n\n> Why has the algorithmic variants in section 5.3 been provided when there are no experiments to back up the claims here?\n\nWe have given a series of increasingly complex algorithmic variants so that it is possible to discuss and analyze the different algorithmic components separately. Our final algorithm (\u0010given explicitly as Alg 7, section D) uses all the enhancements discussed in this section, and these are crucial for making inference possible.\n\n> What kinds of models would experts and non-experts have generated for these scenarios and how do the highest weight models compare to those?\n\nWe were advised that doing this would be considered human subjects research and thus subject to extremely cumbersome legal and IRB requirements. Our hope was that the reader of the paper could largely play this role. For example, in Figures 6-8 (appendix B.1) it is shown that for the rain example, essentially all posterior weight is given to a single model, which we provided in the text as m^(742). Realistically, we believe it is reasonable to expect the reader of the paper to look at the given informal text and this model and judge for themselves if the final model model should be considered reasonable. The same is true for the other models.\n\nHow well calibrated are the posteriors?\n\nWe ask for a bit more detail on what is requested here. (We may be able to provide it.) Typical notions of calibration (as, e.g., in a brier score decomposition) require doing inference a number of times and comparing the predictions to observed outcomes. This would not be practical given the expense of this method. For the city temperature problem, we actually did provide some small amount of evidence in Figure 30."}, {"review_id": "EPzTCjTrzw", "reviewer": "Reviewer_zbWh", "summary": "This paper proposes a simple but powerful way to turn a plain-language problem description into a full Bayesian analysis without hand-coding a model.\nAn LLM is prompted to emit many candidate formal models in a probabilistic-programming language; each model’s posterior and marginal likelihood are approximated.\nThe authors give a practical algorithm, theoretical justifications and experiments on synthetic dataset to support their proposed idea.", "questions": "1. Is my understanding of \"only need to match target\" correct?\n2. Can the authors provide some discussion on how can one apply their method in practice?\n3. How does the computational cost behave when applying the proposed algorithm?", "limitations": "Yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "## Strength\n1. I think the paper is very well-written, it is very easy to follow what the authors wish to express.\n2. The paper provides a practical way for LLM users to perform principled Bayesian predictions, which leads to nice practical impact.\n3. Theoretically, one only needs to match the target variable, which in my opinion provides strong flexibility.\n\n## Weakness\n1. The experiments are done on synthetic data. While I would assume it is not easy to find a real world data in this case, can the authors provide some discussion on how can one apply their method in practice?\n2. Typo: line 9 (justified \"an\" analyzed).", "was_revised": true, "final_justification": "The authors have addressed my concerns, I remain my score.", "rebuttal": "Thanks for your review!\n\n> Is my understanding of \"only need to match target\" correct?\n\nIt is true that the user only needs to specify the *target* variables of interest. While running, the LLM can and does invent other random variables. This is extremely convenient since most applied statistics problems involve, e.g., hierarchical modeling and partial pooling. But here the LLM takes care of all that for you. (Please let us know if we misunderstand the question.)\n\n> Can the authors provide some discussion on how can one apply their method in practice?\n\nSure. In practice, all the user needs to do is provide a model text (as a .txt file) and a data file (in our implementation, as a .json file). Then, our implementation uses a custom system prompt and sequence of example generations to generate a set of candidate Stan models (See sections F.1, F.2, and G). Then, each of these models is compiled into a Stan executable and (if it compiles) MCMC inference is run. Finally, those MCMC samples are used to construct a Gaussian distribution. That Gaussian distribution is then used for importance-weighted VI to estimate a marginal likelihood. Finally, all the MCMC samples are weighted by the marginal likelihoods and returned to the user.\n\nPractically speaking, running this algorithm is a challenge. We run the LLM on a cluster with an A100 GPU using continuous batching. The generations are saved as text files. Then, separate (CPU) jobs are submitted for each to run MCMC and estimate the marginal likelihoods. The final combination is easy and can be done on any local computer.\n\n> How does the computational cost behave when applying the proposed algorithm?\n\nPlease see section F for specifics on this point. Note also that we also give some ideas for future directions to improve efficiency on lines 303-307."}, {"review_id": "i8SJPwhbTb", "reviewer": "Reviewer_ovpx", "summary": "The paper proposes a procedure by which to use LLMs alongside Stan or another PPL for posing probabilistic inference problems: prompt the language-model with an informal description of the generative model and the target variables for posterior prediction, have it generate probabilistic programs, and then ensemble those probabilistic programs to predict the target variables. They demonstrate that this can work on some example problems.", "questions": "Can the authors systematically characterize which varieties of approximate inference they tried out in order to get the scheme we see in this manuscript?  How do the experimental results vary if hyperparameters or the choice of Monte Carlo algorithms are varied?", "limitations": "I would like to see an answer to my questions above to really get a sense of whether the inference results shown are robust enough to \"deploy\" to people who are not PPL full-timers.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 2, "significance": 4, "originality": 4, "strengths_and_weaknesses": "Strength: manages to still research Bayesian inference while making use of the latest computational capabilities\n\nWeakness: model averaging as a weighting strategy can suffer from high variance and unstable weighting in the case of model misspecification (see: https://proceedings.mlr.press/v238/reichelt24a.html), which is almost sure to happen when converting from natural language to probabilistic program via LLMs.  Self-normalized importance sampling similarly has bias and variance issues of its own.\n\nStrength: the authors created a fresh benchmark of probabilistic inference problems instead of using standard ones that could be in training corpora.", "was_revised": true, "final_justification": "I'm impressed with the author's response to my question and the simplicity of their model averaging method, so I've raised from score from Borderline to Accept.  The computational load of such a procedure is of course quite large, but after all you don't train or use LLMs when you're compute-poor.", "rebuttal": "Thanks for your review!\n\n> model averaging as a weighting strategy can suffer from high variance and unstable weighting in the case of model misspecification\n\nWe believe these issues are captured by our analysis in section 5.1.In particular, see the expression for the asymptotic variance of the estimate in Eq. 9. (Or a more careful analysis in C.2) This variance is determined rather the divergence between the LLM prior p(m|t) and the posterior p(m|x,t). So the main issue, as we see it, is if the LLM prior p(m|t) is reasonably broad. We allude to these issues on line 228-231, though it may be appropriate to give this discussion more emphasis. (Please also see our response to eRLK for more discussion of this point.) \n\n> Can the authors systematically characterize which varieties of approximate inference they tried out in order to get the scheme we see in this manuscript? How do the experimental results vary if hyperparameters or the choice of Monte Carlo algorithms are varied?\n\nSure. For the MCMC step, we immediately settled on using Stan's default algorithm with no modifications. We never experimented with any variants of this. The only question was how to estimate the marginal likelihood. To do this, we experimented a little bit with doing importance-weighted VI using normalizing flows. However, we found that this was a mess, as it's very difficult to make stochastic optimization over 10,000 models fully automatic. Thus we discarded it and switched to the current method.\n\nThe final method, while pretty expensive, is absolutely bulletproof and has *very* few hyperparameters to set. All that happens for each model is:\n\n1. Do MCMC using Stan's default algorithm.\n2. Compute the mean and covariance of those MCMC samples.\n3. Create a Gaussian with that mean and covariance.\n4. Repeat N times: Draw M samples from the Gaussian, and estimate the importance-weighted ELBO.\n\n(See Algorithm 9 in section F)\n\nThere are literally no hyperparameters in this algorithm other than N or M. We used M=25 and N=10,000, simply on the basis that they seemed reasonable. We never experimented with tuning either of them.\n\nProvided that Stan's MCMC method succeeds, this method is extraordinarily robust. After creating it on small test problems, we only ran it once to generate all the results. There was no secret fiddling. There may be some bias to the marginal likelihoods if the posterior for some model is poorly approximated by a Gaussian and importance-weighted VI cannot correct this. However, notice that if this happens, the failure is \"gentle\"—the resulting ELBO will be lower than the marginal likelihood, meaning that model gets lower posterior weight than it deserves. In this way, the method is very \"safe\"."}, {"review_id": "sdTa8MPzXY", "reviewer": "Reviewer_SVR6", "summary": "This manuscript introduces Large Language Bayes (LLB), a method for automating Bayesian modeling from natural language descriptions using a large language model (LLM).\nThe primary objective of the paper is to generate multiple candidate models based on a user's prompt. Then, a final prediction is computed by creating a weighted average of the posteriors from each model. The weight of each model is determined by its marginal likelihood, which measures how well observed data is fitted.", "questions": "1. When sampling PPL models, are the output generations individualized, meaning they each have their own “message thread”?\n    * If yes, how did the authors deal with the LLM generating/suggesting the same model multiple times? Any statistics on this? Wouldn’t this behavior possibly be a bias built into the LLM? Do you have any thoughts on how to handle this when writing the model with the PPL?\n    * If not, wouldn’t prior model output influence subsequent ones? That is, creating a chain of biases.\n2. How would this method compare with traditional Bayesian learning, where the optimization problem is to find the model parameters from a distribution that could, potentially, have a prior?\n3. Any thoughts on using the approximate ~p(z|x,m) as a score for dynamically prompting the LLM? In an iterative approach, one could use the posterior approximation as a guide for the LLM prompt template, aiming for more accurate models.", "limitations": "- Reproducibility of results is dependent on LLM size/quality and prompting\n- Inference demands limit the method to smaller use-cases", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 2, "originality": 2, "strengths_and_weaknesses": "- Strengths\n    - The paper presents an accessible tool for Bayesian modelling using LLMs through a theoretical analysis.\n    - Experiment results are encouraging.\n- Weaknesses\n    - LLM model sampling appears to be highly sensitive to multiple factors, including prompting, sequencing, and model. The manuscript could mitigate practical concerns related to sampling sensitivity with further experiments. (See Questions).\n    - The result in Theorem 2 - addressing the posterior approximation quality - isn’t feasible in practice, given the ample search space of models.", "was_revised": true, "final_justification": "The authors have addressed my concerns in their rebuttal, although the response didn't address my point with the formality it would require. However, after reading other reviewers' points of view, I believe my concern is outweighed by some of the good strengths highlighted by my peers. I am changing my recommendation to a weak accept to reflect the discussion in the rebuttal phase.", "rebuttal": "Thanks for your review!\n\n> When sampling PPL models, are the output generations individualized, meaning they each have their own “message thread”?\n\nThis is a good question. Yes, each sampled model had its own \"thread\", as mathematically required for model to be a valid IID sample from p(m|t).\n\n> If yes, how did the authors deal with the LLM generating/suggesting the same model multiple times? Any statistics on this? Wouldn’t this behavior possibly be a bias built into the LLM? Do you have any thoughts on how to handle this when writing the model with the PPL?\n\nIf the LLM wanted to suggest the same model multiple times, it was free to do so. We simply treated those as separate samples and processed them independently. This is correct, mathematically, but is potentially computationally wasteful. (It's silly to run MCMC multiple times on the same model.)\n\nNote that it's extremely rare for the LLM to generate exactly the same model. However—unfortunately—it's possible to create *equivalent* models, and these are somewhat challenging to detect. For example, equivalent models might use different variable names, or declare those variables in a different order. For these reasons we unfortunately can't easily offer any hard statistics. However, our qualitative impression is that this happens fairly often for the very simple models (e.g. coins) with the more complex models, a large fraction of the generations are unique. (This can be seen, to some degree, by looking at the distribution of estimated marginal likelihoods in the detailed experimental results in section B.)\n\n> How would this method compare with traditional Bayesian learning, where the optimization problem is to find the model parameters from a distribution that could, potentially, have a prior?\n\nWe're very sorry, but we can't quite understand this question. We're confused by the description, which doesn't correspond to traditional Bayesian learning as we understand it. We'd like to reply, so we'd be grateful for a bit more detail on what's being asked here.\n\n> Any thoughts on using the approximate ~p(z|x,m) as a score for dynamically prompting the LLM? In an iterative approach, one could use the posterior approximation as a guide for the LLM prompt template, aiming for more accurate models.\n\nAnother good question! We have considered several variants of this idea, but haven't validated one yet.\n\nHere's one option. Imagine you generate one model at a time and do inference on it before generating the next model. One could conceivably take some of the top scoring models and provide them as additional examples as part of the prompt. Done naively, this would, of course, bias the results. However, one could potentially correct this with importance sampling—after the chain of thought text, generate the model with the extra examples, compute the probability of that generation, and then also compute the probability using the LLM without the extra example models. In this way, one could hope that the model would converge to the true results more efficiently. We caution that we haven't actually done this yet. But we believe this general kind of strategy (along with the ideas we allude to on lines 303-307) may greatly improve the method."}]

Author rebuttal:

Meta-review: This article combines large language models (LLMs) and probabilistic programming to define a joint distribution over models, latent variables, and data. The LLM is used to sample candidate models from informal descriptions, while the probabilistic program performs per-model posterior inference and marginal likelihood estimation. The approach instantiates Bayesian model averaging (BMA): the posterior is approximated via self-normalized importance sampling over models, relying on (i) sampling models with the LLM, (ii) estimating per-model posteriors using Stan or similar, and (iii) approximating per-model marginal likelihoods via variational inference. Experiments are conducted on synthetic data.

After discussion, the paper remains borderline, with a majority leaning toward acceptance. 

Strengths: 

* Clear writing.
* Novel combination of LLMs with probabilistic programs.
* Model-agnostic approximate inference.
* Generality of the framework.

Limitations:

* Experiments are limited to synthetic data.
* Computational cost is high, limiting applicability.
* The importance sampling-based BMA is relatively simple and may suffer from high-variance estimates.
* Some inaccuracies in the discussion of the BMA literature.

The authors’ response satisfactorily addressed most concerns; however, issues regarding the treatment of the BMA literature remain and should be corrected in the final version.
Decision: Accept (poster)