Symbolic Density Estimators for Unnormalized Distributions

Vikas Kanaujia; Riyansha Singh; Shashank Sharma; Vipul Arora

cs.LGOct 7, 2026

Symbolic Density Estimators for Unnormalized Distributions

Authors: Vikas Kanaujia, Riyansha Singh, Shashank Sharma, Vipul Arora

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

Abstract

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(\mathbf{x})$ can be factorized with a $2\times 2$ patch. (b) While inference with SR, different patch sizes provide different inductive biases. (c) Illustration for redundancy in representation.

\theta\leftarrow\theta-\eta\nabla_{\theta}\mathcal{L}_{M}(\theta)

Algorithm 1 Symbolic Function Estimation for Unnormalized Distributions

		MSE( $\downarrow$ )	$R^{2}(\uparrow)$	RCE( $\downarrow$ )	ST( $\downarrow$ )	P( $\uparrow$ )	R( $\uparrow$ )	F1( $\uparrow$ )	Estimated $\hat{h}(x_{p(l)})$
$\Sigma=\begin{bmatrix}1.0&0.0\\ 0.0&1.0\end{bmatrix}$	MeSSY	0.0114	0.99	0.0433	1	0.80	1.00	0.89	$0.48x_{1}^{2}-{\color[rgb]{0.5,0.5,0.5}0.02x_{1}x_{2}}+0.50x_{2}^{2}+1.39x_{1}-1.41x_{2}$
	NF + PySR	0.0001	0.99	0.0017	0	1.00	1.00	1.00	$0.5x_{1}^{2}+0.5x_{2}^{2}+1.5x_{1}-1.49x_{2}$
	NF + EQL	0.0011	0.99	0.0033	0	1.00	1.00	1.00	$0.5x_{1}^{2}+0.5x_{2}^{2}+1.52x_{1}-1.5x_{2}$
	SM + PySR	0.0013	0.99	0.0150	0	1.00	1.00	1.00	$0.51x_{1}^{2}+0.49x_{2}^{2}+1.53x_{1}-1.5x_{2}$
	SM + EQL	0.0014	0.99	0.0230	0	1.00	1.00	1.00	$0.49x_{1}^{2}+1.47x_{1}+0.52x_{2}^{2}-1.52x_{2}$
$\Sigma=\begin{bmatrix}1.0&0.5\\ 0.5&1.0\end{bmatrix}$	MeSSY	0.0493	0.95	0.1199	0	1.00	1.00	1.00	$0.72x_{1}^{2}-0.85x_{1}x_{2}+0.63x_{2}^{2}+3.38x_{1}-3.21x_{2}$

Table 1: Results and estimated expressions

\hat{h}([x_{1},x_{2}];\lambda)

for a 2-D Gaussian distribution using various methods. The mean vector is fixed at

\mu=[-1.5,1.5]^{T}

, with two different covariance matrices,

\Sigma_{1}

and

\Sigma_{2}

, as detailed in the table. The target functions are: (A)

0.5x_{1}^{2}+0.5x_{2}^{2}+1.5x_{1}-1.5x_{2}

, and (B)

0.67x_{1}^{2}-0.67x_{1}x_{2}+0.67x_{2}^{2}+3.0x_{1}-3.0x_{2}

. Incorrect terms are shown in gray. NF = Flow model, SM = Score matching.

$d$	Model	MSE( $\downarrow$ )	$R^{2}(\uparrow)$	RCE( $\downarrow$ )	ST( $\downarrow$ )	P( $\uparrow$ )	R( $\uparrow$ )	F1( $\uparrow$ )	Estimated Equation $\hat{h}(x_{p(l)})$
2	MeSSY	0.02	0.98	0.2396	2	0.67	1.00	0.80	$0.96x_{1}^{4}+{\color[rgb]{0.5,0.5,0.5}0.1x_{1}^{3}}-5.77x_{1}^{2}-0.83x_{1}+0.61x_{2}^{2}-{\color[rgb]{0.5,0.5,0.5}0.23x_{2}}$
	NF + PySR	0.00	0.99	0.0008	0	1.00	1.00	1.00	$x_{1}^{4}-6.02x_{1}^{2}-0.5x_{1}+0.5x_{2}^{2}$
	NF + EQL	0.00	0.99	0.0196	0	1.00	1.00	1.00	$0.99x_{1}^{4}-5.95x_{1}^{2}-0.49x_{1}+0.52x_{2}^{2}$
	SM + PySR	0.00	0.99	0.0142	0	1.00	1.00	1.00	$0.98x_{1}^{4}-5.9x_{1}^{2}-0.51x_{1}+0.5x_{2}^{2}$
	SM + EQL	0.00	0.99	0.0617	0	1.00	1.00	1.00	$0.99x_{1}^{4}-5.9x_{1}^{2}-0.59x_{1}+0.48x_{2}^{2}$
4	NF + PySR	0.00	0.97	0.0296	0	1.00	1.00	1.00	$1.02x_{1}^{4}-6.11x_{1}^{2}-0.52x_{1}+0.52x_{2}^{2}$

Table 2: Results for Many Well distribution with various dimensions,

d=2,4,8,16,32,64

. The true function is

h(x_{p_{l}})=x_{1}^{4}-6x_{1}^{2}-0.5x_{1}+0.5x_{2}^{2}

. Estimated terms that do not match the true function are shown in gray. NF = Flow model, SM = Score matching.

Model	MSE( $\downarrow$ )	$R^{2}(\uparrow)$	RCE( $\downarrow$ )	ST( $\downarrow$ )	P( $\uparrow$ )	R( $\uparrow$ )	F1( $\uparrow$ )	Estimated $\hat{h}(x_{p(l)})$
MeSSY	-	-	-	-	-	-	-	-
NF + PySR	4.83	0.92	0.1127	0	1.00	1.00	1.00	$-0.63\cos{\left(0.9x_{0}-0.9x_{1}\right)}-0.63\cos{\left(x_{0}-x_{2}\right)}$
DSM + PySR	42.22	0.30	0.6200	2	0.50	0.50	0.50	$-0.27\cos(x_{0}-x_{1})-{\color[rgb]{0.5,0.5,0.5}0.04x_{0}x_{3}}+{\color[rgb]{0.5,0.5,0.5}0.04x_{3}^{2}}$
NF + EQL	6.74	0.87	0.0141	4	0.67	1.00	0.80	$-0.72\cos(x_{0}-x_{1})-0.72\cos(x_{0}-x_{2})+{\color[rgb]{0.5,0.5,0.5}0.11\cos(x_{0}-x_{3})}+{\color[rgb]{0.5,0.5,0.5}0.19\cos^{2}(x_{3})}-{\color[rgb]{0.5,0.5,0.5}0.28\cos(x_{0})+{\color[rgb]{0.5,0.5,0.5}0.32\cos(x_{3})}}$
DSM + EQL	47.72	0.21	N/A	2	0.33	0.50	0.40	$-{\color[rgb]{0.5,0.5,0.5}0.46\cos(0.34x_{0}-0.59x_{2}+0.30x_{3}}+{\color[rgb]{0.5,0.5,0.5}0.64\cos(0.48x_{0}-0.61x_{2}-0.33x_{3}))}$

Table 3: Results for XY Model. The true function is

h(x_{p(l)})=-0.71\cos(x_{0}-x_{1})-0.71\cos(x_{0}-x_{2})

. Estimated terms that do not match the true function are shown in gray. NF = Flow model, DSM = Denoising score model.

	Model	MSE( $\downarrow$ )	$R^{2}(\uparrow)$	RCE( $\downarrow$ )	ST( $\downarrow$ )	P( $\uparrow$ )	R( $\uparrow$ )	F1( $\uparrow$ )	Estimated $\hat{h}(x_{p(l)})$
(A)	MeSSY	-	-	-	-	-	-	-	-
	NF + PySR	0.02	0.99	0.0175	0	1.00	1.00	1.00	$3.91x_{0}^{4}-1.97x_{0}x_{1}-1.97x_{0}x_{2}$
	NF + EQL	0.29	0.99	0.1583	4	0.43	1.00	0.60	$3.42x_{0}^{4}-2.33x_{0}x_{1}-2.33x_{0}x_{2}+{\color[rgb]{0.5,0.5,0.5}0.29x_{0}x_{1}^{3}}+{\color[rgb]{0.5,0.5,0.5}0.40x_{1}^{3}x_{3}}-{\color[rgb]{0.5,0.5,0.5}0.13x_{0}x_{3}}+{\color[rgb]{0.5,0.5,0.5}0.49x_{1}x_{2}}$
	DSM + PySR	0.01	0.99	0.0108	0	1.00	1.00	1.00	$3.95x_{0}^{4}-1.98x_{0}x_{1}-1.98x_{0}x_{2}$
	DSM + EQL	0.06	0.99	0.0358	1	0.75	1.00	0.86	$3.65x_{0}^{4}-2.01x_{0}x_{1}-2.03x_{0}x_{2}+{\color[rgb]{0.5,0.5,0.5}0.12x_{0}^{2}}$
(B)	MeSSY	-	-	-	-	-	-	-	-

Table 4: Results for

\phi^{4}

theory with various

\lambda_{1},\lambda_{2}

settings. The true expressions are: (A) for

\lambda_{1}=4,\lambda_{2}=-4

h(x_{p(l)})=4x_{0}^{4}-2x_{0}x_{1}-2x_{0}x_{2}

, and (B) for

\lambda_{1}=4,\lambda_{2}=1

h(x_{p(l)})=4x_{0}^{4}+5x_{0}^{2}-2x_{0}x_{1}-2x_{0}x_{2}

. Incorrect terms are highlighted in gray. NF = Flow model, DSM = Denoising score model.

	$\|p\|$	MSE( $\downarrow$ )	$R^{2}(\uparrow)$	RCE( $\downarrow$ )	ST( $\downarrow$ )	P( $\uparrow$ )	R( $\uparrow$ )	F1( $\uparrow$ )	Estimated $\hat{h}(x_{p(l)})$
PySR	$2\times 2$	0.003	0.99	0.0241	0	1.00	1.00	1.00	$4.01x_{0}^{4}+5.17x_{0}^{2}-2.06x_{0}x_{1}-2.06x_{0}x_{2}$
	$3\times 3$	0.006	0.99	0.0199	0	1.00	1.00	1.00	$4.11x_{0}^{4}+5.06x_{0}^{2}-2.04x_{0}x_{1}-2.04x_{0}x_{2}$
	$4\times 4$	0.008	0.99	0.0158	0	1.00	1.00	1.00	$3.90x_{0}^{4}+5.14x_{0}^{2}-1.99x_{0}x_{1}-1.99x_{0}x_{2}$
EQL	$2\times 2$	0.07	0.99	0.0516	4	0.50	1.00	0.67	$3.69x_{0}^{4}+5.27x_{0}^{2}-2.01x_{0}x_{2}-2.14x_{0}x_{1}-{\color[rgb]{0.5,0.5,0.5}0.68x_{1}x_{2}^{3}}-{\color[rgb]{0.5,0.5,0.5}0.27x_{2}^{3}x_{3}}+{\color[rgb]{0.5,0.5,0.5}0.14x_{0}x_{3}}+{\color[rgb]{0.5,0.5,0.5}0.45x_{1}x_{2}}$
	$3\times 3$	0.98	0.98	0.0330	4	0.50	1.00	0.67	$3.86x_{0}^{4}+5.21x_{0}^{2}-2.1x_{0}x_{1}-1.99x_{0}x_{2}+{\color[rgb]{0.5,0.5,0.5}0.42x_{1}x_{2}}-{\color[rgb]{0.5,0.5,0.5}1.01x_{0}^{3}x_{2}}-{\color[rgb]{0.5,0.5,0.5}0.14x_{0}^{3}x_{1}}-{\color[rgb]{0.5,0.5,0.5}0.13x_{0}^{3}x_{5}}$
	$4\times 4$	0.11	0.99	0.0409	7	0.36	1.00	0.53	$3.61x_{0}^{4}+4.97x_{0}^{2}-2.09x_{0}x_{1}-2.03x_{0}x_{2}+{\color[rgb]{0.5,0.5,0.5}0.11x_{0}^{3}x_{3}}-{\color[rgb]{0.5,0.5,0.5}0.12x_{0}^{3}x_{8}}-{\color[rgb]{0.5,0.5,0.5}0.1x_{0}x_{9}}+{\color[rgb]{0.5,0.5,0.5}0.28x_{1}x_{4}}+{\color[rgb]{0.5,0.5,0.5}0.27x_{0}x_{5}}+{\color[rgb]{0.5,0.5,0.5}0.14x_{0}x_{2}}+{\color[rgb]{0.5,0.5,0.5}0.13x_{0}x_{8}}$

Table 5: Results for estimating

\hat{h}(x_{p(l)})

for

\phi^{4}

theory with different patch sizes

|p|

for flow model. True

h(x_{p(l)})=4x_{0}^{4}+5x_{0}^{2}-2x_{0}x_{1}-2x_{0}x_{2}

. Incorrect terms are shown in gray.

Figure 3: Renormalization results for non-perturbative region, $\lambda_{1}=4$ and $\lambda_{2}=1$ .

	32x32		16x16		8x8
Observable	% Overlap( $\uparrow$ )	EMD ( $\downarrow$ )	% Overlap( $\uparrow$ )	EMD ( $\downarrow$ )	% Overlap( $\uparrow$ )	EMD ( $\downarrow$ )
Magnetization, M	98.116	$2.28\times 10^{-4}$	98.845	$9.09\times 10^{-5}$	98.759	$8.79\times 10^{-5}$
Absolute Magnetization, \|M\|	98.197	$2.06\times 10^{-4}$	98.794	$8.69\times 10^{-5}$	98.749	$8.30\times 10^{-5}$
Magnetic Susceptibility	98.244	$3.81\times 10^{-3}$	99.187	$4.47\times 10^{-4}$	99.187	$1.13\times 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 ( $\sigma$ )	Mean	Std Dev ( $\sigma$ )	Mean	Std Dev ( $\sigma$ )
$32\times 32$	Coarsened Data	$-1.82\times 10^{-6}$	0.0066	0.0053	0.0039	0.0441	0.0609
$32\times 32$	Generated Data	$9.09\times 10^{-5}$	0.0068	0.0055	0.0041	0.0481	0.0681
$16\times 16$	Coarsened Data	$-1.82\times 10^{-6}$	0.0066	0.0053	0.0039	0.0110	0.0152
$16\times 16$	Generated Data	$-3.36\times 10^{-5}$	0.0067	0.0053	0.0040	0.0115	0.0161
$8\times 8$	Coarsened Data	$-1.82\times 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_{\theta})$	MSE $(H,\hat{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; $\lambda=1.4^{-1}$	9.79	4.83	6.74
$\phi^{4};\lambda_{1}=4,\lambda_{2}=-4$	1.43	0.02	0.29

Table 8: Comparison of MSE

(H,H_{\theta})

, where

H_{\theta}

is computed by flow model, and MSE

(H,\hat{H})

, where

\hat{H}

is estimated by SR methods (PySR and EQL) using flow model outputs.

Appendix figures & tables11 assets

Supplementary material from the paper’s appendix.

Appendix

Experiment	Operator Library	Layers	Patch Size	Range of Interaction	Exact Terms Present
Multivariate Gaussian	Const(2), Id(2), Cos(1), Square(2), Exp(2), Mul(2), Pow4(2)	2/3	–	–	Yes
Many-Well	Const(2), Id(2), Square(3), Exp(2), Mul(1), Pow4(2)	2/3	$2\times 1$	$2\times 1$	Yes
$\phi^{4}$	Const(1), Id(2), Square(2), Exp(1), Mul(2), Pow4(2)	2	$2\times 2$ , $3\times 3$ , $4\times 4$	$2\times 2$	Yes
XY	Const(1), Id(1), Cos(1), Sin(1), Square(1), Mul(2)	2	$2\times 2$ , $3\times 3$	$2\times 2$	Yes

Appendix

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	$\{+,-,*,\mathrm{pow2}\}$	30	30	$2\times 1$	$2\times 1$	Yes
$\phi^{4}$	$\{+,-,*,\mathrm{pow2},\mathrm{pow3},\mathrm{pow4}\}$	30	30	$2\times 2$ , $3\times 3$ , $4\times 4$	$2\times 2$	Yes
XY	$\{+,-,*,\sin,\cos\}$	30	30	$2\times 2$ , $3\times 3$	$2\times 2$	Yes

Appendix

Table 10: PySR configurations and inductive biases used across benchmarks.

No. of Samples	MSE( $\downarrow$ )	$R^{2}(\uparrow)$	RCE( $\downarrow$ )	P( $\uparrow$ )	R( $\uparrow$ )	F1( $\uparrow$ )	Estimated Equation
10K	0.47	0.994	0.09	1	1	1	$3.5x_{0}^{4}-1.85x_{0}x_{1}-1.85x_{0}x_{2}$
50K	0.07	0.999	0.03	1	1	1	$4.2x_{0}^{4}-2.04x_{0}x_{1}-2.04x_{0}x_{2}$
100K	0.01	0.999	0.02	1	1	1	$3.9x_{0}^{4}-1.95x_{0}x_{1}-2.00x_{0}x_{2}$

Appendix

Table 11: Effect of sample size on model performance.

	MSE( $\downarrow$ )	$R^{2}(\uparrow)$	RCE( $\downarrow$ )	ST( $\downarrow$ )	P( $\uparrow$ )	R( $\uparrow$ )	F1( $\uparrow$ )	Estimated Equation $\hat{h}(x_{p(l)})$
1	0.00	0.999	0.024	0	1.00	1.00	1.00	$4.01x_{0}^{4}+5.17x_{0}^{2}-2.06x_{0}x_{1}-2.06x_{0}x_{2}$
2	0.02	0.999	0.063	0	1.00	1.00	1.00	$4.76x_{0}^{4}+4.76x_{0}^{2}-1.97x_{0}x_{1}-2.00x_{0}x_{2}$
3	0.01	0.999	0.136	0	1.00	1.00	1.00	$3.72x_{0}^{4}+7.27x_{0}^{2}-2.02x_{0}x_{1}-2.02x_{0}x_{2}$
4	0.04	0.998	0.075	1	0.75	0.75	0.75	$5.93x_{0}^{2}-2.04x_{0}x_{1}-2.04x_{0}x_{2}+{\color[rgb]{0.5,0.5,0.5}4.21x_{0}^{6}}$
5	0.11	0.996	0.042	0	1.00	1.00	1.00	$3.82x_{0}^{4}+5.01x_{0}^{2}-2.12x_{0}x_{1}-2.12x_{0}x_{2}$
6	0.05	0.998	0.016	1	0.75	0.75	0.75	$5.01x_{0}^{2}-2.01x_{0}x_{1}-2.08x_{0}x_{2}+{\color[rgb]{0.5,0.5,0.5}3.65x_{0}^{6}}$

Appendix

Table 12: Symbolic estimation results for the

\phi^{4}

theory across 20 different random seed settings using NF + PySR. The ground-truth local Hamiltonian for

\lambda_{1}=4

and

\lambda_{2}=1

h(x_{p(l)})=4x_{0}^{4}+5x_{0}^{2}-2x_{0}x_{1}-2x_{0}x_{2}

. Incorrect terms are highlighted in gray. NF denotes the flow-based model.

	MSE( $\downarrow$ )	$R^{2}(\uparrow)$	RCE( $\downarrow$ )	ST( $\downarrow$ )	P( $\uparrow$ )	R( $\uparrow$ )	F1( $\uparrow$ )	Estimated Equation $\hat{h}(x_{p(l)})$
1	0.07	0.99	0.058	4	0.50	1.00	0.67	$3.69x_{0}^{4}+5.27x_{0}^{2}-2.01x_{0}x_{2}-2.14x_{0}x_{1}-{\color[rgb]{0.5,0.5,0.5}0.68x_{1}x_{2}^{3}}-{\color[rgb]{0.5,0.5,0.5}0.27x_{2}^{3}x_{3}}+{\color[rgb]{0.5,0.5,0.5}0.14x_{0}x_{3}}+{\color[rgb]{0.5,0.5,0.5}0.45x_{1}x_{2}}$
2	0.054	0.998	0.067	11	0.27	1.00	0.42	$3.46x_{0}^{4}+4.78x_{0}^{2}-2.14x_{0}x_{1}-2.01x_{0}x_{2}-{\color[rgb]{0.5,0.5,0.5}0.13x_{0}x_{1}^{3}-0.32x_{0}x_{2}^{3}+\cdots\text{Other terms}}$
3	0.048	0.998	0.089	17	0.19	1.00	0.32	$3.32x_{0}^{4}+5.08x_{0}^{2}-2.07x_{0}x_{1}-2.28x_{0}x_{2}+{\color[rgb]{0.5,0.5,0.5}0.06x_{0}x_{3}^{3}+0.23x_{0}x_{1}^{3}+\cdots\text{Other terms}}$
4	0.054	0.098	0.065	9	0.31	1.00	0.47	$3.46x_{0}^{4}+4.78x_{0}^{2}-2.14x_{0}x_{1}-2.04x_{0}x_{2}-{\color[rgb]{0.5,0.5,0.5}0.32x_{0}x_{2}^{3}-0.13x_{0}x_{1}^{3}+\cdots\text{Other terms}}$
5	0.041	0.998	0.066	24	0.14	1.00	0.25	$3.61x_{0}^{4}+5.06x_{0}^{2}-2.26x_{0}x_{1}-2.05x_{0}x_{2}-{\color[rgb]{0.5,0.5,0.5}0.16x_{1}^{2}x_{2}^{2}+0.28x_{0}^{2}x_{1}^{2}+\cdots\text{Other terms}}$

Appendix

Table 13: Symbolic estimation results for the

\phi^{4}

theory across 5 different random seed settings using NF + EQL. The true expressions is

h(x_{p(l)})=4x_{0}^{4}+5x_{0}^{2}-2x_{0}x_{1}-2x_{0}x_{2}

. Incorrect terms are highlighted in gray.

	MSE( $\downarrow$ )	$R^{2}(\uparrow)$	RCE( $\downarrow$ )	ST( $\downarrow$ )	P( $\uparrow$ )	R( $\uparrow$ )	F1( $\uparrow$ )	Estimated Equation $\hat{h}(x_{p(l)})$
1	0.72	0.985	0.279	0	1.00	1.00	1.00	$1.58x_{1}^{4}-7.27x_{1}^{2}-0.65x_{1}+0.51x_{2}^{2}$
2	0.46	0.990	0.126	0	1.00	1.00	1.00	$1.37x_{1}^{4}-6.64x_{1}^{2}-0.51x_{1}+0.50x_{2}^{2}$
3	0.32	0.993	0.242	1	0.80	1.00	0.89	$1.68x_{1}^{4}-7.46x_{1}^{2}-0.52x_{1}+0.50x_{2}^{2}$
4	0.28	0.994	0.160	0	1.00	1.00	1.00	$1.46x_{1}^{4}-6.96x_{1}^{2}-0.50x_{1}+0.51x_{2}^{2}$
5	0.40	0.991	0.179	0	1.00	1.00	1.00	$1.44x_{1}^{2}-6.92x_{1}^{2}-0.45x_{1}+0.51x_{2}^{2}$
6	0.42	0.991	0.189	0	1.00	1.00	1.00	$1.49x_{1}^{4}-7.02x_{1}^{2}-0.55x_{1}+0.50x_{2}^{2}$

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(x_{p_{l}})=x_{1}^{4}-6x_{1}^{2}-0.5x_{1}+0.5x_{2}^{2}

. Incorrect terms are highlighted in gray.

	MSE( $\downarrow$ )	$R^{2}(\uparrow)$	RCE( $\downarrow$ )	ST( $\downarrow$ )	P( $\uparrow$ )	R( $\uparrow$ )	F1( $\uparrow$ )	Estimated Equation $\hat{h}(x_{p(l)})$
1	0.57	0.997	0.583	3	0.57	1.00	0.73	$1.03x_{1}^{4}-6.01x_{1}^{2}-1.64x_{1}+0.49x_{2}^{2}-{\color[rgb]{0.5,0.5,0.5}0.22x_{1}^{5}}+{\color[rgb]{0.5,0.5,0.5}0.64x_{1}^{3}}+{\color[rgb]{0.5,0.5,0.5}0.23x_{1}x_{2}}$
2	0.92	0.971	1.452	2	0.67	1.00	0.80	$0.89x_{1}^{4}-5.65x_{1}^{2}+2.30x_{1}+0.48x_{2}^{2}+{\color[rgb]{0.5,0.5,0.5}0.14x_{1}^{5}-1.09x_{1}^{3}}$
3	1.78	0.953	0.868	4	0.50	1.00	0.67	$0.51x_{1}^{4}-5.14x_{1}^{2}+0.89x_{1}+0.47x_{2}^{2}+{\color[rgb]{0.5,0.5,0.5}0.16x_{1}^{6}-0.12x_{1}^{5}-0.24x_{1}^{3}-0.14x_{1}x_{2}}$
4	0.57	0.976	1.517	2	0.67	1.00	0.80	$0.81x_{1}^{4}-5.52x_{1}^{2}+2.38x_{1}+0.48x_{2}^{2}+{\color[rgb]{0.5,0.5,0.5}0.16x_{1}^{5}-1.16x_{1}^{3}}$
5	2.01	0.949	1.570	3	0.57	1.00	0.73	$0.73x_{1}^{4}-5.46x_{1}^{2}+2.44x_{1}+0.48x_{2}^{2}+{\color[rgb]{0.5,0.5,0.5}0.19x_{1}^{5}-1.23x_{1}^{3}+0.13x_{1}x_{2}}$
6	2.58	0.928	1.205	2	0.67	1.00	0.80	$1.33x_{1}^{4}-6.49x_{1}^{2}+1.69x_{1}+0.48x_{2}+{\color[rgb]{0.5,0.5,0.5}0.13x_{1}^{5}-0.89x_{1}^{3}}$

Appendix

Table 15: Symbolic estimation results for the MW-64 benchmark across 8 random seed settings using SM + EQL. The true local Hamiltonian is

h(x_{p_{l}})=x_{1}^{4}-6x_{1}^{2}-0.5x_{1}+0.5x_{2}^{2}

. Incorrect terms are highlighted in gray.

Method	MSE( $\downarrow$ )	$R^{2}(\uparrow)$	RCE( $\downarrow$ )	ST( $\downarrow$ )	P( $\uparrow$ )	R( $\uparrow$ )	F1( $\uparrow$ )	Estimated Equation $\hat{h}(x_{p(l)})$
NF + PySR	0.7711	0.9704	0.1428	0	1.00	1.00	1.00	$4.73x_{12}^{4}+4.73x_{12}^{2}-0.83x_{2}x_{12}-0.84x_{10}x_{12}-0.85x_{12}x_{14}-0.86x_{12}x_{22}$
DSM + PySR	0.8372	0.9680	0.0570	0	1.00	1.00	1.00	$3.40x_{12}^{4}+5.16x_{12}^{2}-0.96x_{2}x_{12}-0.95x_{10}x_{12}-0.97x_{12}x_{14}-0.96x_{12}x_{22}$
NF + EQL	13.8171	0.4225	1.5385	1467	0.003	1.00	0.005	$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}}$
DSM + EQL	23.8519	0.0031	0.9181	4	0.50	0.67	0.57	$0.03x_{12}^{4}-0.12x_{2}x_{12}-0.10x_{10}x_{12}-0.10x_{12}x_{14}-{\color[rgb]{0.5,0.5,0.5}0.16x_{7}x_{12}+0.03x_{11}x_{12}-0.05x_{12}x_{13}+0.06x_{12}x_{17}}$

Appendix

Table 16: Results for a synthetic long-range interaction system. The true expression is

h(x_{p(l)})=4x_{12}^{4}+5x_{12}^{2}-x_{2}x_{12}-x_{10}x_{12}-x_{12}x_{14}-x_{12}x_{22}

Incorrect terms are highlighted in gray. NF = Flow model, DSM = Denoising score model.

Method	MSE( $\downarrow$ )	$R^{2}$ ( $\uparrow$ )	Estimated $\hat{h}(x_{p(l)})$
NF + PySR	0.0085	0.99	$3.89x_{0}^{4}+5.15x_{0}^{2}-2.00x_{0}x_{1}-2.00x_{0}x_{2}$
NF + EQL	2.2872	0.90	$5.10x_{0}^{2}-2.16x_{0}x_{1}-1.75x_{0}x_{2}+{\color[rgb]{0.5,0.5,0.5}2.62x_{0}x_{3}}+{\color[rgb]{0.5,0.5,0.5}0.25x_{1}x_{2}}+\cdots+{\color[rgb]{0.5,0.5,0.5}\text{higher-order polynomial terms}}+\cdots+{\color[rgb]{0.5,0.5,0.5}\text{nested }\sin,\cos\text{ terms}}$

Appendix

Table 17: Symbolic estimation on the scalar

\phi^{4}

benchmark under primitive library misspecification. Explicit polynomial primitives are removed from the symbolic library, leaving only arithmetic and trigonometric operators.

Method	MSE( $\downarrow$ )	$R^{2}$ ( $\uparrow$ )	Estimated $\hat{h}(x_{p(l)})$
NF + PySR	8.84	0.81	$-0.03x_{0}^{2}x_{1}^{2}-0.03x_{0}^{2}x_{2}^{2}+0.03x_{0}x_{1}^{3}+0.02x_{0}x_{2}^{3}+0.02x_{0}^{3}x_{2}+0.15x_{0}^{2}+0.56x_{1}^{2}-0.43x_{0}x_{1}-0.30x_{0}x_{2}-{\color[rgb]{0.5,0.5,0.5}0.07x_{1}^{3}-0.07x_{0}x_{1}^{2}+0.15x_{0}^{2}x_{1}}+{\color[rgb]{0.5,0.5,0.5}0.15x_{1}}$
NF + EQL	7.13	0.85	$0.46x_{0}^{2}+0.14x_{1}^{2}+0.17x_{2}^{2}-0.37x_{0}x_{1}-0.24x_{0}x_{2}+{\color[rgb]{0.5,0.5,0.5}0.14x_{1}+0.11x_{0}}$

Appendix

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	$\mathrm{MSE}(H,H_{\theta})$	$\mathrm{MSE}(H,\hat{H})$
$\phi^{4};\lambda_{1}=4,\lambda_{2}=1$	NF + PySR	0.22	0.01
$\phi^{4};\lambda_{1}=4,\lambda_{2}=1$	NF + EQL	0.22	2.29
XY; $\lambda=1.4^{-1}$	NF + PySR	9.79	8.84
XY; $\lambda=1.4^{-1}$	NF + EQL	9.79	7.13

Appendix

Table 19: Error propagation analysis under primitive library misspecification.

\mathrm{MSE}(H,H_{\theta})

denotes the error of the learned flow model relative to the true Hamiltonian, while

\mathrm{MSE}(H,\hat{H})

denotes the error of the estimated symbolic approximation.

Explore similar work

Symbolic Density Estimation for Discrete Distributions
May 20, 2026Ziwen Liu, Meng LiDensity Estimation
Rice University Houston, TX
Deep Variational Inference Symbolic Regression
May 1, 2026James Butterworth, Gevik Grigorian, Alejandro DiazDelaOUncertainty Quantification Variational Inference
Clinical Operational Research Unit University College London London, United Kingdom · Department of Computer Science University College London London, United Kingdom
Bayesian Symbolic Regression with Entropic Reinforcement Learning
Aug 10, 2026Oussama Boussif, Mohammed Mahfoud, Younesse Kaddar +6Bayesian Inference Symbolic Regression
Mila – Québec AI Institute · Université de Montréal · Independent +5

May 20, 2026, cs.LGJ/K move · Enter open · S save

Symbolic Density Estimation for Discrete Distributions

Ziwen Liu, Meng Li

Rice University Houston, TX

Density Estimation

Open paper AskPDF arXiv

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.

Abstract

Explore similar work

Symbolic Density Estimation for Discrete Distributions

Deep Variational Inference Symbolic Regression

Bayesian Symbolic Regression with Entropic Reinforcement Learning