cs.LGOct 19, 2026

Bayesian Model Discovery with Language Models: Closing the Loop between Hypotheses and Experiments

Authors: JooKyung Yuh, Byungmin Park, Eunseo Choi, Jinheon Baek, Junmo Cho, Kyu Hyun Choi, Sung Ju Hwang, Kyong Hwan Jin, 

Organizations: Korea University · KAIST · KETI

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

Abstract

Scientific model discovery requires proposing candidate models of an unknown system and designing experiments that discriminate among them. However, LLM-based approaches typically rely on in-context learning to predict outcomes and select experiments, lacking explicit model evaluation to systematically guide hypothesis search, while Bayesian approaches enable principled model comparison but are confined to a predefined candidate set which results in limited search spaces. To this end, we introduce BMD, which closes the loop between LLM-driven hypothesis generation and Bayesian experimental design. Specifically, an LLM proposes model families, while a numerical engine fits their parameters, weighs them by model evidence, and averages their predictions. Then, BMD chooses the next experimental conditions by maximizing the joint expected information gain over families and their parameters, while feeding the updated beliefs and residual diagnostics back to the LLM to guide its next proposals. Across four scientific discovery benchmarks, BMD achieves lowest final prediction error than in-context learning and LLM-based scientific agents while recovering exact functional forms more often, with ablations supporting numerical feedback and the joint objective.

Figures & tables

Figure 1: Overview of Bayesian Model Discovery. Given the observations Dt\mathcal{D}_{t} , the numerical engine fits each candidate family and returns the resulting beliefs and residual diagnostics to the LLM, which proposes new families. These proposals are validated, fitted, and compared with the existing families, after which the updated pool selects the next experimental conditions xt+1\mathbf{x}_{t+1} by maximizing the joint expected information gain.
BenchmarkScientific settingDesign →\rightarrow observationEvaluation target
ChemBenchEnzyme kineticsReaction conditions →\rightarrow rateNoise-free observation mean (nMSLE)
BoxingGymEcological/statistical processesTask inputs →\rightarrow noisy responsesNoise-free response mean (nMSE)
ForceBenchForce and motion lawsConditions →\rightarrow trajectoriesShifted-design trajectories (nMSE)
NewtonBenchCounterfactual physical lawsSettings →\rightarrow physical responsesNoise-free governing law (nMSLE)
Table 1: Overview of the four scientific discovery benchmarks.
Figure 2: Predictive performance. Lower is better. Lines show geometric-mean prediction error, aggregated over seeds within each task and then across tasks. Shaded bands indicate seed dispersion in log space. The x-axis shows the total number of observations. AutoSciLab curves start later because additional observations are required before its first executable model is available. Appendix A details the evaluation protocol, aggregation, and exclusions.
Figure 3: Strict recovery and ground-truth containment. Each curve shows the fraction of runs in which the selected top-1 model satisfies the corresponding criterion. Strict recovery (solid lines) requires structural equivalence to the ground-truth functional form. Containment (dotted lines) counts more general model families that can represent the ground truth as a special case. For example, a quadratic family contains a linear family when its quadratic coefficient is set to zero. Thus, containment does not necessarily imply recovery of the exact functional form. Appendix C details the evaluation criteria, including numerical checks.
 ComponentsRelative error ( A=1A=1 , ↓\downarrow )GT identification (%, ↑\uparrow )
ModelExplicitSMC fitEvidenceChemBoxingForceNewtonOverallChemBoxingForceNewtonOverall
AXXX1.0001.0001.0001.0001.000–––––
BOXX0.9170.8680.2300.7770.61398.918.550.090.079.3
COOX0.3511.3880.1740.2190.36998.918.550.090.079.3
DOOO0.3490.6670.0360.2470.21496.763.083.378.684.8
Table 2: The role of numerical model evaluation. The methods progressively introduce (B) explicit model-based prediction using an LLM-selected family and parameters, (C) Bayesian parameter fitting via SMC for the family selected by B, and (D) evidence-based selection instead of LLM selection. Prediction errors are normalized by direct ICL (A) within each benchmark; Overall is the geometric mean across benchmarks.
Figure 4: Exploring the model pool. Rows follow first-admission order; shading gives normalized evidence weights. Red marks ground-truth-equivalent models, blue other active models, and white absence from the pool. Black dots mark the highest-weight model. Annotations show candidate functional forms; X1X_{1} is the anonymized input, and letters denote coefficients.
Figure 7
Appendix figures & tables9 assets

Supplementary material from the paper’s appendix.

Appendix

 ChemBenchBoxingGymForceBenchNewtonBench
TasksMichaelis–Menten, product inhibition, ping-pong; each easy/medium/hardDeath, dugongs, hyperbolic, peregrines, survivalGravity, Yukawa, Coulomb (easy), oscillator, extra dimensionsMagnetic force, Fourier, Snell, radioactive decay, Malus, sound speed, Hooke, Bose–Einstein, heat transfer (medium)
BMD runs27151527
Recovery runs2712924
n0n_{0} / additional / final nn3 / 30 / 332 / 30 / 322 / 16 / 184 / 30 / 34
Design space7 continuous reaction conditionsTask-specific; hyperbolic has 3 integer inputs, survival a discrete patient menu2 continuous parameters in [0.1,10]2[0.1,10]^{2}Task-specific continuous inputs
Prediction scoredAfter proposalsBefore proposalsBefore proposalsAfter proposals
Appendix
Table 3: Main experimental tasks and budgets. Structural recovery uses only targets representable in the interface; predictive evaluation retains all tasks. Each task has three seeds.
 ChemBenchBoxingGymForceBenchNewtonBench
Initial data3 fixed designs2 random designs2 trajectories4 random designs
Coefficient priorDefaultDefault; hyperbolic: wideSigned N(0,1)\mathcal{N}(0,1) ; positive LogNormal (0,1)(0,1)Wide
Observation familyFixed LogNormalFixed per taskFixed trajectory likelihoodFixed LogNormal
Archive reviewEvery 5 rounds to 30Every 5 rounds to 30Rounds 5, 10, 15Every 5 rounds to 30
EIG draws128 posterior draws128; count tasks use outcome tables128 outer and 128 inner draws128 posterior draws
Search starts5 uniform5 best candidate designs; survival enumerated5 uniform5 uniform
Appendix
Table 4: Benchmark-specific settings for the main BMD runs. Common proposal and SMC settings are given in the text. Appendices E and F specify analysis-specific differences.
Benchmark (metric)ABCD
ChemBench (nMSLE)8.47 ⁣× ⁣10−48.47\!\times\!10^{-4}7.77 ⁣× ⁣10−47.77\!\times\!10^{-4}2.97 ⁣× ⁣10−42.97\!\times\!10^{-4}2.96 ⁣× ⁣10−42.96\!\times\!10^{-4}
BoxingGym (nMSE)4.17 ⁣× ⁣10−24.17\!\times\!10^{-2}3.62 ⁣× ⁣10−23.62\!\times\!10^{-2}5.79 ⁣× ⁣10−25.79\!\times\!10^{-2}2.78 ⁣× ⁣10−22.78\!\times\!10^{-2}
ForceBench (nMSE)1.00 ⁣× ⁣10−31.00\!\times\!10^{-3}2.30 ⁣× ⁣10−42.30\!\times\!10^{-4}1.74 ⁣× ⁣10−41.74\!\times\!10^{-4}3.58 ⁣× ⁣10−53.58\!\times\!10^{-5}
NewtonBench (nMSLE)9.34 ⁣× ⁣10−69.34\!\times\!10^{-6}7.26 ⁣× ⁣10−67.26\!\times\!10^{-6}2.05 ⁣× ⁣10−62.05\!\times\!10^{-6}2.31 ⁣× ⁣10−62.31\!\times\!10^{-6}
Overall (geomean)7.58 ⁣× ⁣10−47.58\!\times\!10^{-4}4.65 ⁣× ⁣10−44.65\!\times\!10^{-4}2.80 ⁣× ⁣10−42.80\!\times\!10^{-4}1.62 ⁣× ⁣10−41.62\!\times\!10^{-4}
Appendix
Table 5: Absolute prediction errors for the shared-candidate study in Table 2 , before normalization by A. Entries are geometric means over seeds within worlds and then over worlds; the overall row is the reported geometric mean across benchmarks. Lower is better.
BenchmarkB/CD
ChemBench89/9087/90
BoxingGym5/2717/27
ForceBench15/3025/30
NewtonBench63/7055/70
Overall172/217184/217
Appendix
Table 6: Top-1 GT identification in the fixed-pool study (correct / evaluated runs). B and C select the same family. A has no structural output.
BundleContentsVisible in
CommonObservations, task interface, candidate expressions, coefficient roles, grammar, structural duplicate history, syntax/repair status, and eligible retired identifiersAll variants
BeliefLog evidence, model weights, best-model marker, parameter posterior summaries, pre-observation predictive log scores, and archived evidence valuesb=1b=1
DiagnosticsStandardized predictive residuals for every active model at the acquired observations, identified by modeld=1d=1
Appendix
Table 7: Information shown to the proposer in the feedback ablation. Numerical inference remains active in all four variants.
ObjectiveChemBenchBoxingGymForceBench
Joint EIG1.47 ⁣× ⁣10−51.47\!\times\!10^{-5}0.01018.8 ⁣× ⁣10−88.8\!\times\!10^{-8}
Model-only1.93 ⁣× ⁣10−51.93\!\times\!10^{-5}0.01118.2 ⁣× ⁣10−68.2\!\times\!10^{-6}
Parameter-only1.63 ⁣× ⁣10−51.63\!\times\!10^{-5}0.01683.6 ⁣× ⁣10−63.6\!\times\!10^{-6}
Random2.98 ⁣× ⁣10−52.98\!\times\!10^{-5}0.05703.6 ⁣× ⁣10−53.6\!\times\!10^{-5}
EPIG—0.01329.5 ⁣× ⁣10−79.5\!\times\!10^{-7}
Appendix
Table 8: Earlier auxiliary acquisition comparison, separate from Figure 6 . Entries are final geometric-mean errors; ChemBench uses n=33n=33 , BoxingGym n=32n=32 , and ForceBench n=18n=18 . The feedback and initialization configuration is described in the text.
TaskPrimary EIGPrimary randomPrimary ICLObs. EIGObs. randomObs. ICL
c0 easy2.13 ⁣× ⁣10−42.13\!\times\!10^{-4}2.14 ⁣× ⁣10−42.14\!\times\!10^{-4}2.08 ⁣× ⁣10−42.08\!\times\!10^{-4}1.69 ⁣× ⁣10−61.69\!\times\!10^{-6}1.11 ⁣× ⁣10−51.11\!\times\!10^{-5}1.96 ⁣× ⁣10−41.96\!\times\!10^{-4}
c0 medium2.18 ⁣× ⁣10−42.18\!\times\!10^{-4}2.79 ⁣× ⁣10−42.79\!\times\!10^{-4}2.29 ⁣× ⁣10−52.29\!\times\!10^{-5}2.88 ⁣× ⁣10−62.88\!\times\!10^{-6}1.36 ⁣× ⁣10−51.36\!\times\!10^{-5}2.45 ⁣× ⁣10−42.45\!\times\!10^{-4}
c0 hard2.17 ⁣× ⁣10−42.17\!\times\!10^{-4}2.46 ⁣× ⁣10−42.46\!\times\!10^{-4}2.99 ⁣× ⁣10−52.99\!\times\!10^{-5}5.99 ⁣× ⁣10−65.99\!\times\!10^{-6}8.62 ⁣× ⁣10−68.62\!\times\!10^{-6}2.81 ⁣× ⁣10−42.81\!\times\!10^{-4}
c2 easy1.95 ⁣× ⁣10−41.95\!\times\!10^{-4}2.01 ⁣× ⁣10−42.01\!\times\!10^{-4}3.65 ⁣× ⁣10−43.65\!\times\!10^{-4}5.78 ⁣× ⁣10−65.78\!\times\!10^{-6}6.26 ⁣× ⁣10−66.26\!\times\!10^{-6}6.62 ⁣× ⁣10−56.62\!\times\!10^{-5}
c2 medium2.01 ⁣× ⁣10−42.01\!\times\!10^{-4}2.25 ⁣× ⁣10−42.25\!\times\!10^{-4}1.72 ⁣× ⁣10−41.72\!\times\!10^{-4}8.72 ⁣× ⁣10−78.72\!\times\!10^{-7}1.22 ⁣× ⁣10−51.22\!\times\!10^{-5}2.02 ⁣× ⁣10−42.02\!\times\!10^{-4}
c2 hard1.94 ⁣× ⁣10−41.94\!\times\!10^{-4}1.85 ⁣× ⁣10−21.85\!\times\!10^{-2}2.29 ⁣× ⁣10−42.29\!\times\!10^{-4}6.71 ⁣× ⁣10−66.71\!\times\!10^{-6}1.82 ⁣× ⁣10−21.82\!\times\!10^{-2}1.57 ⁣× ⁣10−41.57\!\times\!10^{-4}
Appendix
Table 9: Earlier ChemBench diagnostic nMSLE at n=63n=63 on 256 evaluation inputs. These are not per-task values for the current main figure.
  EntryStablePrimary error at admissionFinal error
PolicySeedroundroundBeforeAfterPrimaryObservation
EIG31998.18 ⁣× ⁣10−18.18\!\times\!10^{-1}1.64 ⁣× ⁣10−41.64\!\times\!10^{-4}2.72 ⁣× ⁣10−42.72\!\times\!10^{-4}2.90 ⁣× ⁣10−62.90\!\times\!10^{-6}
EIG37778.48 ⁣× ⁣10−18.48\!\times\!10^{-1}4.84 ⁣× ⁣10−44.84\!\times\!10^{-4}2.86 ⁣× ⁣10−42.86\!\times\!10^{-4}1.29 ⁣× ⁣10−51.29\!\times\!10^{-5}
EIG4118183.12 ⁣× ⁣10−33.12\!\times\!10^{-3}4.65 ⁣× ⁣10−54.65\!\times\!10^{-5}2.78 ⁣× ⁣10−42.78\!\times\!10^{-4}1.04 ⁣× ⁣10−51.04\!\times\!10^{-5}
Random31991.31 ⁣× ⁣10−11.31\!\times\!10^{-1}9.29 ⁣× ⁣10−59.29\!\times\!10^{-5}2.76 ⁣× ⁣10−42.76\!\times\!10^{-4}1.05 ⁣× ⁣10−51.05\!\times\!10^{-5}
Random375175.65 ⁣× ⁣10−15.65\!\times\!10^{-1}5.27 ⁣× ⁣10−55.27\!\times\!10^{-5}2.29 ⁣× ⁣10−42.29\!\times\!10^{-4}1.39 ⁣× ⁣10−51.39\!\times\!10^{-5}
Appendix
Table 10: Earlier ping-pong-easy fitted histories on the official 1,000-input set. Entry and stable dominance concern the broader primary-family classifier. All errors are nMSLE. Final values use n=63n=63 ; admission rounds count observations acquired after n0=3n_{0}=3 .
DifficultyControlHistoriesMinimumMedianMaximum
EasyConstants estimated277.20 ⁣× ⁣10−77.20\!\times\!10^{-7}1.20 ⁣× ⁣10−51.20\!\times\!10^{-5}6.02 ⁣× ⁣10−56.02\!\times\!10^{-5}
EasyConstants supplied271.22 ⁣× ⁣10−91.22\!\times\!10^{-9}8.21 ⁣× ⁣10−68.21\!\times\!10^{-6}4.08 ⁣× ⁣10−54.08\!\times\!10^{-5}
HardConstants estimated279.19 ⁣× ⁣10−39.19\!\times\!10^{-3}3.11 ⁣× ⁣10−23.11\!\times\!10^{-2}6.82 ⁣× ⁣10−26.82\!\times\!10^{-2}
HardConstants supplied273.85 ⁣× ⁣10−73.85\!\times\!10^{-7}9.54 ⁣× ⁣10−69.54\!\times\!10^{-6}4.44 ⁣× ⁣10−54.44\!\times\!10^{-5}
Appendix
Table 11: Previously completed supplied-primary-form controls on frozen c5 histories, evaluated on 1,000 inputs. All entries are primary nMSLE.
More in this topic

Scientific Discovery

Scientific discovery uses computational methods to find hypotheses, laws, mechanisms, and materials from data, literature, and experiments.

See all 86 papers

Explore similar work

CardsList

  1. Model Discovery Agent: LLM-assisted Bayesian experiment design for data-efficient discovery of mechanistic world models

    Aug 10, 2026Kevin MurphyDynamical SystemsAI Agents for Scientific Discovery

    Dept. Computer Science Univ. British Columbia, Canada.

  2. From Prompts to Posteriors: Can LLM Agents Perform Bayesian Inference in Context?

    Oct 21, 2026Linus Filbry, Luckeciano C. Melo, Yarin GalBelief Updating in LLMsBayesian Inference

    OATML, University of Oxford

  3. Large Discovery Models: Empirically-grounded Model-Based Open-Ended Search

    Aug 16, 2026Zhongwei Yu, Yan Song, Xue Yan +9Bayesian OptimizationLarge Language Model-Guided Optimization

    The Hong Kong University of Science and Technology (GZ) · University College London · AI Lab, The Yangtze River Delta +3