Characterizing LLM performance via a Bayesian lens

1 Introduction

Large language models (LLMs) represent a significant leap in artificial intelligence, demonstrating remarkable capabilities in diverse applications, from specialized domains such as medicine [1] to broader information processing tasks [2]. Despite fast-evolving LLM capabilities, their associated evaluation methodology has not kept up with the need for nuanced performance analysis. Current standard practices often rely on aggregate performance metrics, such as point estimates of accuracy on benchmark datasets [3, 4]. While such evaluation metrics are easy to interpret, they obscure critical aspects of model behavior by failing to adequately account for either the inherent stochasticity in LLM response generation or the substantial heterogeneity in difficulty levels across different questions or tasks within benchmarks.

Furthermore, recent years have witnessed the meteoric rise of agentic AI which has become one of the most pivotal applications of LLMs. In agentic AI, an LLM is tasked with choosing a tool call from multiple tools available, essentially amounting to a multiple choice question. Because the results of the tool call are used as input in the next LLM call, variation in tool selection cascades across the iterative agent loop. In many cases, there is not one single “correct” tool call but there could be multiple options that are equally reasonable depending on the situation. Given this context, the ability to accurately and holistically characterize the behavior of LLMs not only in terms of accuracy but also capturing their intra-question stochasticity and inter-question heterogeneity becomes even more critical.

This urgency for more holistic evaluation approaches of LLMs is even more pressing in high stakes domains such as healthcare and real time control systems [5, 6]. Uncertainty Quantification (UQ) can also be important within model deployment, e.g., within the context of active learning or when adding a rejection option to a classification task when uncertainty is too high. Although the probabilistic underpinnings of LLMs are well-known [7], and response variations can arise even under ostensibly deterministic settings due to computational and implementation factors [8, 9], contemporary evaluation frameworks often lack the statistical sophistication to properly disentangle and quantify these effects [10].

There have been attempts to remedy this limitation. Beyond using simple traditional frequentist approaches, e.g., reporting standard deviation and/or estimated confidence intervals along with the chosen performance metric [10], prior work on UQ has often relied on Bayesian or ensemble-based techniques. This includes Monte Carlo dropout, which performs stochastic forward passes to assess output variability [11] and deep ensembles where multiple models are trained independently [12]. Recent studies have also categorized sources of variance in terms of aleotoric (irreducible; due to the intrinsic noise in the data itself) vs. epistemic (reducible; due to ignorance) uncertainty [13] or in terms of uncertainty dimensions specific to LLMs like in the work of Liu et al. [14] where uncertainty was quantified in terms of input, reasoning, parameter, and prediction. Atf et al. [15] also presented a domain-specific review of UQ within the context of medicine where different UQ approaches were described and categorized but not quantified. Liang et al. [4] presented a thorough framework for holistic evaluation of LLMs ranging from accuracy to bias and toxicity, but uncertainty was not the focus.

While several studies have used Bayesian techniques to quantify uncertainty in deep learning, they have often been centered on modeling the weights themselves as probability distributions (Bayesian Neural Networks; BNN) [16–18]. The computational cost of exact Bayesian inference over model weights becomes intractable at the scale of modern LLMs [19], which typically employ billions and in some cases trillions of parameters [20]. Even though techniques have been suggested to deal with this problem [11, 19], the approach of modeling weights as distributions is not possible in the case of closed-source models, which is relevant since many of the most commonly used state-of-the-art models are at present closed-source.

In the approach presented here, we avoid this problem by modeling the binary correct/incorrect output of the LLMs as a stochastic variable sampled from a binomial distribution, the parameter of which (likelihood of answering correctly) is in turn modeled as a stochastic variable sampled from a Beta distribution with two hyperparameters [21, 22]. This hierarchical BetaBinomial Bayesian modeling framework gives us the means to quantify uncertainty in a computationally inexpensive way that is model-agnostic and can be applied to closed-source models as well. The core advantage of this approach is its ability to simultaneously estimate distinct facets of performance, namely average accuracy, response consistency (stochasticity), and inter-question heterogeneity, while providing a principled quantification of uncertainty through full posterior distributions for these interpretable characteristics [11, 23, 24].

This paper makes the following contributions: We present the mathematical formulation and implementation of this hierarchical Bayesian framework. We demonstrate its application by evaluating a diverse set of contemporary LLMs on standard benchmarks, including under semantic perturbations designed to probe robustness. Our analysis, detailed in subsequent sections, reveals nuanced performance profiles, quantifies model-specific sensitivities to prompt variations, and uncovers distinct behavioral patterns linked to model architecture and task domain. By offering a richer, statistically grounded perspective on LLM evaluation, this work provides an essential tool for moving beyond simplistic leaderboard rankings toward the development and deployment of more dependable and thoroughly understood artificial intelligence systems. The end-to-end evaluation pipeline is schematically illustrated in Figure 1.

2 Method

Evaluating LLMs robustly requires disentangling multiple sources of variation [4, 6]. We introduce a hierarchical Bayesian framework [23, 24] tailored for the rigorous assessment and comparison of LLM performance across diverse questions, especially when examining behavior under distinct conditions (e.g., original vs. rephrased prompts). Our Beta-Binomial hierarchical model explicitly addresses three fundamental aspects of LLM evaluation:

• Intra-Question Response Variability: The inherent stochasticity in LLM response generation can lead to varying outcomes even when processing the same question repeatedly. This trial-level variability is captured via the model's Binomial likelihood component. Its expected magnitude, averaged across questions, serves to quantify the LLM's response consistency under each experimental condition.
• Inter-Question Heterogeneity: Questions naturally differ in their intrinsic difficulty or other properties affecting LLM success. We model this by treating the underlying success probability for each question as a draw from a condition-specific Beta distribution. Characterized by its mean and concentration, this distribution quantifies the heterogeneity across the question set, providing insights beyond simple aggregate accuracy.
• Systematic Condition Effects: A key capability of this framework is the explicit modeling and comparison of LLM performance across different conditions. Assessing model sensitivity to prompt variations, such as rephrasing, is vital for understanding robustness. Crucially, the Bayesian approach allows for the precise estimation of the difference between conditions in key performance metrics.

We analyze paired experimental data, the same questions under two conditions, while modeling each condition with its own hierarchical structure. This structure leverages partial pooling across questions, yielding regularized and stable parameter estimates for both individual questions and population characteristics [22]. Bayesian inference provides full posterior distributions for comprehensive uncertainty quantification. Critically, we estimate the posterior distribution of the difference between condition parameters, enabling rigorous assessment of intervention effects, such as question rephrasing impact on LLM performance [3].

LLMs were used exclusively for generating rephrased benchmark questions under controlled prompts. No LLMs were used to generate the statistical methodology, analysis code or results.

2.1 Data and notation

Consider a set of N unique questions, indexed by i = 1, …, N. Each question i is presented under two conditions, m = 1 (original) and m = 2 (rephrased). For question i under condition m, let yij(m)∈{0,1} denote the correctness (1 for correct, 0 for incorrect) of the response on the j-th repeated trial, where j=1,…,Ji(m). The total number of trials for question i in condition m is Ji(m). The total number of correct responses is given by Equation 1:

yi(m)=∑j=1Ji(m)yij(m), m∈{1,2}.(1)

The observed data consist of the set of pairs of counts {(Ji(1),yi(1)),(Ji(2),yi(2))} for i = 1, …, N.

2.2 Likelihood: Beta-Binomial construction

For a specific question i under condition m, let θi(m)∈[0,1] represent the underlying probability that the LLM provides a correct response to a single trial. Assuming conditional independence of trials given θi(m), the number of correct responses yi(m) follows a Binomial distribution (Equation 2):

yi(m)∣θi(m),Ji(m)~Binomial(Ji(m),θi(m)).(2)

This Binomial likelihood models the stochasticity inherent in the LLM's response generation process for a fixed question i and condition m over Ji(m) independent trials. To model question heterogeneity (i.e., variability in θi(m) across questions i within condition m), we assume these probabilities are drawn from a common distribution specific to the condition. We place a Beta prior on each θi(m) (Equation 3):

θi(m)∣αm,βm~Beta(αm,βm),(3)

where the hyperparameters (α_*m* > 0, β_*m* > 0) are shared across all questions i within condition m. Integrating out the latent probability θi(m) yields the marginal likelihood for the observed count yi(m), which follows a Beta-Binomial distribution [21], often employed for modeling Binomial data exhibiting overdispersion relative to a simple Binomial model with a fixed probability:

yi(m)∣Ji(m),αm,βm~BetaBinomial(Ji(m);αm,βm), m∈{1,2}.(4)

The probability mass function (PMF) for this distribution is

P(yi(m)∣Ji(m),αm,βm)=(Ji(m)yi(m))×B(αm+yi(m),βm+Ji(m)-yi(m))B(αm,βm)(5)

where B(·, ·) denotes the Beta function. This hierarchical formulation induces partial pooling: the hyperparameters α_*m*, β_*m* capture the central tendency and dispersion of the correctness probabilities across all questions within condition m, while allowing individual θi(m) values to vary [22].

This hierarchical structure connects our framework to Bayesian Item Response Theory (IRT) and random-item-effects psychometric models, which similarly treat item difficulty as a latent quantity drawn from a population distribution [25, 26]. The concentration parameter κ_*m* plays a role analogous to the spread of item-difficulty parameters in IRT. However, a fundamental difference is that IRT is designed to separate person-ability from item properties across a population of distinct examinees, whereas our setting involves a single fixed model queried repeatedly on each item. The dominant source of variation in our data, trial-to-trial response stochasticity, thus has no counterpart in standard IRT and no identifiable person-ability dimension. We model it explicitly through the Binomial likelihood and summarize it via γ̄. Our mean-concentration Beta-Binomial parameterization is consequently chosen to yield three directly interpretable and identifiable quantities (μ, κ, γ̄) under this repeated-trial design, while remaining computationally inexpensive and model-agnostic for closed-source LLMs.

2.3 Mean–concentration hyperpriors

For enhanced interpretability and prior specification, we parameterize the Beta distribution's parameters (α_*m*, β_*m*) for each condition m∈{1, 2} in terms of its mean μ_*m*∈[0, 1] and concentration κ_*m* > 0:

αm=μmκm and βm=(1-μm)κm.(6)

This mean-concentration parameterization [23, Chapter 5] is equivalent to the standard (α_*m*, β_*m*) parameterization, with the inverse relationship given by Equation 7:

μm=αmαm+βm, κm=αm+βm.(7)

For each condition m:

• μ_*m* represents the expected accuracy across the population of questions in that condition, E[θi(m)]=μm.
• κ_*m* functions as a concentration parameter governing the dispersion of the question-level probabilities θi(m) around the mean μ_*m*, with Var[θi(m)]=μm(1-μm)κm+1. Because this variance depends on both κ_*m* and μ_*m*, κ_*m* maps monotonically to inter-question heterogeneity only at a fixed μ_*m*. A larger κ_*m* implies lower dispersion and a smaller κ_*m* implies higher dispersion, but comparisons of heterogeneity across conditions with differing μ should reference the Beta variance directly rather than κ alone.

Note that κ_*m* is a population-level concentration parameter governing the dispersion of question difficulties, and is distinct from a Beta-Binomial smoothing count added to the per-question trials. Its effect is therefore on the spread of the θi(m) across questions rather than on the effective number of trials. For example, at μ = 0.86, κ = 0.03 corresponds to SD[θ_*i*]≈0.34 and a near-bimodal difficulty profile, whereas κ = 0.5 corresponds to SD[θ_*i*]≈0.28 and a substantially more homogeneous one.

We specify hyperprior distributions for μ_*m* and κ_*m*. Assuming prior independence between conditions and between parameters within a condition, we use standard weakly informative priors:

• For the mean accuracy: μ_*m*~Beta(1, 1), which is equivalent to Uniform(0, 1).
• For the concentration: We place a prior on the logarithm, logκm~Normal(μlogκ=0,σlogκ2=1.52) (i.e., variance 2.25, standard deviation 1.5). This choice is common practice [24] and avoids restricting κ_*m* unduly while discouraging extremely small or large values a priori.

Thus, the joint hyperprior is

p(μ1,κ1,μ2,κ2)=p(μ1)p(κ1)p(μ2)p(κ2).

2.4 Posterior distribution

Applying Bayes' theorem, the joint posterior distribution of the hyperparameters, given the observed data D={(Ji(1),yi(1)), (Ji(2),yi(2))}i=1N, is proportional to the product of the likelihoods and the priors:

p(μ1,κ1,μ2,κ2∣D)∝ ∏m=12∏i=1NP(yi(m)∣Ji(m),αm(μm,κm),βm(μm,κm)) ×p(μ1)p(κ1)p(μ2)p(κ2)(8)

where P(yi(m)∣·) is the Beta-Binomial PMF given in Equation 5, with α_*m* and β_*m* substituted using Equation 6, and p(μ_*m*), p(κ_*m*) are the specified hyperprior densities.

2.5 Bayesian inference via MCMC

Since the posterior distribution (Equation 8) is typically intractable analytically, we employ Markov Chain Monte Carlo (MCMC) methods for inference [27] using PyMC [28]. We denote the index for MCMC samples (draws) by k.

• Initialization: Four chains are initialized. While Maximum A Posteriori (MAP) estimates can inform initialization, random initialization within the parameter space or based on data summaries is often preferred to better explore potential multimodality.
• Sampling: Samples (μm(k),κm(k)) are drawn from the posterior p(μ₁, κ₁, μ₂, κ₂∣D) using an efficient gradient-based sampler, such as the No-U-Turn Sampler (NUTS) [29], available in probabilistic programming frameworks like PyMC [28] or Stan [30]. A sufficient number of warm-up (50,000) iterations are discarded.
• Post-Warm-up Sampling: A large number of posterior samples (200,000) are collected after the warm-up across multiple chains to ensure stable estimates of posterior summaries.
• Convergence Diagnostics: Convergence is assessed rigorously. Checks include monitoring for sampler divergences, ensuring the Gelman–Rubin diagnostic R^ is close to 1.0 (e.g., R^<1.01) for all parameters [31, 32], and verifying sufficiently large effective sample sizes (ESS) for all parameters of interest [23, Chapter 11]. The convergence diagnostic across all models and datasets is presented in our posterior summaries (see Supplementary Tables S1–S4).

We emphasize that MAP estimation is used only to initialize the sampler and is not used for inference. All reported quantities are posterior summaries derived from the post-warm-up MCMC samples.

2.6 Posterior summaries and derived measures

The MCMC procedure yields samples {(μm(k),κm(k))∣m=1,2}k=1K from the joint posterior distribution p(μ₁, κ₁, μ₂, κ₂∣D), where K is the total number of post-warm-up samples. From these samples, we can derive posterior distributions (represented by the collection of samples) for several key quantities. When discussing these quantities as random variables under the posterior, we omit the index k; when referring to values derived from specific MCMC samples, we include the index k.

• Population Mean Accuracy *μ^(m)*: Represents the expected correctness probability across the population of questions for condition m. Higher values signify better average performance. The MCMC samples μm(k) constitute its posterior distribution p(μ^(*m*)∣D), quantifying the uncertainty associated with this population average.
• Population Concentration *κ^(m)*: The concentration parameter for condition m, inversely related to the inter-question heterogeneity (variance in difficulty). Increasing values of *κ^(m) correspond to decreasing heterogeneity, meaning individual question probabilities θi(m) are increasingly concentrated around the population mean μ^(m). Conversely, as κ^(m*) → 0, heterogeneity increases, indicating greater dispersion of question-level probabilities θi(m). The samples κm(k) form its posterior distribution p(κ^(*m*)∣D), reflecting uncertainty about the degree of concentration.
• Question-Level Correctness Probability θi(m): The posterior distribution for the underlying correctness probability θi(m) for a specific question i under condition m. Given the MCMC samples of the hyperparameters (μm(k),κm(k)), which determine αm(k)=μm(k)κm(k) and βm(k)=(1-μm(k))κm(k) at iteration k, the conditional posterior for θi(m) follows from the standard Beta-Binomial conjugate update (Equation 9): Drawing samples θi(m)(k) from this distribution for each MCMC iteration k allows the construction of the empirical marginal posterior p(θi(m)∣D) for each item i and condition m. Comparing p(θi(1)∣D) and p(θi(2)∣D) enables the identification of specific questions most affected by the condition change.
• Expected Per-Trial Intra-Question Variance γi(m): The posterior expectation of the variance associated with a single trial for question i under condition m, denoted Var(yij(m)∣θi(m))=θi(m)(1-θi(m)). This variance arises from the inherent stochasticity in the LLM's response generation for a given question. Formally (Equation 10): Samples vi(m)(k)=θi(m)(k)(1-θi(m)(k)), derived from each posterior sample θi(m)(k), represent the marginal posterior distribution p(θi(m)(1-θi(m))∣D). The mean of this distribution estimates γi(m), reflecting the expected single-trial response inconsistency for question i.
• Mean Expected Total Intra-Question Variance γ̄(m): The average, across all N questions in condition m, of the expected total variance (Ji(m)γi(m)) of the success count yi(m). This metric quantifies the overall magnitude of response stochasticity within questions for condition m. It is computed as shown in Equation 11:

A value of γ̄(m) close to zero indicates high average response consistency (low stochasticity) across the question set, whereas larger values indicate substantial average stochasticity.

Since the per-trial variance θ_*i*(1−θ_*i*) is maximized at θ_*i* = 0.5, γ̄(m) carries an expected dependence on the mean accuracy μ_*m*, with intermediate accuracies contributing larger intra-question variance. This dependence is not deterministic in μ_*m*, however, because γ̄(m) is a functional of the entire posterior distribution of the θi(m) and therefore also reflects their dispersion, which is controlled by κ_*m*. Two conditions with the same μ but differing θ dispersion will in general have different γ̄. We consequently interpret γ̄ jointly with μ and κ.

2.7 Posterior predictive checks (PPC)

We assessed model adequacy using posterior predictive checks (PPC) [23, 33]. Replicated datasets {y~i(m)(k)} were simulated from the posterior predictive distribution by drawing from the likelihood using posterior samples of (αm(k),βm(k)). Summary statistics calculated from these replicated datasets were compared to those from the observed data yi(m) to diagnose potential model misspecification. Across all models and datasets, inspection of the posterior predictive checks confirmed model fidelity, with the observed data consistently falling within the range of high probability density of the simulated replicated datasets. Illustrative examples of these PPC plots for each model-benchmark pairing can be found in Supplementary material (e.g., Supplementary Figure S7, and subsequent related figures in Sections 2.2–2.5).

In summary, our hierarchical Bayesian framework employs a Beta-Binomial likelihood (Equation 4), parameterized via interpretable mean and concentration (Equation 6), to rigorously model both LLM response stochasticity and inter-question heterogeneity [23, 24]. Model adequacy is assessed through posterior predictive checks. This integrated methodology delivers the holistic probabilistic performance profiles essential for robust LLM evaluation beyond simplistic point estimates.

3 Experimental setup

3.1 Large language models

We evaluated a diverse set of LLM models:

• Claude-Sonnet-3.5 [34]
• Gemini-Flash-2.0 [35]
• GPT-4o [36]
• O1 [37]
• O1-mini [38]
• O3-mini [39]
• DeepSeek-V3 [40]
• DeepSeek-V3-0324 [41]
• DeepSeek-R1 [42]
• QwQ-32B [43]
• Llama-3.3-70B-Instruct [44]

All models were accessed through the OpenRouter stateless API using versioning to ensure full reproducibility. No additional conversation context was provided between each API call. Model configurations and system prompts were pinned and are publicly accessible at https://github.com/pfizer-opensource/llm-uncertainty. The “reasoning vs. non-reasoning” grouping is based on OpenRouter model definitions.

3.2 Datasets

We curated a set of 100 multiple-choice questions sampled from four widely-used benchmarks:

• MMLU [3]
• MedMCQA [45]
• CommonsenseQA (CSQA) [46]
• MathQA [47]

This collection spans factual knowledge, medical reasoning, commonsense inference, and mathematical problem-solving, allowing us to analyze heterogeneous question difficulty and domain-specific variability.

Each question was evaluated in two forms:

• Exact version—the original benchmark phrasing.
• Rephrased version—a semantically equivalent variant generated using controlled LLM transformations.

Rephrased questions were manually reviewed to ensure semantic equivalence. To capture response stochasticity, each LLM was queried for 100 independent trials per question per condition (exact and rephrased), resulting in 20,000 model outputs per benchmark per model. All models were evaluated on an identical set of question pairs for strict comparability.

3.3 Prompt synthesis

LLMs were used to generate semantically equivalent rephrased versions of benchmark questions under controlled prompts. No LLMs were used to generate or modify the statistical methodology, analysis code, or conclusions

Rephrased versions were generated with a fixed, controlled instruction rather than at the model's discretion, using Claude 3.5 Haiku. The instruction directed the model to preserve the exact meaning and all specific content, including numbers, percentages, and other quantitative or factual details, while varying only the surface phrasing. The exact prompt used was:

“Rephrase the following question in a different way while preserving its exact meaning. Keep the rephrased version clear and concise. Only return the rephrased question, nothing else. Make sure to not miss any important details or context especially numbers, percentages, and other specific information mathematical or non-mathematical.”

The objective was to perturb linguistic form while holding semantic content constant, so that any measured change in performance reflects sensitivity to phrasing rather than to a change in what the question asks. Every rephrased item was manually reviewed to confirm semantic equivalence to its original, and items that altered meaning were regenerated. We deliberately did not condition the rephrasing on psychometric or survey methodology constructs, since the goal was to isolate robustness to naturalistic linguistic variation rather than to engineer difficulty.

Each evaluation followed a unified prompt format. Five in-context examples were provided to illustrate the answer structure and improve comprehension. These examples:

• were drawn from the same benchmark as the target question,
• were fixed for all models and all trials,
• were identical for exact and rephrased variants.

This ensures that performance differences arise from question semantics–not from prompt variability.

3.4 Evaluation parameters

All models were queried with:

• temperature = 0,
• random seed = 42,
• OpenRouter default max tokens, top-k, and top-p values.

We note that T = 0 corresponds to greedy decoding, which is deterministic in principle but not in practice. Residual response variability arises from sources independent of the sampling temperature, including non-associative floating-point arithmetic and hardware-level nondeterminism in GPU execution, load-dependent batching and routing at the serving layer, and, for mixture-of-experts architectures, batch-dependent expert selection [8, 9]. Because models were queried through a stateless hosted API, these effects are fully present in our data. Fixing T = 0 serves two purposes. It provides the most conservative, maximally deterministic setting, so that any observed stochasticity reflects intrinsic model and serving behavior rather than deliberate sampling noise, and it holds decoding fixed so that differences in variability across conditions can be attributed to semantic rephrasing rather than to sampling temperature. A single common temperature is further required for strict cross-model comparability. The substantial and condition-dependent intra-question variability we observe despite this deterministic setting confirms that T = 0 does not suppress the phenomenon of interest but rather bounds it from below.

3.5 Computational resources

All evaluations were conducted on an AWS EC2 g5.2xlarge instance (8 vCPUs, 32 GB RAM, 1 × NVIDIA A10G GPU). Computational cost was modest. The evaluation for all models had a runtime of approximately 5 h. Bayesian inference and statistical analyses were executed using PyMC and ArviZ on the same machine.

4 Results

Applying our framework to LLM-generated answers to the MMLU question dataset, we characterized LLM performance across three key dimensions: average accuracy (μ), inter-question homogeneity (κ), and intra-question response variability (γ̄), revealing distinct effects of question rephrasing. Posterior summaries for the examined LLMs are presented in Tables 1, 2.

Table 1. Hierarchical Bayesian posterior summaries for LLM performance on the MMLU benchmark under the Exact question condition.

Reported parameters include mean accuracy (μ), inter-question homogeneity (κ), and expected intra-question variance (γ̄), with posterior means and 95% Highest Density Intervals (HDIs). Reasoning models appear in the top block, followed by non-reasoning models.

Table 2. Hierarchical Bayesian posterior summaries for LLM performance on the MMLU benchmark under the Rephrased question condition.

Reported parameters include mean accuracy (μ), inter-question homogeneity (κ), and expected intra-question variance (γ̄), summarized as posterior means with 95% Highest Density Intervals (HDIs). Reasoning models are listed first, followed by non-reasoning models.

Throughout, all population parameters (μ, κ, γ̄) and their between-condition differences are reported as posterior means with 95% Highest Density Intervals rather than as point estimates. Because we retain the full posterior, we additionally report the posterior distribution of the difference between conditions for each parameter, enabling direct probabilistic assessment of the effect of rephrasing. Credibility of an effect is judged by whether the 95% HDI of the corresponding difference excludes zero, or equivalently by separation of the per-condition HDIs.

We note that κ and γ̄ are mathematically linked to μ and are therefore interpreted conditionally rather than as independent robustness dimensions. Posterior contrasts indicate that their observed divergences across conditions are not trivially explained by accuracy shifts alone. We accordingly interpret κ as a concentration parameter and quantify condition-level differences in question-difficulty dispersion using the μ-adjusted Beta variance rather than κ in isolation. Fixed in-context exemplars were reused across conditions to isolate the effect of semantic rephrasing.

Question rephrasing consistently reduced mean accuracy (μ) across all evaluated models Tables 1, 2). This performance decline was statistically credible, evidenced by distinct downward shifts in the 95% Highest Density Intervals (HDIs) between conditions. The extent of this reduction varied by architecture: the reasoning model QwQ 32B exhibited a modest drop [Δμ≈0.06; decreasing from 0.85 [0.79, 0.90] to 0.79 [0.73, 0.85]], whereas the non-reasoning Gemini 2.0 Flash suffered a more pronounced decline [Δμ≈0.11; decreasing from 0.86 [0.79, 0.92] to 0.75 [0.69, 0.81]; Figure 2a]. Comprehensive posterior summaries and standard descriptive statistics for all datasets are provided in Supplementary Tables S1–S8, respectively.

The effect of rephrasing on inter-question homogeneity κ differed notably between model groups. Several non-reasoning models showed a significant rise in κ when rephrasing the questions, such as Gemini 2.0 Flash where κ increased from 0.03 [0.01, 0.06] to 0.50 [0.34, 0.66] (Figure 2a), indicating that rephrasing substantially homogenized question difficulty for them. In contrast, reasoning models such as O3 mini and Deepseek R1 exhibited no credible change in homogeneity, with κ changing negligibly from 0.34 [0.19, 0.50] to 0.35 [0.23, 0.48] in the case of O3 mini (Figure 2c) and from 0.36 [0.18, 0.55] to 0.37 [0.23, 0.52] in the case of Deepseek R1, suggesting their perception of relative question difficulty was more robust to the semantic perturbation of rephrasing as compared to the non-reasoning models (Tables 1, 2).

Across all evaluations, rephrasing substantially increased intra-question variability (γ̄), indicative of heightened response stochasticity. This shift was marked by clearly separated 95% HDIs between conditions for every model (Tables 1, 2). Notably, the magnitude of this increase was consistently lower for reasoning models compared to their non-reasoning counterparts. To illustrate this architectural contrast, the reasoning model Deepseek R1 saw γ̄ rise by approximately 1.9 [from 2.37 [2.21, 2.55] to 4.25 [4.04, 4.46]], whereas Gemini 2.0 Flash experienced a significantly sharper surge of approximately 5.1 [from 0.49 [0.46, 0.52] to 5.58 [5.33, 5.83]; Figure 2a].

To confirm that these γ̄ increases are not a mechanical consequence of the concurrent decrease in μ, we compared the observed posterior shift in γ̄ against the change predicted by the μ shift alone. The residual increase in γ̄ beyond this expectation was credible across models, indicating that rephrasing amplifies intra-question stochasticity over and above the effect attributable to reduced accuracy.

These results suggest that while rephrasing invariably drives up response inconsistency, reasoning models are relatively more resilient to this effect. Crucially, this distinction was obscured when relying on standard descriptive statistics (Supplementary Tables S5–S8). While raw standard deviations generally increased with rephrasing, they failed to capture the consistent, architecture-dependent stratification revealed by the Bayesian parameters.

Figure 2 illustrates the impacts of question rephrasing by contrasting population parameter distributions (Figures 2a, c) and corresponding item-level accuracy (θ_*i*, Figures 2b, d) for Gemini 2.0 Flash (non-reasoning) and O3 mini (reasoning). While both models exhibit the common trends of decreased mean accuracy (μ) and substantially increased intra-question variability (γ̄) upon rephrasing, their change in inter-question homogeneity (κ) diverges significantly. Gemini 2.0 Flash (Figure 2a) shifts from extremely low homogeneity (κ≈0.03) toward moderate homogeneity (κ≈0.5), indicating significant homogenization of question difficulty. This is reflected at the item level (Figure 2b), where the θ_*i* estimates for exact questions (blue) is concentrated around either 0 or 1, consistent with very low κ, while the rephrased estimates (orange) are more dispersed. Conversely, O3 mini (Figure 2c) maintains stable moderate homogeneity (κ≈0.35) across both conditions, corroborated by the similar spread of θ_*i* estimates in Figure 2d for exact (blue) and rephrased (orange) conditions. Despite O3 mini's stable homogeneity profile, both models show increased response stochasticity, evidenced by the rightward shift in the γ̄ distributions (Figures 2a, c) and visually mirrored by the generally wider 95% HDIs for the rephrased θ_*i* estimates (Figures 2b, d, orange bars). This multi-level visualization underscores how population parameters aggregate item-level phenomena and reveals quantitatively different model sensitivities to semantic perturbations, particularly regarding the homogenization of task difficulty. Comprehensive plots of this nature, including Posterior Predictive Checks (PPC), question-level forest plots, and parameter comparison distributions for all evaluated models across all four datasets, are available in Supplementary material, Sections 2.2–2.5 (Supplementary Figures S5–S136).

Analysis of the MATHQA benchmark revealed responses to rephrasing partly mirroring MMLU, yet exhibiting critical divergence concerning inter- and intra-question heterogeneity (Figure 3). Models displayed high baseline accuracy (μ_*exact* > 0.89), which consistently decreased upon rephrasing (median Δμ = −0.047). Concurrently, baseline intra-question stochasticity (γ̄exact, initially spanning [0.50, 3.97]) was universally amplified (median Δγ̄=+1.15). The primary divergence occurred in inter-question homogeneity: whereas baseline homogeneity κ varied widely, spanning [0.09, 1.07], rephrasing induced a bifurcation strongly correlated with model type. Most reasoning models shifted toward decreased homogeneity (median Δκ = −0.102). In striking contrast, non-reasoning models predominantly experienced increased homogeneity (median Δκ = +0.124), perceiving the questions as more homogeneous post-rephrasing (e.g., gpt_4o, claude35). This pronounced domain and model-type-dependency in homogeneity response highlights the complex nature of LLM output robustness. Similar visualizations for other datasets are provided in Supplementary Figures S1–S4, illustrating benchmark-specific nuances. The distributional shifts between Exact and Rephrased conditions across all models, quantified via normalized Bhattacharyya distances, are summarized in Figure 4.

5 Discussion

As the use of LLMs continues to grow explosively across virtually all facets of modern society, their robust evaluation has become more critical than ever. While first-order accuracy metrics provide an easily interpretable measure of LLM performance, they overlook the inherent output variability of these models [3, 4, 6, 48, 49]. Given that contemporary applications increasingly employ agentic workflows, where LLM outputs are chained together and output variability therefore compounds, rigorous characterization of this variability should be considered as important as the average output accuracy itself.

By simultaneously modeling output variability as intra-question stochasticity (γ̄) and inter-question homogeneity (κ) in a hierarchical Bayesian Beta-Binomial framework, we obtain full posterior distributions for mean accuracy μ, as well as κ, and γ̄, enabling a thorough characterization of this variability [5, 23, 24, 50]. Crucially, our output-based modeling is computationally inexpensive and model-agnostic, making it readily applicable to models of any size, including closed-source LLMs.

It is useful to position our framework against existing uncertainty-aware evaluation approaches. Methods such as Monte Carlo dropout [11], Bayesian Neural Networks [19], and deep ensembles [12] quantify uncertainty by accessing or perturbing model internals or by training multiple models. These approaches are not applicable to closed-source, API-only LLMs, which represent much of current state-of-the-art practice, because weights, internal stochastic layers, and retraining are unavailable. Our output-based Beta-Binomial formulation deliberately avoids this requirement and is therefore model-agnostic. Relative to frequentist reporting of standard deviations and confidence intervals, which is applicable in our setting, our results (Section 4) show that such descriptive statistics obscure the architecture-dependent differences in inter-question homogeneity (κ) and intra-question stochasticity (γ̄) that the Bayesian decomposition makes visible. A direct empirical comparison against internals-based methods, restricted to open-weight models where those methods can be applied, is a natural direction for future work.

Applying this framework to publicly available question datasets revealed some critical vulnerabilities: semantic perturbations like question rephrasing consistently reduced mean accuracy (μ) and, importantly, markedly amplified intra-question stochasticity (γ̄) across all tested LLMs and benchmarks (Table 2, Figures 2, 3) [51]. This highlights a performance brittleness concealed by single-trial evaluations and underscores the challenge of response inconsistency [10, 52] (see detailed results for all models and benchmarks in Supplementary material).

Our framework also uncovered model- and domain-specific responses to prompt variations by examining inter-question homogeneity κ. On MMLU, rephrasing led non-reasoning models to perceive questions as more homogeneous (increased κ), while reasoning models showed stable homogeneity (Table 2, Figure 2). Intriguingly, this trend flipped on MATHQA (Figure 3): rephrasing typically decreased homogeneity for reasoning models but increased it for non-reasoning ones, highlighting that model architecture and dataset traits jointly influence robustness [53]. This observation suggests that underlying architectural or algorithmic differences interact profoundly with domain characteristics to shape responses to semantic variation. Consequently, evaluating inter-question homogeneity (κ) provides a crucial diagnostic dimension beyond aggregate accuracy (μ) and average stochasticity (γ̄).

6 Limitations and future work

While insightful, the proposed framework has limitations. Bayesian MCMC inference, though rigorous, is more computationally intensive than point estimation, potentially affecting scalability [27]. Priors, though weakly informative, can influence posterior inference [23]. The model assumes question exchangeability, which may not fully capture benchmark structures (e.g., topic clusters). Our empirical validation, while illustrative, used a specific set of models, benchmarks, and one perturbation type (rephrasing). Moreover, the Beta-Binomial approach is tailored for questions with a binary right or wrong answer, not for evaluating open-ended responses. We evaluated all models under a single fixed decoding temperature to ensure comparability. Characterizing how the accuracy, homogeneity, and stochasticity profiles evolve across a range of temperatures, including model-specific recommended settings for reasoning models, is a promising extension.

Future work can extend this research by applying the framework to more diverse benchmarks (e.g., HellaSwag, BoolQ) [54, 55], and perturbation types (e.g., adversarial, multilingual) [56, 57]. Investigating richer hierarchical structures incorporating question-level features could yield deeper insights. In particular, incorporating IRT-style item covariates (e.g., a discrimination parameter) or a model-ability layer for joint comparison of multiple LLMs would bridge our approach with established psychometric modeling. A systematic empirical benchmarking of our framework against internals-based uncertainty methods, feasible only for open-weight models, would further clarify the trade-offs between model-agnostic and internals-dependent approaches. Exploring lighter inference methods like variational inference [58] could enhance scalability. Adapting this framework for other response types (e.g., categorical, continuous) is another promising direction. Critically, future work should aim to connect parameters (μ,κ,γ̄) directly to downstream performance and failure modes in real-world applications, especially agentic systems.

7 Conclusion

In conclusion, we have presented a hierarchical Bayesian framework that provides a necessary advancement over traditional point-estimate evaluations, offering a statistically principled and interpretable lens through which to understand LLM performance. By dissecting performance into accuracy, stochasticity, and heterogeneity, and quantifying the associated uncertainties, it enables a comprehensive assessment of LLM capabilities and robustness. This richer, multi-faceted perspective is critical for guiding the development of the increasingly complex landscape of generative and agentic AI.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fams.2026.1854495/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Ronald Wesonga, Sultan Qaboos University, Oman

Reviewed by: Manoj Kumar, Sultan Qaboos University, Oman

Kyosuke Bunji, Kobe University, Japan