Which Method?

Recommendations and Defaults

We have run the most comprehensive simulation study to date on confidence interval methods for LLM eval data, and with few exceptions, we have validated these results on real eval data as well. Based on this research, we recommend the following methods for calculating confidence intervals in LLM eval contexts. We also performed a smaller simulation study on p-value methods, and report recommendations as well. These recommendations are implemented as the method="auto" defaults in evalstats, which adapt to data type and sample size based on which methods performed best in our simulations. Note that the contents of this section are subject to change as we gather more evidence for different situations.

Recommended Confidence Intervals Methods (By Data Type and Type of Comparison)

Recommended p-value methods by analysis type

Sample Size Guidance

N < 15: Do not report statistics. The CI width will be so large as to be meaningless, and the sample is extremely unlikely to be representative.

15 ≤ N < 30: Report statistics with a prominent warning. Results should be treated as exploratory only.

N ≥ 30: CIs and point estimates of true performance begin to become more reliable.

Simulation Study: Which CI method should I use in X situation?

We ran comprehensive simulation studies over eval-like distributions to compare methods for calculating CIs head-to-head— to the best of our knowledge, the most comprehensive test of CI methods ever conducted for small-sample AI evals data. We compared almost all known CI methods (and even a few ones invented with the help of AI) to see which ones are best for different data types, sample sizes, and estimands. We also measured the computational cost of each method, since some of the more complex methods can be very expensive to run at scale. Finally, we validated the results on real eval data by subsampling from real benchmark outputs, to double-check how well the simulation findings hold up in practice.

How to measure CI "performance"?

A confidence interval (CI) is only useful if it's trustworthy, tight, and fast to compute. We measure CI quality along three dimensions:

• Coverage: If you ask a method to build a 95% CI, it should actually contain the true value 95% of the time across many repeated experiments. A method with under-coverage (say, 88%) gives intervals that are too narrow and misleadingly confident. A method with over-coverage (say, 99%) is overly conservative—it works, but it may be too cautious. We target nominal coverage of 95% for all simulations.
• Width: Once a method achieves good coverage, narrower intervals are better. Width measures how "confident" the CI is: all else equal, the method that produces the tightest interval while still hitting 95% coverage wins.
• Cost: how long the method takes to run. Simple analytic formulas (like the Wilson interval for binary scores) compute in microseconds. Resampling methods like the bootstrap require thousands of resampled datasets and can take seconds or more, especially as the sample size increases. Cost matters when you're running evaluations over many models, prompts, or benchmarks simultaneously.

In practice, coverage comes first—a cheap, narrow interval that's wrong 20% of the time is useless. Among methods that achieve adequate coverage, we prefer narrower intervals, and among those, faster ones.

Methodology

We rigged realistic distributions, randomly sampled from them, then measured how close each CI method gets to the true value across repeated trials. Two estimands were tested: single-sample CIs (does the interval contain the true mean?) and pairwise difference CIs (does the interval contain the true difference?). Sample sizes: N = 10, 15, 20, 30, 40, 50, 60, 70, 80, 90, 100. Bootstrap and Bayesian draws: 10,000. Monte-Carlo trials: 2,000 per condition. We also validated on real LLM eval data by subsampling from 1) OpenEvals outputs across 15 unique model × benchmark pairs, for the single-run case, and 2) outputs of six models across nine benchmarks run 5 times per benchmark, for the multi-run case. Simulations are implemented in simulations/ folder in the project repository.

Distributions Tested

Across our synthetic simulations, we tested the following four data types (our simulations using real data test only the first two, binary and continuous):

• Binary scores ({0, 1}): Bernoulli with p ∈ {0.02, 0.05, 0.1, 0.3, 0.5, 0.7, 0.9, 0.95, 0.98} (9 scenarios)
• Continuous scores (0–1 floats): Beta families — uniform Beta(1,1), U-shaped Beta(0.5,0.5), right-skewed Beta(2,8), left-skewed Beta(8,2), moderate-skew Beta(2,5), extreme-right Beta(0.35,6), extreme-left Beta(6,0.35), near-boundaries Beta(0.3,0.3), near-center Beta(6,6), mixture of Beta(0.5,4) and Beta(5.5,1.2); plus a logit-normal (μ=−0.35, σ=1.35), zero-inflated (70% spike at 0, Beta(2,4) body), and one-inflated (70% spike at 1, Beta(4,2) body) (13 scenarios)
• Likert scores (integers 1–5): latent-normal model rounded and clipped to 1–5, covering uniform (high variance), skewed-low, skewed-high, bimodal (peaks at extremes), center-peaked, near-floor, near-ceiling, polarized, and flat-middle (9 scenarios)
• Grades (0–100 floats): clipped-normal scenarios — symmetric N(50,20), high-scoring N(75,15), low-scoring N(35,20), ceiling-heavy N(88,10), floor-heavy N(12,10), very-high N(92,7), very-low N(8,7), high-variance N(50,34); plus a 3-component truncated-normal mixture, a heavy-tailed t(df=3) distribution, a zero-spiked (40% spike at 0, N(45,20) body), and a hundred-spiked (40% spike at 100, N(65,18) body) (12 scenarios)

Single-Sample Coverage (target: 0.95)

The following image and table presents the results of the Monte Carlo simulation. The table breaks down the final average coverage rates at each sample size, with notes on performance. (Coverage rate = fraction of trials where the CI contains the true population mean. Higher is better; under-coverage means the method is over-confident.)

Coverage target is 0.95. ▼ marks under-target coverage. ★ marks preferred methods. Methods with a scope note (binary only / continuous [0, 1] only) are not applied to other data types; their numbers are therefore not directly comparable to general-purpose methods. Bootstrap-family methods are averaged across all data types. MCSE values are negligible (around ~0.005 or lower), so they are omitted from this display.

Single-Sample Coverage on Real LLM Eval Data (target: 0.95)

To complement the synthetic simulation, we also evaluated CI methods on real LLM evaluation data by randomly drawing samples from evaluation score outputs in the OpenEvals dataset, across 15 unique model × benchmark pairs, for variety. Coverage was measured by repeatedly subsampling from each benchmark's full score distribution and checking whether the CI contained the true (full-population) mean. This tests how methods behave on actual eval data shapes, which may differ from the idealized distributions in the synthetic study.

Coverage target is 0.95. ▼ marks under-target coverage. ★ marks preferred methods. Methods with a scope note (binary only / continuous [0, 1] only) are not applied to other data types. Results are averaged across 15 model × benchmark pairs and all sample sizes.

Pairwise Comparison Coverage (target: 0.95)

For paired comparisons, coverage measures whether the CI on the difference contains the true difference. This is the more practically important statistic for comparing two prompts or models.

Pairwise scenarios are parameterised by two factors that are swept jointly: the intraclass correlation (ICC), which captures how much variance is attributable to the input (prompt/example) rather than the model, and the standardised effect size (Cohen's d), which sets the true difference between the two candidates. The official test sweeps ICC ∈ {0.05, 0.15, 0.30, 0.50} and d ∈ {0.0, 0.2, 0.4} (the d = 0 condition measures Type I error rate; the non-null conditions measure coverage under a real effect). ICC is implemented differently per data type: a Beta-hierarchical model for binary scores, a Gaussian-noise perturbation of a Beta base for continuous scores, and a latent-normal model for Likert and grades. All pairwise scenarios also use the full expanded scenario suite, so the distributions tested are the same rich set described above.

Coverage target is 0.95. ▼ marks under-target coverage. ★ marks preferred methods. Results are averaged across all eval types, scenarios, and sample sizes. Binary-only methods (Newcombe, Tango, Bayesian paired/independent) are shown for diagnostic reference and are not directly comparable to general-purpose methods.

So Did Tango Score Interval Just Outperform the Bayesian Paired Method?

Based on the synthetic simulation alone, it certainly looks that way. In that study, Tango (0.955) edges out the Bayesian paired method (0.958) in coverage while being nearly 100× faster and without the catastrophic N≈200+ failure mode. The synthetic distributions, however, are drawn from idealized parametric families. The natural question is: does this hold on real LLM eval data, where item correlations, difficulty clustering, and lopsided score distributions are facts of life rather than simulation choices?

Pairwise Coverage on Real LLM Eval Data (Binary)

To check, we ran a pairwise coverage simulation on OpenEvals data, subsampling from pairs of models that share items for the same benchmark, across multiple benchmarks, and measuring whether each method's CI on the paired difference contained the true (full-population) difference. The results reveal a discrepancy that the synthetic study masked:

Coverage target is 0.95. ▼ marks under-target coverage. ★ marks preferred methods. Results are averaged across all benchmarks, model pairs, and sample sizes. Methods shown for all types (bootstrap, t-interval, etc.) are included for reference; in practice only binary-specific methods are applied to binary data.

Tango is still good, but it slightly under-covers at small N<50—so we can't entirely reproduce the strong result. Why?

Why does Tango under-cover at small N on real data?

When we looked more closely at which model pairs drive the under-coverage, two patterns emerge:

• Dominated pairs: one model scores near 0% on a benchmark while the other scores, say, 80%. The discordant pairs (inputs where the two models disagree) are too lopsided for Tango's score statistic to accumulate enough signal, and the CI ends up slightly too narrow.
• Jointly sparse pairs: both models score very low, e.g. both around 2%. Again, discordant pairs are too rare, and the CI ends up slightly too narrow.

These extreme situations are very common in our sample of OpenEval's benchmark, since we sample a diverse set of models at different parameter scales (e.g., a 1B model versus a 70B one), where the strong model dominates. Our synthetic simulation represented these extreme situations, but in the end, their results are averaged with more "well-behaved" pairs across a wide range of distributions.

In fact, we can see how it might be an artifact of the OpenEval sampling by running the same simulation on different LLM evals data—taking a single-run of the multi-run benchmark data we collected (explained below). When we do that, we see Tango again matches Bayesian paired:

CI Coverage vs. Computation Cost — Pairwise Binary (Real LLM Evals, Single-Run Subsample)

Figure: Coverage vs. computation cost for pairwise binary CI methods on real LLM eval data (single-run subsample of multi-run benchmarks), comparing Tango and Bayesian paired across sample sizes.

This pattern could also be an artifact of how the real-data simulation is constructed. Because we subsample from large benchmarks — and large benchmarks contain correlated items and difficulty clusters — a random draw of N = 20 items can easily land in a sparse or lopsided regime that would be unusual in a carefully curated, representative eval set. In a smaller, purpose-built evaluation that attempts reasonable coverage of difficulty, the dominance and sparsity problems would largely disappear.

Finally, the situations where Tango under-performs are also extreme cases where CIs don't tell us much —if one model completely dominates in performance (0% vs 80%), merely eyeballing mean differences may be enough. Similarly, if both models perform very poorly (2% vs 2%), statistics are similarly not that informative. So if a CI method is "slightly off", it may not be a big deal in practice.

That said, evalstats cannot base recommendations on "it's usually fine." Since the failure mode exists and is observable on real benchmark data at small N, our recommendation changes:

• N < 50: use Bayesian paired with 10,000 resamples for binary pairwise CIs. It provides better small-sample coverage in the real-data regime, including the sparse/dominated cases where Tango struggles. We sacrifice speed considerably, but it's safer for extreme cases that may very well appear in practice on poorly constructed evaluation sets or "David vs. Goliath" scenarios.
• N ≥ 50: switch to Tango score interval. At this point Tango reaches its target coverage, is dramatically faster, produces tighter CIs than Bayesian paired, and sits closer to the nominal α level rather than over-covering. Critically, it also does not deteriorate at large N the way the Bayesian paired method does.

This threshold-based switch is what evalstats implements under method="auto" for binary pairwise comparisons.

Multi-Run Coverage (Hierarchical Data)

Increasingly, it has become good practice to collect multiple "runs" per input, to account for the inherent stochasticity in LLM outputs. These multi-run scenarios introduce a second level of variance — run stochasticity within each input — that single-run methods ignore or may not even be able to handle. We ran a separate Monte Carlo simulation (200 MC reps; 10% of the single-run study) over the same distribution suite, sweeping both N (number of inputs) and f_run (runs per input), to compare "flat" methods that just use the mean of runs or take only the first run, against nested resampling methods. For a harder challenge that better matches real LLM eval data (where each "item" is really of different difficulty), run-level noise in these simulations is heteroscedastic: run noise variance differs across inputs.

The key finding is that taking the per-input mean across runs and running standard methods over those means performs reasonably well — Wilson, Tango, and NIG all achieve near-target coverage this way. Input-level clustering is automatically absorbed when the mean is treated as a single observation, so these methods are not as badly hurt as truly flat approaches that treat every run as an independent data point. However, tighter CI widths can be achieved with methods that explicitly account for within-item run variance: by decomposing total variance into a between-item term and a within-item Monte Carlo noise term (as in tango_multirun_mmnt and the nested resampling variants), the interval can be narrowed without sacrificing coverage. The practical upshot is that running on means is a safe baseline, but if you have enough runs per input the multi-run-aware methods will give you meaningfully narrower intervals at no extra compute cost.

CI Coverage vs. Computation Cost — Point Estimates

Figure: Coverage vs. computation cost for point-estimate (mean) CI methods in the multi-run setting, averaged across eval types, run-noise fractions, and sample sizes.

Coverage target is 0.95. ▼ marks severe under-coverage. ★ marks preferred methods. 200 MC reps per condition. "Flat" methods treat multiple runs as independent observations and ignore run-level variance entirely. Nested methods resample inputs, then runs within inputs.

Multi-Run CI Method Performance for Pairwise Comparisons

CI Coverage vs. Computation Cost — Pairwise Comparisons

Figure: Coverage vs. computation cost for pairwise CI methods in the multi-run setting, averaged across eval types, run-noise fractions, and sample sizes.

Coverage target is 0.95. ▼ marks under-target coverage. ★ marks preferred methods. 200 MC reps per condition. Nested methods resample inputs then runs within inputs for each resample step, greatly increasing compute cost but properly accounting for run-level variance.

Multi-Run Coverage with Real LLM Evals Data

Since we had run the other cases on real LLM eval data, we wanted to check if the multi-run results held up there as well. Unfortunately, we could not use OpenEvals (or any other public dataset) because they do not have multi-run data. (Only DOVE does, but they sample only 100 items per benchmark, which is too small for a coverage simulation.) Thus, we had to run benchmarks ourselves on a sample of models. For this test, we only ran benchmarks that produce binary scores, for two reasons: the Tango under-coverage problem was more pronounced there and we wanted to see if it held up on real data; and it is harder to find numeric score benchmarks.

We ran the following benchmarks and models using the Inspect AI benchmark harness and the OpenRouter API. This range of models was chosen as a good sample of model diversity in terms of provider and size:

• Benchmarks: MMLU, ARC, HellaSwag, GSM8K, WinoGrande, TruthfulQA, BoolQ, PIQA, BBQ
• Models: meta-llama/llama-3.1-8b-instruct, mistralai/ministral-8b-2512, google/gemma-3n-e4b-it, ibm-granite/granite-4.1-8b, qwen/qwen3.5-35b-a3b, openai/gpt-4o-mini

As you can see, the above results largely hold up on real data for the binary case. The one difference is that many methods are brought down a little, which is why we recommend NIG-nested here instead of NIG on means.

CI Coverage vs. Computation Cost — Point Estimates (Real LLM Evals)

Figure: Coverage vs. computation cost for point-estimate CI methods in the multi-run setting, measured on real LLM eval data.

CI Coverage vs. Computation Cost — Pairwise Comparisons (Real LLM Evals)

Figure: Coverage vs. computation cost for pairwise CI methods in the multi-run setting, measured on real LLM eval data.

Comparison with Bayesian Methods for Binary Scores

To address the ICML 2025 position paper arguing that practitioners should avoid CLT-style / frequentist intervals for LLM evals in favor of Bayesian beta-posterior methods ("Don't Use the CLT in LLM Evals With Fewer Than a Few Hundred Datapoints"), we ensured the simulation included the exact Bayesian methods the authors advocate, implemented using their bayes_evals library. Our simulation results put some of their recommendations into question. For the binary scores 0/1, we find:

• Wilson score interval lands on target coverage (0.951 synthetic; 0.955 real data) at much cheaper computational cost, so the paper's preferred Bayesian alternative is not needed here.
• For the pairwise comparisons, the story is nuanced. In the synthetic simulation, Tango (0.955) matches Bayesian paired (0.958) in coverage while being nearly 100× faster. But on real LLM eval data — where dominated pairs and jointly-sparse model comparisons are common — Tango under-covers at N < 50, caused by extreme situations of sparse data regimes. Bayesian paired remains better in these small-sample setups. At N ≥ 50, Tango reaches target coverage, is dramatically faster, produces tighter intervals, and does not deteriorate at large N.
• In the new multi-run simulation, the NIG-nested variant achieves the highest point-estimate coverage overall (0.966), well above the authors' Bayesian independent samples method (0.920). NIG-nested is itself a Bayesian method — placing a Normal-Inverse-Gamma prior on (μ, σ²) and returning a marginal Student-t credible interval — but uses a Gaussian likelihood rather than a Binomial/Beta model. The result shows that the authors' specific Beta-posterior model is not uniquely suited to the hierarchical binary case; a different conjugate Gaussian model substantially outperforms it.

In the course of this analysis, we also uncovered a severe problem in the authors' library that makes it unsuitable for practical usage on benchmark data: around N≈200+, the bayes_evals Bayesian paired method begins to fail catastrophically with overconfident CIs that continue to decline as N increases. (We ran additional simulations up to N=1000 confirming this pattern. The reason why seems to be that the importance sampling procedure becomes unstable at larger N.) Because of this failure mode, Bayesian paired is not recommended for N ≥ 100 — and it should never be used beyond N ≈ 200. For small samples (N < 50), it retains a narrow use-case: real-data simulations show it handles dominated and jointly-sparse pairs better than Tango in that regime. For larger samples, and as the default once N reaches 50, practitioners could switch to the Tango and Wilson intervals (and their multi-run analogues), which are fast, well-calibrated, and stable. For multi-run data, which was not considered in the paper, practitioners should use NIG-nested for point estimates and Tango's multi-run variant for pairwise comparisons.

Is BCa bootstrap a good method for estimating CIs?

If you chat with an LLM, it might start to go off about how BCa is the gold standard. Don't listen to it. BCa is actually one of the worst methods in most of our simulations. We aren't the only skeptics, either. One simulation study found:

This is consistent with our own simulation results on both synthetic and real LLM eval data.

Which p-value Method Works Best?

We also ran another Monte Carlo simulation study to compare p-value methods for pairwise testing across LLM eval data types and sample sizes, implemented in simulations/sim_compare_pvalues.py. The criterion here is Type I error control under the null (target: α = 0.05) plus consistency across data distributions and N.

Wilcoxon Signed Ranks is our recommended default for pairwise comparisons for those adopting a NHST framework. In this simulation, it is slightly conservative, consistently performant across sample sizes, data types, and distributions, and does not show pathological behavior where Type I error is inflated above the nominal threshold (α = 0.05). Wilcoxon is a non-parametric test that makes fewer assumptions about the data distribution, which likely contributes to its robustness in this context. It is also widely implemented, well-understood, and computationally efficient, making it a practical choice for developers.

That said, the Sign Test is also a reasonable choice for safety-critical settings where false positives are especially costly, because it is much more conservative than Wilcoxon.

For comparing 3+ models/prompts, we recommend the Friedman test as the default, with Wilcoxon signed ranks as the post-hoc pairwise test. Friedman is a non-parametric test for comparing more than two related samples, it doesn't assume normality, and it's robust to homoscedasticity (variance of residuals). Friedman tests are already used in some ML contexts for comparing models, as advocated by Demsar (2006). Note that evalstats does not output Friedman by default, since it does not foreground a NHST-style analysis pipeline, but it is available as an option for users who want to use it.

For multiple comparisons correction, we recommend the Benjamini-Hochberg procedure (correction="fdr_bh") to those in lower-stakes contexts, as it which performed comparably to other methods in our simulation, but had the greatest power among methods that maintained Type I error control. In extremely safety-critical contexts, the Bonferroni-Holm correction (correction="holm") is more appropriate for its strong control over false positives, albeit at the cost of increased false negatives.