Organizations: Department of Electrical Engineering IIT Kanpur · Department of Computer Science and Engineering IIT Kanpur · Department of Physics IIT Kanpur · Department of Electrical Engineering KU Leuven
Estimating the symbolic or analytical form of probability density functions (PDFs) from observed samples is a fundamental challenge in statistical and computational modelling. This process is critical for deriving interpretable and generalizable relationships characterizing the underlying phenomenon. Traditionally, this estimation depends strongly on domain expertise and prior field-specific knowledge, with experts selecting appropriate functional forms or parametric families based on empirical evidence and theoretical understanding. The coefficients of these forms are then typically determined through parameter estimation. In this paper, we develop a framework for estimating symbolic expressions of unnormalized distributions from observed samples using domain-specific prior knowledge, such as the range of interactions and a predefined set of primitive functions. We integrate deep generative models with symbolic regression (SR), incorporating inductive biases, such as factorizing large distributions, to keep the problem tractable. The deep generative models we examine include likelihood-based models, viz., flow models, and score-based models. Experiments show the effectiveness of the proposed framework for estimating density functions for multivariate toy distributions as well as lattices from computational physics, namely, XY model and φ4 theory. When applied to the renormalization problem in φ4 theory, the proposed framework estimates compact symbolic approximations of the hamiltonian function at different scales directly from samples, yielding expressions that may be challenging to derive using traditional perturbative or analytic approaches in nonperturbative settings.
Figures & tables
Figure 1: Block diagram of the proposed framework. (A) Flow-based estimation. (B) Score-based estimation. (C) Integration-free score-based estimation.
Figure 2: Lattice Factorization. (a) True H(x) can be factorized with a 2×2 patch. (b) While inference with SR, different patch sizes provide different inductive biases. (c) Illustration for redundancy in representation.
θ←θ−η∇θLM(θ)
Algorithm 1 Symbolic Function Estimation for Unnormalized Distributions
Table 1: Results and estimated expressions h^([x1,x2];λ) for a 2-D Gaussian distribution using various methods. The mean vector is fixed at μ=[−1.5,1.5]T , with two different covariance matrices, Σ1 and Σ2 , as detailed in the table. The target functions are: (A) 0.5x12+0.5x22+1.5x1−1.5x2 , and (B) 0.67x12−0.67x1x2+0.67x22+3.0x1−3.0x2 . Incorrect terms are shown in gray. NF = Flow model, SM = Score matching.
Table 2: Results for Many Well distribution with various dimensions, d=2,4,8,16,32,64 . The true function is h(xpl)=x14−6x12−0.5x1+0.5x22 . Estimated terms that do not match the true function are shown in gray. NF = Flow model, SM = Score matching.
Table 3: Results for XY Model. The true function is h(xp(l))=−0.71cos(x0−x1)−0.71cos(x0−x2) . Estimated terms that do not match the true function are shown in gray. NF = Flow model, DSM = Denoising score model.
Table 4: Results for ϕ4 theory with various λ1,λ2 settings. The true expressions are: (A) for λ1=4,λ2=−4 , h(xp(l))=4x04−2x0x1−2x0x2 , and (B) for λ1=4,λ2=1 , h(xp(l))=4x04+5x02−2x0x1−2x0x2 . Incorrect terms are highlighted in gray. NF = Flow model, DSM = Denoising score model.
Table 5: Results for estimating h^(xp(l)) for ϕ4 theory with different patch sizes ∣p∣ for flow model. True h(xp(l))=4x04+5x02−2x0x1−2x0x2 . Incorrect terms are shown in gray.
Figure 3: Renormalization results for non-perturbative region, λ1=4 and λ2=1 .
32x32
16x16
8x8
Observable
% Overlap( ↑ )
EMD ( ↓ )
% Overlap( ↑ )
EMD ( ↓ )
% Overlap( ↑ )
EMD ( ↓ )
Magnetization, M
98.116
2.28×10−4
98.845
9.09×10−5
98.759
8.79×10−5
Absolute Magnetization, |M|
98.197
2.06×10−4
98.794
8.69×10−5
98.749
8.30×10−5
Magnetic Susceptibility
98.244
3.81×10−3
99.187
4.47×10−4
99.187
1.13×10−4
Table 6: Quantitative comparison of physical observables computed from coarse-grained data and samples generated from the estimated symbolic Hamiltonians at different renormalization scales. We use 100,000 samples to compute the observables.
Magnetization, M
Absolute Magnetization, |M|
Magnetic Susceptibility
Mean
Std Dev ( σ )
Mean
Std Dev ( σ )
Mean
Std Dev ( σ )
32×32
Coarsened Data
−1.82×10−6
0.0066
0.0053
0.0039
0.0441
0.0609
Generated Data
9.09×10−5
0.0068
0.0055
0.0041
0.0481
0.0681
16×16
Coarsened Data
−1.82×10−6
0.0066
0.0053
0.0039
0.0110
0.0152
Generated Data
−3.36×10−5
0.0067
0.0053
0.0040
0.0115
0.0161
8×8
Coarsened Data
−1.82×10−6
0.0066
0.0053
0.0039
0.0028
0.0038
Table 7: Observable statistics of the coarse-grained data and generated samples. Mean and standard deviation are reported for magnetization, absolute magnetization, and magnetic susceptibility across different lattice sizes.
Dataset
MSE (H,Hθ)
MSE (H,H^)
PySR
EQL
Many Well, N=2
0.01
0.00
0.00
Many Well, N=4
0.05
0.00
0.00
Many Well, N=8
1.60
0.00
0.23
XY; λ=1.4−1
9.79
4.83
6.74
ϕ4;λ1=4,λ2=−4
1.43
0.02
0.29
Table 8: Comparison of MSE (H,Hθ) , where Hθ is computed by flow model, and MSE (H,H^) , where H^ is estimated by SR methods (PySR and EQL) using flow model outputs.
Table 9: EQL configurations and inductive biases used across benchmarks. Numbers in parentheses denote the number of neurons associated with each symbolic operator. The final column indicates whether the primitive library contains the exact functional terms present in the target Hamiltonian.
Experiment
Operator Library
Population
Complexity
Patch Size
Range of Interaction
Exact Terms Present
Multivariate Gaussian
{+,−,∗}
10
20
–
–
No
Many-Well
{+,−,∗,pow2}
30
30
2×1
2×1
Yes
ϕ4
{+,−,∗,pow2,pow3,pow4}
30
30
2×2 , 3×3 , 4×4
2×2
Yes
XY
{+,−,∗,sin,cos}
30
30
2×2 , 3×3
2×2
Yes
Appendix
Table 10: PySR configurations and inductive biases used across benchmarks.
No. of Samples
MSE( ↓ )
R2(↑)
RCE( ↓ )
ST( ↓ )
P( ↑ )
R( ↑ )
F1( ↑ )
Estimated Equation
10K
0.47
0.994
0.09
0
1
1
1
3.5x04−1.85x0x1−1.85x0x2
50K
0.07
0.999
0.03
0
1
1
1
4.2x04−2.04x0x1−2.04x0x2
100K
0.01
0.999
0.02
0
1
1
1
3.9x04−1.95x0x1−2.00x0x2
Appendix
Table 11: Effect of sample size on model performance.
Table 12: Symbolic estimation results for the ϕ4 theory across 20 different random seed settings using NF + PySR. The ground-truth local Hamiltonian for λ1=4 and λ2=1 is h(xp(l))=4x04+5x02−2x0x1−2x0x2 . Incorrect terms are highlighted in gray. NF denotes the flow-based model.
Table 13: Symbolic estimation results for the ϕ4 theory across 5 different random seed settings using NF + EQL. The true expressions is h(xp(l))=4x04+5x02−2x0x1−2x0x2 . Incorrect terms are highlighted in gray.
MSE( ↓ )
R2(↑)
RCE( ↓ )
ST( ↓ )
P( ↑ )
R( ↑ )
F1( ↑ )
Estimated Equation h^(xp(l))
1
0.72
0.985
0.279
0
1.00
1.00
1.00
1.58x14−7.27x12−0.65x1+0.51x22
2
0.46
0.990
0.126
0
1.00
1.00
1.00
1.37x14−6.64x12−0.51x1+0.50x22
3
0.32
0.993
0.242
1
0.80
1.00
0.89
1.68x14−7.46x12−0.52x1+0.50x22
4
0.28
0.994
0.160
0
1.00
1.00
1.00
1.46x14−6.96x12−0.50x1+0.51x22
5
0.40
0.991
0.179
0
1.00
1.00
1.00
1.44x12−6.92x12−0.45x1+0.51x22
6
0.42
0.991
0.189
0
1.00
1.00
1.00
1.49x14−7.02x12−0.55x1+0.50x22
Appendix
Table 14: Symbolic estimation results for the MW-64 benchmark across 18 random seed settings using SM + PySR. The true local Hamiltonian is h(xpl)=x14−6x12−0.5x1+0.5x22 . Incorrect terms are highlighted in gray.
Table 15: Symbolic estimation results for the MW-64 benchmark across 8 random seed settings using SM + EQL. The true local Hamiltonian is h(xpl)=x14−6x12−0.5x1+0.5x22 . Incorrect terms are highlighted in gray.
0.26x_{12}^{4}+1.62x_{12}^{2}+2.07x_{2}x_{12}+1.47x_{10}x_{12}-2.11x_{12}x_{14}-1.97x_{12}x_{22}+\cdots\text{{\color[rgb]{0.5,0.5,0.5}Higher order terms}}
Table 16: Results for a synthetic long-range interaction system. The true expression is h(xp(l))=4x124+5x122−x2x12−x10x12−x12x14−x12x22 Incorrect terms are highlighted in gray. NF = Flow model, DSM = Denoising score model.
Table 17: Symbolic estimation on the scalar ϕ4 benchmark under primitive library misspecification. Explicit polynomial primitives are removed from the symbolic library, leaving only arithmetic and trigonometric operators.
Table 18: Symbolic estimation results for the XY model under primitive library misspecification. The symbolic regressor approximates the cosine interaction indirectly through polynomial compositions.
Method
MSE(H,Hθ)
MSE(H,H^)
ϕ4;λ1=4,λ2=1
NF + PySR
0.22
0.01
NF + EQL
2.29
XY; λ=1.4−1
NF + PySR
9.79
8.84
NF + EQL
7.13
Appendix
Table 19: Error propagation analysis under primitive library misspecification. MSE(H,Hθ) denotes the error of the learned flow model relative to the true Hamiltonian, while MSE(H,H^) denotes the error of the estimated symbolic approximation.
Discrete probability laws underpin statistical modeling, yet the catalog of interpretable distributions has expanded only gradually through centuries of case-by-case mathematical derivations. We introduce symbolic density estimation (SDE), an unsupervised framework that automatically recovers closed-form probability mass functions by composing elementary analytic operations within a structured search space. Our method integrates domain-specific structural priors with evolutionary search and a validity-aware inference stage, and it extends to richer distribution families such as zero inflation and finite mixtures. To support systematic evaluation and future research, we contribute a benchmark dataset spanning a broad collection of commonly used discrete distributions. The proposed algorithm recovers all benchmark families with accurate parameter estimates. A real data application shows that it identifies concise and interpretable mixture models that improve goodness-of-fit over standard models.
Symbolic regression discovers explicit, interpretable equations without assuming a functional form in advance. A Bayesian approach strengthens this through probability distributions over candidate expressions, thus quantifying uncertainty in the presence of noisy and limited data. Deep Symbolic Regression (DSR) uses a neural network to generate symbolic expressions, but it is designed to identify a single best-fitting expression rather than infer a posterior distribution over models. We introduce Deep Variational Inference Symbolic Regression (DVISR), a variational Bayesian extension of DSR. DVISR replaces the original reward with the integrand of the evidence lower bound. It also extends the network architecture to output distributions over constants within expressions, enabling posterior inference over both expression trees and their associated constants. We show that DVISR can recover the true posterior in simple settings, both with and without constant tokens, and we examine how its performance changes as the size of the expression space increases. These results position DVISR as a step toward scalable Bayesian symbolic regression with uncertainty over full symbolic models.
James Butterworth, Gevik Grigorian, Alejandro DiazDelaO
Clinical Operational Research Unit University College London London, United Kingdom · Department of Computer Science University College London London, United Kingdom
Symbolic regression is the problem of finding an algebraic expression describing a stochastic dependence of a target variable on a set of inputs. Unlike forms of regression that fit parameters assuming a fixed model structure, symbolic regression is a search problem over the space of expressions, represented, for example, as abstract syntax trees using a library of operators. Symbolic regression is typically used in settings with limited, noisy data in the natural sciences. However, searching for a single best-fitting expression fails to capture the epistemic uncertainty about the expression, which motivates a Bayesian perspective that enables uncertainty quantification and specification of natural priors to constrain the search space. In this work, we propose ERRLESS (Entropy-Regularized Reinforcement Learning for Expression Structure Sampling), a scalable approach for sampling from the posterior distribution over expressions given data using maximum-entropy reinforcement learning. ERRLESS learns a neural policy that constructs expressions sequentially by building up their abstract syntax trees. At convergence, the policy samples expressions from the posterior. At test time, expressions can be sampled by rollouts of this policy. We demonstrate that ERRLESS achieves competitive results on the Feynman benchmark while producing short and interpretable expressions. Additionally, we demonstrate that the mean of the posterior predictive approximated by ERRLESS achieves a high coefficient of determination (R2) compared to an SMC baseline, highlighting the benefits of the Bayesian perspective in symbolic regression.
Oussama Boussif, Mohammed Mahfoud, Younesse Kaddar +6
Mila – Québec AI Institute · Université de Montréal · Independent +5