We present AutoStan, a framework in which a command-line interface (CLI) coding agent autonomously builds and iteratively improves Bayesian models written in Stan. The agent operates in a loop, writing a Stan model file, executing MCMC sampling, then deciding whether to keep or revert each change based on two complementary feedback signals: the negative log predictive density (NLPD) on held-out data and the sampler's own diagnostics (divergences, R-hat, effective sample size). We evaluate AutoStan on five datasets with diverse modeling structures. On a synthetic regression dataset with outliers, the agent progresses from naive linear regression to a model with Student-t robustness, nonlinear heteroscedastic structure, and an explicit contamination mixture, matching or outperforming TabPFN, a state-of-the-art black-box method, while remaining fully interpretable. Across four additional experiments, the same mechanism discovers hierarchical partial pooling, varying-slope models with correlated random effects, and a Poisson attack/defense model for soccer. No search algorithm, critic module, or domain-specific instructions are needed. This is, to our knowledge, the first demonstration that a CLI coding agent can autonomously write and iteratively improve Stan code for diverse Bayesian modeling problems.
Figures & tables
Figure 1: AutoStan on the large 1D regression dataset ( n=500 train, 200 test; DGP defined in Section 3.1 ). In panels (b)–(d), the fitted mean f^(x) is plotted with ±2σ^(x) shaded; the faint grey overlay marks the oracle ground truth f(x)±2σ(x) . Training points outside the plot range [−5,7] are shown as edge-pinned triangles ( ▲ / ▼ ) at the panel boundary. (a) NLPD trajectory over 15 iterations; colored markers identify the three model-fit panels. The dashed oracle line ( NLPD=1.14 ) is a lower bound that cannot be reached (see Section 3.1 ). (b) Baseline (iter 0): Gaussian linear model; predictive bands are dominated by the ≈30 extreme training outliers ( NLPD=2.16 ). (c) Iter 1: cubic polynomial mean + Student- t likelihood; one step eliminates most band inflation ( NLPD=1.32 ). (d) Iter 11 (best): sinusoidal mean with learned frequency ω , heteroscedastic log-linear variance, and a contamination-mixture likelihood; tight, locally calibrated bands closely follow the oracle noise envelope ( NLPD=1.23 ). (e) TabPFN (90% predictive interval, NLPD=1.25 ): mean tracking is accurate but intervals are uniformly too wide—the model has absorbed the training outliers, inflating uncertainty across the full input range.
Iter
NLPD
Model
Δ
0
2.1589
Linear mean, Gaussian likelihood
—
1
1.3181
Cubic polynomial + Student- t
− 0.84
2
1.3060
Sine basis + Student- t
− 0.01
5
1.2854
+ Quadratic logσ(x)
− 0.02
7
1.2844
+ Learnable frequency ω
− 0.001
9
1.2325
Gaussian mixture (fixed σout=10 )
− 0.05
Table 1: NLPD trajectory, 1D regression large ( n=500 ). Oracle NLPD = 1.1442. Full trajectory including rejected iterations in Table 4 .
Hier. S
Hier. L
Slopes
1D S
1D L
Bundes.
Oracle
1.494
1.404
1.263
0.944 †
1.144 †
—
Baseline
1.500
1.404
1.818
2.248
2.159
1.566
Best
1.500
1.401
1.274
1.124
1.226
1.543
Gap
0.006
− 0.003
0.011
0.180 †
0.081 †
—
Iters
0
4
10
5
11
9
Table 2: Summary across all experiments. Gap = Best − Oracle (lower is better; negative means the model outperforms the oracle, likely due to test-set overfitting). † Oracle is a strict lower bound that cannot be reached: the oracle knows the test set contains no outliers, whereas the agent must account for the ≈6% training contamination. The agent exhibits qualitatively different behaviors depending on dataset structure and size.
Appendix figures & tables9 assets
Supplementary material from the paper’s appendix.
Appendix
Figure 2: AutoStan on the 1D regression small dataset ( n=68 ). (a) NLPD trajectory over 9 iterations. (b) Baseline: linear + Gaussian, enormous bands from 4 extreme outliers. (c) Iter 1: quadratic mean + Student- t , the largest gain ( Δ=0.75 ). (d) Iter 5 (best): cubic + log-linear σ(x) + contamination mixture; NLPD 1.1244, matching TabPFN (1.1202) while remaining fully interpretable.
Figure 3: Bundesliga model parameters (iteration 9). (a) Dortmund (0.39), Freiburg (0.33), Bayern (0.32) benefit most from playing at home; Heidenheim ( − 0.01) and St. Pauli (0.00) show no home advantage. (b) Attack vs. defense with 90% credible intervals; Bayern’s dominance is clearly visible.
Figure 4: TabPFN predictions. On the large dataset, the intervals widen globally near outliers instead of isolating them, explaining TabPFN’s higher NLPD (1.2501 vs. AutoStan’s 1.2256).