{"id":"d4b2e194-08f2-4c6a-a242-a7895a6cab34","arxiv_id":"2502.04990","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":0,"one_line_summary":"Precomputing sufficient statistics in Stan speeds up Bayesian sampling for several popular models, with gains that grow as the number of observations increases.","lead":"This paper shows that precomputing sufficient statistics makes Bayesian MCMC in Stan much faster for linear regression, mixed effects, and factor models. The authors provide Stan code and benchmarks where the speed gap over brms and rstanarm grows with dataset size.","discovery_kind":"new_application","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Wall-clock speedups for fixed draw counts ignore mixing differences; the mixed-effects comparison even changes parameterization, so per-ESS gains are unverified.","rationale":"The paper's technical core is credible: the sufficient-statistics log-likelihood derivations in Sections 3.1-3.3 are algebraically sound (apart from typos such as d vs p in Section 3.3 and the prior-scale mismatch between the text and Listing A.3), the code is provided, and the posterior overlays suggest correct targets. The remaining doubt is whether the headline speedups measure accurate inference rather than merely faster iterations. The reader identified this as the weakest assumption, and I agree: fixed-draw wall-clock time conflates per-iteration cost with total cost to accuracy. The mixed-effects comparison strengthens the concern because it changes parameterization between the baseline and the proposed implementation, so even the per-draw mixing rates are not obviously comparable. This is an addressable empirical gap, not a flaw in the sufficient-statistics idea itself; a matched-ESS comparison would settle it. Since the reader's CONDITIONAL verdict already reflects this concern, no change in verdict is needed.","tokens_in":21892,"tokens_out":6081,"duration_ms":69153,"concrete_test":"Re-run the Section 3.1-3.4 timing comparisons with each implementation run until all scalar parameters of interest reach a target effective sample size (e.g., ESS > 400 with R-hat < 1.01), then compare wall-clock time to reach that target ESS, reporting ESS per second for every method. If the sufficient-statistics implementations no longer dominate after this accuracy-matched comparison, the fixed-draw speedups in Figures 1-4 do not support the central claim.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The paper's central empirical claim is that sufficient statistics make Bayesian computation faster, but the speed comparisons in Figures 1-4 measure wall-clock time for a fixed number of posterior draws, not time to a given accuracy. This assumes all implementations mix at comparable rates per draw. That assumption is particularly insecure in the mixed-effects example: the vectorised baseline (Listing A.5) uses centered random effects, z_1 ~ N(0, sd_1), while the sufficient-statistics implementation (Listing A.7) uses the non-centered parameterization r_1_1 = sd_1 * z_1 with z_1 ~ N(0,1). These parameterizations have different HMC geometry and can require very different effective sample sizes for the same posterior accuracy. No ESS, R-hat, or Monte Carlo standard error is reported for any timing run; Figures A.1-A.4 show only marginal posterior overlays for one setting, which do not establish comparable convergence or mixing. If the sufficient-statistics chains need more effective samples per draw, the reported speedups shrink proportionally. The central claim 'faster Bayesian computation' therefore rests on an unverified equal-mixing assumption. A secondary fairness issue is that Section 3.1 states the common prior is sigma ~ t3(0, 3.7), but the brms code in Listing A.3 uses t3(0, 3.2), so the linear-regression baselines may not even target identical posteriors.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes that, when a Stan model's likelihood can be written in terms of precomputed sufficient statistics, each MCMC likelihood evaluation becomes independent of the number of observations n, and that this yields large practical speedups over existing Stan implementations. The authors demonstrate the idea on Gaussian linear regression with non-conjugate priors, linear mixed effects models, factor analysis models, and a Poisson regression model where only partial sufficient statistics are available. They compare their implementations against vectorised Stan code, brms, and rstanarm, measuring wall-clock time for a fixed number of posterior draws. Code for all implementations is provided in the appendix and online.","tokens_in":22156,"tokens_out":4277,"duration_ms":48166,"significance":"If the central claim holds, the paper offers a genuinely useful and easy-to-adopt practical improvement for Stan users: the sufficient-statistic reformulations are simple, the code is provided, and the potential savings grow with n. The strengths of the paper are that the Stan code is complete and reproducible, the comparisons use standard external baselines, and the algebraic likelihoods in the code appear correct. However, the headline conclusion that the method is 'faster Bayesian computation' is currently supported only by wall-clock timings for fixed draw counts, without evidence that the compared chains achieve comparable accuracy per draw. Because of this, the practical significance is conditional on additional convergence and mixing diagnostics.","major_comments":[{"comment":"The central empirical claim is measured by wall-clock time to produce a fixed number of posterior draws, not by time to a given estimation accuracy. No effective sample size, R-hat, or Monte Carlo standard error is reported for any timing run, so the reported speedups assume that all implementations mix at comparable rates per draw. This assumption is especially insecure in the mixed-effects example: the vectorised baseline in Listing A.5 uses the centered parameterization z_1 ~ N(0, sd_1), while the sufficient-statistics implementation in Listing A.7 uses the non-centered parameterization r_1_1 = sd_1 * z_1 with z_1 ~ N(0,1). These parameterizations have different HMC geometry and can require substantially different numbers of effective samples. Please report ESS and convergence diagnostics for every implementation and scenario, and ideally compare time to a target ESS or to a target Monte Carlo error rather than time for a fixed draw count.","section":"Sections 3.1-3.4, Figures 1-4"},{"comment":"The priors used by the baselines are not identical to the priors stated in the text. Section 3.1 states that σ ~ t3(0, 3.7) for all models, but the brms code in Listing A.3 uses σ ~ t3(0, 3.2), while the vectorised and sufficient-statistics code use σ ~ t3(0, 3.7). This means the brms baseline targets a slightly different posterior, and the posterior overlay in Figure A.1 cannot establish that the four implementations sample from the same distribution for σ. Please either use the same prior in all implementations or explicitly document the discrepancy and assess whether it affects the timing comparisons.","section":"Section 3.1, Listings A.2-A.4"},{"comment":"Section 3.2 states that the common prior is σ, σ_u ~ t3(0, 3.7), but Listings A.5, A.6, and A.7 all use t3(0, 2.5) for both σ and sd_1. The paper acknowledges that rstanarm uses a different prior on σ_u, but the code-to-text mismatch for the other baselines is unexplained. Since the mixed-effects comparison also changes the parameterization of the random effects between the vectorised and sufficient-statistics implementations, the per-draw speedup cannot be interpreted without first establishing that the targets and the mixing behavior are comparable.","section":"Section 3.2, Listings A.5-A.7"},{"comment":"The displayed sufficient-statistic likelihood equations contain errors that make them difficult to verify against the code. In Section 3.2, the expression contains the term '- n⊤1 / 2 log(2πσ^2)', which appears to be a typographical corruption of a factor involving n (or the vector of group sizes); it should be written as -(sum_j n_j)/2 * log(2πσ^2). In Section 3.3, the factor-model log-likelihood is written as n/2 (-d log(2π) + |Ω| - tr(SΩ)); this should be -(n p/2) log(2π) + (n/2) log |Ω| - (n/2) tr(SΩ), i.e., the dimension is p, not d, and the determinant term is missing the logarithm. The Stan code in Listings A.9 and A.10 contains the correct expression, but the displayed equations need to be corrected so that the mathematical claim matches the implementations.","section":"Sections 3.2 and 3.3"}],"minor_comments":[{"comment":"In the displayed expression for the log-likelihood, the sufficient statistic is written as s(θ), but it should be s(y), since it is a function of the data only.","section":"Section 2.1"},{"comment":"The Poisson regression data block declares vector[N] Y, although the response is a count variable. Using array[N] int Y would match the data type and avoid an implicit real conversion; the subsequent sufficient-statistic algebra is correct either way.","section":"Section 3.4, Listing A.13"},{"comment":"There are several typographical issues: the author name 'Bürkner' appears as 'Burkner' or 'B¨ urkner' in the text and references, the Poisson distribution is typeset as 'P oi' in Sections 3.4 and A.4, and the phrase 'rstanarm' is sometimes split across a line break. These are cosmetic but should be fixed.","section":"Throughout"},{"comment":"The caption says the figure overlays posteriors from 'different Stan implementations', but the legend shows only suff-stat, Wood, and vect; please make the caption consistent with the three implementations actually compared for the factor model.","section":"Figure A.3"}],"recommendation":"major_revision","confidential_remarks":"The paper addresses a practical and currently underserved aspect of Stan usage, and the provided code is a real resource for practitioners. The main risk is not the sufficient-statistics algebra but the fairness of the empirical speed comparison: fixed-draw wall-clock times without ESS or convergence diagnostics cannot establish 'faster Bayesian computation' when the compared models use different parameterizations and, in places, different priors. These issues are fixable within the scope of the paper, so I see major revision rather than rejection as the appropriate outcome."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Short version: this is a useful, honest engineering paper. The sufficient-statistics trick is classical, but the authors actually implement it in Stan for four model classes, ship the code, and benchmark against brms, rstanarm, and vectorised Stan. The speedups at large n are real, and the posterior overlays give reasonable evidence that the implementations target the same distribution. The math checks out modulo typos.\n\nThe new thing is not the algebra; it is the demonstration that the current defaults in major R packages leave a lot of performance on the table. That is worth publishing.\n\nSoft spots, in order of importance. First, the benchmarks time a fixed number of draws and do not report ESS or R-hat. For the linear, factor, and Poisson examples the competing implementations use the same parameterization, so equal mixing per draw is a safe assumption. But the mixed-effects comparison is confounded: the vectorised baseline (Listing A.5) uses a centered parameterization, while the sufficient-statistics version (Listing A.7) uses non-centered. That alone can change HMC efficiency. The brms comparison is fair (both non-centered), but the vectorised comparison is not. The paper should add ESS/R-hat and ideally re-run the vectorised baseline in non-centered form.\n\nSecond, there are several typos and inconsistencies between the text and the code. Section 3.1 states the prior is sigma ~ t3(0, 3.7) but Listing A.3 has 3.2; Section 3.2 text repeats 3.7 while Listings A.5-A.7 use 2.5. The factor-model log-likelihood in Section 3.3 writes -d log(2pi) where it should be -p log(2pi), and the linear-regression log-likelihood drops the -n/2 log(2pi) constant without comment. None of these affect the conclusions, but a reader cannot resolve them without reading the code.\n\nThird, the paper never quantifies the overhead of computing the sufficient statistics in transformed data. For the mixed-effects example that is an O(n) pass, so the \"constant in n\" claim is loose; the figures show this, but the text should be explicit.\n\nThe Poisson example is honestly partial and shows modest gains, which is good.\n\nBottom line: this deserves peer review. The central claim holds, the code is there, and the issues are fixable. I would send it out, asking the referees to focus on the benchmarking protocol.","headline":"A practical, honest demonstration that sufficient statistics make Stan much faster for several standard models; the benchmarks need diagnostics and a fairer mixed-effects comparison, but the core claim holds.","tokens_in":22649,"tokens_out":4250,"would_cite":false,"duration_ms":45292,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["62F15","62J05","62H25","65C05"],"pacs":[],"model":"deepseek-v4-flash","headline":"Precomputing sufficient statistics makes each MCMC likelihood evaluation independent of the number of observations, cutting runtime dramatically in common Bayesian models.","keywords":["sufficient statistics","exponential family","Bayesian computation","Markov chain Monte Carlo","Stan","probabilistic programming","hierarchical models","factor analysis"],"falsifier":"Run the same four implementations with convergence diagnostics reported per unit time: compute effective sample size per second and require all chains to pass $\\hat{R}<1.01$. If the sufficient-statistics implementations need substantially more draws to reach the same effective sample size, the wall-clock speedups shrink or disappear; this can be checked with the supplied Stan code on the paper's simulated data.","tokens_in":21708,"feed_emoji":"⚡","tokens_out":7502,"duration_ms":72731,"temperature":0.7,"pith_summary":"The paper argues that a classical idea—writing an exponential-family likelihood in terms of its sufficient statistics—has been left unused by leading probabilistic programming tools, and that using it produces large practical speedups in Bayesian computation. For Gaussian linear regression, hierarchical mixed effects models, and factor analysis, the likelihood is first compressed into a few precomputed summaries such as $X^\\top X$, $X^\\top y$, and $y^\\top y$; each MCMC iteration then evaluates the log-likelihood without looping over the $n$ data points. The paper reports that this implementation produces the same posteriors as current Stan-based packages while taking a fraction of the wall-clock time, with the gain growing as $n$ grows. A fourth example, Poisson regression, shows that even when the likelihood can only be partially rewritten in terms of sufficient statistics, moderate gains remain. The broader point is that the bottleneck in modern Bayesian software is often the unnecessary per-observation cost of likelihood evaluation, not the sampling algorithm itself.","feed_headline":"Precomputed summaries keep Bayesian sampling cost flat as n grows","feed_subtitle":"Rewriting exponential-family likelihoods with precomputed summaries outperforms brms, rstanarm, and vectorized Stan code.","key_machinery":"The load-bearing object is the sufficient-statistics reparameterisation of an exponential-family likelihood: if $f(y;\\theta)=h(y)\\exp\\{\\eta(\\theta)^\\top T(y)-A(\\theta)\\}$, then the joint log-likelihood is $\\eta(\\theta)^\\top s(y)-nA(\\theta)$ plus a data-only constant, with $s(y)=\\sum_i T(y_i)$. Precomputing $s(y)$ in the transformed data block of Stan moves the $O(n)$ sum out of the MCMC loop; for linear regression the summaries are $S_{xx}=X^\\top X$, $S_{yx}=y^\\top X$, and $S_{yy}=y^\\top y$, for mixed effects models they include per-group counts and sums, and for factor models the summary is the sample covariance $S=y^\\top y/n$. An optional Woodbury identity computes the precision matrix $\\Omega=(\\Lambda\\Lambda^\\top+\\mathrm{diag}(\\psi))^{-1}$ in $O(d^3)$ rather than $O(p^3)$ operations, which matters when the number of factors $d$ is much smaller than the observation dimension $p$.","core_discovery":"On the paper's own terms, the discovery is that for any model whose likelihood is an exponential family, the log-likelihood for $n$ iid observations can be written as $\\eta(\\theta)^\\top s(y) - nA(\\theta)$, where $s(y)=\\sum_i T(y_i)$ is the vector of sufficient statistics. Since $s(y)$ depends only on the data, it can be computed once in Stan's transformed data block. After that, each evaluation of the likelihood costs the same regardless of $n$, because all $n$-dependent work is absorbed into the precomputed summaries. The paper supplies Stan programs that do this for linear regression with Student-$t$ and Cauchy priors, linear mixed effects models, and factor models (using the sample covariance matrix $S=Y^\\top Y/n$ and, in one variant, a Woodbury decomposition of the precision matrix), and shows the sampled posteriors match those of brms, rstanarm, and a vectorized implementation while the runtime either stays flat or grows far more slowly with $n$. For Poisson regression, where the term $\\sum_i \\exp(x_i^\\top\\beta)$ cannot be precomputed, the paper still reports gains from precomputing $X^\\top y$ and vectorizing the remaining sum.","pith_inferences":["My inference: the same precomputation applies to any other exponential-family likelihood with non-conjugate or hierarchical priors—logistic, gamma, negative binomial, and multinomial regression—so the demonstrated speedups should generalise beyond the three examples.","My inference: the one-time cost of forming $X^\\top X$ and $X^\\top y$ scales like $O(np^2)$ in dense form, so the method's advantage is largest when $n$ is large and $p$ is moderate; for very high $p$, that precomputation, not the MCMC likelihood evaluation, could become the bottleneck.","My inference: a natural stress test is to compare effective sample size per second, not wall-clock time per number of draws, because the sufficient-statistics reformulation changes the HMC target's gradient and could in principle alter mixing; if mixing per draw is worse, part of the reported speedup would be offset.","My inference: the partial-sufficient-statistics idea in Poisson regression suggests a general recipe for GLMs where only an additive nonlinear term resists summarisation—precompute what can be summarised and vectorise the rest—which could extend to other nonlinear link functions."],"forward_implications":["For Gaussian linear regression with non-conjugate priors, time to produce 5000 post-warm-up draws stays essentially constant as $n$ rises from 100 to 10,000, while brms, rstanarm, and vectorised code all slow down considerably.","For mixed effects models, precomputed per-group summaries make the sufficient-statistics implementation scale more slowly with $n$ and keep it faster than brms and rstanarm as the number of groups $J$ varies.","For factor models, replacing the raw data by the sample covariance matrix $Y^\\top Y/n$ speeds up sampling, and adding the Woodbury decomposition yields further gains that grow with observation dimension $p$.","Even when only partial sufficient statistics exist, as in Poisson regression, the paper reports that the implementation using $X^\\top y$ plus vectorised $\\sum_i \\exp(x_i^\\top\\beta)$ is the fastest of the four compared methods.","If such implementations were folded into brms and rstanarm, users would get these speedups without changing how they write models."],"supporting_citations":[{"why":"Defines Stan, the probabilistic programming language whose implementations the paper speeds up.","marker":"Carpenter et al., 2017"},{"why":"Supplies the brms implementation that serves as a baseline in the runtime comparisons.","marker":"Bürkner, 2017"},{"why":"Supplies the rstanarm implementation that serves as a baseline in the runtime comparisons.","marker":"Goodrich et al., 2024"},{"why":"Provides the vectorised reference implementation from the Stan User Guide used as another baseline.","marker":"Stan Development Team, 2024"},{"why":"Provides the original Stan factor model which the paper extends and compares against.","marker":"Farouni, 2015"},{"why":"Supplies the microbenchmark package used to measure all wall-clock runtimes.","marker":"Mersmann, 2024"},{"why":"Source of the Woodbury decomposition used to invert the factor-model covariance more cheaply.","marker":"Ghahramani et al., 1996"}],"fun_headline_variants":["Sufficient statistics slash Bayesian sampling time in Stan","Precompute summaries, sample fast: sufficient stats in Stan","Flat-cost Bayesian inference via sufficient statistics","Sufficient statistics speed up Stan for large data","Constant-cost Bayesian sampling with sufficient statistics"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The speedups are measured as wall-clock time for a fixed number of posterior draws, and the comparison assumes all implementations converge and mix at comparable rates per draw; if the sufficient-statistics chains need effectively more iterations to reach the same accuracy—or if the differing prior on $\\sigma_u$ in the rstanarm mixed-effects model biases that baseline—the reported gains would change.","fun_headline_variants_meta":{"raw":{"variants":["Sufficient statistics slash Bayesian sampling time in Stan","Precompute summaries, sample fast: sufficient stats in Stan","Flat-cost Bayesian inference via sufficient statistics","Sufficient statistics speed up Stan for large data","Constant-cost Bayesian sampling with sufficient statistics"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000283,"raw_usage":{"total_tokens":1711,"prompt_tokens":1021,"completion_tokens":690,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":637,"completion_tokens_details":{"reasoning_tokens":621}},"tokens_in":637,"tokens_out":690,"duration_ms":8024,"temperature":1.0,"reasoning_tokens":621,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-08T20:41:36.764073+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the same four implementations with convergence diagnostics reported per unit time: compute effective sample size per second and require all chains to pass $\\hat{R}<1.01$. If the sufficient-statistics implementations need substantially more draws to reach the same effective sample size, the wall-clock speedups shrink or disappear; this can be checked with the supplied Stan code on the paper's simulated data.","supporting_citations":[{"cited_title":", author Gelman, A","cited_arxiv_id":null,"evidence_quote":"Defines Stan, the probabilistic programming language whose implementations the paper speeds up."},{"cited_title":", author Gabry, J","cited_arxiv_id":null,"evidence_quote":"Supplies the rstanarm implementation that serves as a baseline in the runtime comparisons."},{"cited_title":"title Stan functions reference","cited_arxiv_id":null,"evidence_quote":"Provides the vectorised reference implementation from the Stan User Guide used as another baseline."},{"cited_title":", year 2015","cited_arxiv_id":null,"evidence_quote":"Provides the original Stan factor model which the paper extends and compares against."},{"cited_title":", year 2024","cited_arxiv_id":null,"evidence_quote":"Supplies the microbenchmark package used to measure all wall-clock runtimes."},{"cited_title":", author Hinton, G.E","cited_arxiv_id":null,"evidence_quote":"Source of the Woodbury decomposition used to invert the factor-model covariance more cheaply."}],"review_version":1}