{"id":"ccd12c05-10f0-44bd-a007-93f6aa636dba","arxiv_id":"2411.09851","paper_version":4,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":5.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":3,"one_line_summary":"Symbolic regression with a re-optimization step automates parametric modeling of binned high-energy physics data and provides uncertainty estimates.","lead":"SymbolFit is a new software framework that uses symbolic regression to automatically search for closed-form functions that fit binned data, and then estimates uncertainties on the fitted parameters. The authors demonstrate it on five real CMS datasets and several toy examples, showing fits comparable to the hand-crafted empirical functions used in published LHC searches.","discovery_kind":"new_application","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The uncertainty estimates from the least-squares covariance are never calibrated, so the paper's main added value over plain symbolic regression is unvalidated.","rationale":"The reader's weakest assumption already identifies exactly this concern: the uncertainty bands are derived from a nonlinear least-squares covariance on a frozen functional form, with no coverage calibration. My stress-test pass finds no reason to move away from that judgment. The paper's genuine strengths include a released API, demonstration on five real LHC datasets, blinded signal-injection tests with out-of-sample validation, and the fact that the ROF step yields noticeably better chi2 scores than the raw SR output. These support the claim that SymbolFit can find good point fits. They do not, however, support the claim that the uncertainty estimates are valid for statistical inference, because coverage is a property of the interval, not of the best-fit curve, and no interval property is tested. The absence of a known true distribution for both the toy and real datasets is precisely why a synthetic-truth or bootstrap calibration study is necessary. I therefore keep the reader's CONDITIONAL verdict: the framework is promising and worth presenting, but the uncertainty-calibration claim needs a dedicated validation before the tool is used for production LHC fits. The selection issue is real but secondary: the paper openly presents candidate selection as a user task, whereas the uncertainty estimate is advertised as an automatic output, making calibration the more load-bearing gap.","tokens_in":30566,"tokens_out":3358,"duration_ms":39255,"concrete_test":"For each of the five LHC datasets, take the fitted CMS empirical function (or, where unavailable, a kernel-smoothed template) as the synthetic truth, generate 500 pseudo-experiments with Poisson fluctuations, run SymbolFit with the same configuration, and for each selected candidate compute the empirical coverage of the 68% parameter-sampling band across all bins. If the mean coverage deviates from 68% by more than 10 percentage points, the uncertainty estimates are not calibrated.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim that SymbolFit 'provides uncertainty estimates' rests on the assumption that the covariance matrix from the LMFIT re-optimization fit (Sec. 4, step iv) yields valid parametric uncertainties for the frozen SR functional forms. This is a standard nonlinear least-squares covariance, which is calibrated only if the fitted function is close to the true data-generating process and if parameter nonlinearities or boundary effects are negligible. No coverage test is reported: the paper never checks whether the 68% bands shown in Figs. 7, 8, A3, A4, A8, A12, and A16 actually contain the true distribution at the claimed rate. The toy datasets are explicitly generated without an underlying function, so they cannot serve as ground truth; for the real LHC datasets the true background is unknown. The signal-injection studies validate signal extraction but do not validate background uncertainty coverage, because they assess point estimates of injected-signal parameters rather than the coverage of the background-model bands. If the bands are miscalibrated, downstream inference with Combine or pyhf using these uncertainties would be biased, which undermines the primary motivation stated in Sec. 3: 'parametric models without well-defined uncertainties cannot be used in the statistical inference workflows within HEP.' The ROF procedure is a reasonable heuristic, but the paper provides no evidence of its calibration across the wide range of candidate functions it recommends.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper introduces SymbolFit, a Python framework that combines symbolic regression (via PySR) with a nonlinear least-squares re-optimization step (via LMFIT) to produce parametric fits of binned data together with uncertainty estimates. The framework is demonstrated on several toy datasets and on five real LHC invariant-mass spectra (dijet, diphoton, trijet, paired-dijet, and dimuon), where the resulting functions are compared with the empirical functions used by CMS. In addition, the authors perform signal-injection studies: background models are fitted on pseudodata with the signal region blinded, and then used to extract an injected Gaussian signal from the unblinded spectrum. The central claim is that a simple, largely configuration-free SR setup can automate the manual process of empirical function selection while also providing uncertainty estimates suitable for downstream statistical inference.","tokens_in":30719,"tokens_out":4326,"duration_ms":45406,"significance":"If the claimed capabilities hold, SymbolFit would be a genuinely useful tool for HEP data analysis, potentially reducing the manual effort in background and signal modeling. The paper's strengths include a public API, the use of five real LHC datasets with comparison to official CMS empirical functions, genuine out-of-sample validation through blinded signal injection, and a demonstration on multi-dimensional toy data. The re-optimization step clearly improves the initial SR fits. However, the central added value over plain symbolic regression—uncertainty estimation—is not validated by any coverage or calibration test, which is a substantial gap given the intended use in statistical inference.","major_comments":[{"comment":"The 68% uncertainty bands are obtained by sampling parameters from the LMFIT covariance matrix of the frozen functional form, but the paper never checks whether these bands actually cover the true distribution at the claimed rate. The toy datasets are explicitly described as generated 'without reference to an underlying function' (Sec. 5.1, Fig. 3 caption), so they cannot serve as ground truth for coverage; for the real LHC datasets the true background is unknown. The signal-injection studies validate signal extraction, but they assess point estimates of injected signal parameters and do not test the coverage of the background-model bands themselves. Since the abstract claims that the framework provides 'uncertainty estimates in a single run' and Sec. 3 states that parametric models without well-defined uncertainties cannot be used in HEP statistical inference, the lack of any calibration test for the uncertainty bands is a load-bearing gap. The authors should add ensemble-based coverage checks (for example, fitting to pseudodata generated from a known analytic function and measuring the empirical coverage of the reported 68% intervals, or performing posterior-predictive checks on toy data with a known ground truth).","section":"Sec. 4"},{"comment":"The claim that the same fit configuration can be applied to distributions of different shapes, with only the random seed varied, is contradicted by the experimental setup. For the diphoton and dimuon datasets, the manuscript states that the maximum complexity is reduced from 80 to 20 (A.2 and A.5), and Sec. 5.2 notes 'variations in the maximum complexity values' for Toy Dataset 2. This means the configuration is not actually identical across all demonstrations. Since the abstract and Sec. 3 explicitly rest the flexibility claim on a single simple configuration, the authors should either (i) restrict the claim to the operator set and loss function, treating maxsize as a user-tunable parameter, or (ii) provide at least one additional demonstration where a fully identical configuration, including maxsize, is used across several different distribution shapes.","section":"Sec. 5.2"},{"comment":"The signal extraction validation is reported only qualitatively: the text states that fitted signal parameters are 'within reasonable uncertainties' and shows residual plots, but no quantitative pull distributions, bias values, or coverage statistics are provided. For a method intended to feed into statistical inference, the calibration of the fitted signal parameters (e.g., the mean and width of the pull distribution, or the fraction of fits where the true injected value lies within the reported 1σ interval) should be quantified and compared with the empirical model. This would also provide indirect evidence about whether the uncertainties from the SR models are trustworthy.","section":"Sec. 5.2"}],"minor_comments":[{"comment":"The text refers to 'Fig. A15' for the dijet signal extraction study, but the corresponding figure is Fig. 11; the same figure number is later used for the paired-dijet extraction in Appendix A.4, which will confuse readers.","section":"Sec. 5.2"},{"comment":"The transformation for Toy Dataset 1 is given as f(x) → 165 × f(0.00211(x − 12.5)), but the reader is not told what the variable x represents in the fitted, scaled dataset; defining x as the scaled mass variable explicitly would improve reproducibility.","section":"Sec. 5.1"},{"comment":"The notation for the chi-square statistic is inconsistent across tables and text (for example, 'χ2/NDF' and 'χ2/NDF'); please unify the formatting.","section":"All tables"},{"comment":"The phrase 'ratio of the 68% quantile range to the mean' is ambiguous; it should be clarified as the relative uncertainty, i.e., the width of the band divided by the mean function value.","section":"Figure captions"},{"comment":"The abstract states that the framework provides uncertainty estimates 'in a single run', but the demonstrated workflow uses multiple random seeds and selects a candidate function from each batch; please clarify that a single configuration, run multiple times with different seeds, is what is meant.","section":"Abstract"}],"recommendation":"major_revision","confidential_remarks":"The paper is a potentially useful contribution to the HEP parametric-modeling toolkit, and the authors have made a genuine effort to validate signal extraction with blinded injections. The main concern is that the uncertainty estimates—the central advance over plain symbolic regression—are never calibrated. This is a fixable issue but it is load-bearing; I would not reject the paper outright, but the revision needs to add quantitative coverage studies and temper the universal-configuration claim. The paper might fit better in a journal like JINST than in a pure physics-letters venue, but that is an editorial decision."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"The authors have put together a genuinely useful paper. The core idea is not just symbolic regression, which is well known, but the combination of SR with a re-optimization step that turns raw SR output into a normal parametric model with a covariance matrix, so it can drop into standard HEP statistical tools. That integration, plus the breadth of validation, is the new contribution. They run the same basic PySR configuration on five real LHC datasets and several toy datasets, fit the blinded background with injected signals, and show that the SR models extract signal parameters about as well as the CMS empirical functions. That is real work, and shipping the code makes it reproducible.\n\nThe main soft spot is exactly what your stress-test note says: the uncertainty bands are never calibrated. The green bands in Figs. 7, 8, and the appendices are produced by sampling parameters from the LMFIT covariance. That is a legitimate way to propagate parameter uncertainties if the functional form is correct, but here the functional form is selected post hoc from dozens of candidates, and no coverage test is performed. The toy datasets are deliberately generated without an underlying function, so they cannot serve as ground truth; the real datasets have no known true background. So the claim that SymbolFit 'provides uncertainty estimates' is true in the sense that it produces numbers, but the paper does not show that those numbers behave like 68% intervals. For HEP workflows that use these uncertainties in a likelihood, that matters. A calibration study on toys with known truth, or a comparison to a nonparametric method like Gaussian processes, would close most of the gap.\n\nTwo smaller issues. The selection of the 'SR model 1/2/3' candidates from each batch is not described. That introduces a selection effect that could bias the uncertainty estimates, and it makes the comparison to the CMS empirical model less crisp. Also, the 'same fit configuration' claim is slightly overstated: maxsize is 80 for dijet and 20 for diphoton/dimuon. Minor quibble, not a flaw.\n\nOverall, this is a solid methods paper for HEP analysts who do background modeling. It deserves a serious referee. I would ask for the calibration study and a clearer candidate-selection protocol before publication, but I would not desk-reject it. It is a real contribution, and the authors have done the work to back it up.","headline":"A genuinely useful, well-tested tool for automating parametric fits in HEP, with the main caveat that the uncertainty bands are never calibrated against coverage.","tokens_in":31357,"tokens_out":2754,"would_cite":true,"duration_ms":28499,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"SymbolFit claims that one generic symbolic-regression configuration can automatically produce parametric background models, complete with uncertainties, that match hand-tuned empirical fits across five LHC search spectra.","keywords":["symbolic regression","parametric modeling","background modeling","uncertainty estimation","binned data","high-energy physics","genetic programming","LHC searches"],"falsifier":"Run the framework many times on pseudodata sampled from a known smooth distribution that is deliberately outside the operator set, then count how often the true curve lies inside the reported 68% uncertainty band. Empirical coverage substantially below 68% would show that the covariance-based uncertainty model is miscalibrated, breaking the claim of simultaneous uncertainty estimation.","tokens_in":30290,"feed_emoji":"⚛️","tokens_out":7859,"duration_ms":70824,"temperature":0.7,"pith_summary":"The paper introduces SymbolFit, a framework that replaces the manual, iterative hunt for empirical functions in high-energy physics with symbolic regression, a genetic-programming search over candidate functional forms. Its central claim is that one generic fit configuration, with only the allowed operators and a complexity bound specified, can produce well-fitting background models across five different collision-search spectra (dijet, trijet, paired-dijet, diphoton, dimuon) and several toy datasets, while simultaneously returning parameter uncertainties through a re-optimization step. The authors argue that the results are comparable to the hand-crafted empirical models used in published searches, at a fraction of the manual effort. A sympathetic reader would care because parametric modeling is a bottleneck in many binned-data analyses, and this points to automating it end-to-end.","feed_headline":"One symbolic regression recipe fits five LHC spectra","feed_subtitle":"SymbolFit replaces manual function hunting, returning fits with uncertainties from one run across five search datasets.","key_machinery":"The load-bearing mechanism is genetic-programming symbolic regression, in which candidate functions are expression trees over a fixed operator set and the functional form itself is the search variable, followed by a re-optimization fit in which each discovered function is frozen and its constants re-fitted by nonlinear least squares; the covariance matrix of that fit provides the uncertainty estimates. The same search configuration, an operator set plus nesting constraints and a maximum complexity, is applied across datasets, with only the random seed changed to generate a diverse batch of candidates. A chi-square loss with separate up/down bin uncertainties guides both the search and the re-fit.","core_discovery":"The paper's core claim is that symbolic regression with a fixed, simple configuration can stand in for the traditional empirical procedure of guessing and refining a functional form. For each dataset, a single run of the framework returns a batch of candidate closed-form functions; each function is then frozen and its constant parameters re-fitted by nonlinear least squares, yielding best-fit values, correlations, and uncertainty bands. On blinded pseudodata built from the published dijet spectrum with an injected Gaussian signal, three runs with different random seeds produce background models whose goodness-of-fit and extracted signal parameters are comparable to the empirical model used in the original search. The authors claim the same configuration transfers to other distributions of different shapes without retuning, and that the closed-form output slots directly into standard statistical tools for hypothesis testing.","pith_inferences":["The paper does not calibrate the coverage of its uncertainty bands; an obvious test would be to generate pseudodata from a known smooth distribution, run the framework, and check whether the reported 68% intervals contain the true values at the claimed rate.","Because the search minimizes chi-square on the binned data, the candidate pool may overfit in regions of low statistics; selecting among candidates purely by goodness of fit could introduce selection bias, which the discrete-profiling treatment would mitigate.","A natural extension would be to treat the batch of candidate functions from a single run as a model-averaged background with empirically calibrated coverage, rather than picking one function.","The same recipe could generalize beyond collision physics to any binned data where a smooth empirical description with uncertainties is needed, such as detector corrections in other fields."],"forward_implications":["Background models for resonance searches can be produced automatically from binned data, removing the manual trial-and-error step that currently consumes significant analyst effort.","The same fit configuration can be reused across channels and analyses, so modeling does not need to be re-derived whenever detector conditions or event selections change the data shape.","Because the output is a batch of closed-form functions with uncertainties, the choice of function can be treated as a discrete systematic, and the templates integrate directly with existing hypothesis-testing tools.","Multi-dimensional binned data, such as scale-factor parameterizations in two variables, can be modeled without hand-constructing polynomial ansatze.","A single run routinely produces candidates spanning a range of fit qualities, giving downstream analyses flexibility in trading fit quality against complexity."],"supporting_citations":[{"why":"Supplies the genetic-programming symbolic regression engine that searches over candidate functional forms.","marker":"[5]"},{"why":"Provides the nonlinear least-squares re-optimization fit whose covariance matrix yields the parameter uncertainties and correlations.","marker":"[37]"},{"why":"Supplies the dijet search dataset and the empirical background model that the symbolic-regression results are benchmarked against.","marker":"[25]"},{"why":"Makes the dijet spectrum publicly available as the input data for the main demonstration.","marker":"[42]"},{"why":"Illustrates in the introduction the manual construction of an empirical background function (its Eq. 1) that the framework automates.","marker":"[20]"}],"fun_headline_variants":["One recipe, five spectra: SymbolFit auto-models LHC data","Symbolic regression replaces manual function hunting in LHC fits","Single symbolic regression run fits five LHC search spectra","Automated parametric modeling: SymbolFit fits diverse LHC shapes","SymbolFit: machine-search your fitting function, uncertainties included"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The uncertainty estimate is the covariance matrix of a least-squares re-fit of the frozen functional form, so the entire uncertainty model rests on the assumption that the discovered function family is close enough to the true distribution for that covariance to capture the modeling error; the paper never tests whether the reported 68% bands actually contain the true distribution at the claimed rate.","fun_headline_variants_meta":{"raw":{"variants":["One recipe, five spectra: SymbolFit auto-models LHC data","Symbolic regression replaces manual function hunting in LHC fits","Single symbolic regression run fits five LHC search spectra","Automated parametric modeling: SymbolFit fits diverse LHC shapes","SymbolFit: machine-search your fitting function, uncertainties included"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000441,"raw_usage":{"total_tokens":2249,"prompt_tokens":969,"completion_tokens":1280,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":585,"completion_tokens_details":{"reasoning_tokens":1196}},"tokens_in":585,"tokens_out":1280,"duration_ms":10349,"temperature":1.0,"reasoning_tokens":1196,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-12T20:13:44.701089+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the framework many times on pseudodata sampled from a known smooth distribution that is deliberately outside the operator set, then count how often the true curve lies inside the reported 68% uncertainty band. Empirical coverage substantially below 68% would show that the covariance-based uncertainty model is miscalibrated, breaking the claim of simultaneous uncertainty estimation.","supporting_citations":[{"cited_title":"Search for narrow trijet resonances in proton-proton collisions at $\\sqrt{s}$ = 13 TeV","cited_arxiv_id":"2310.14023","evidence_quote":"Illustrates in the introduction the manual construction of an empirical background function (its Eq. 1) that the framework automates."}],"review_version":1}