Product Monte Carlo: Expectation Under a Product of Densities

Pietro Carlotti; Sifan Liu; Yuling Yao

stat.MEOct 8, 2026

Product Monte Carlo: Expectation Under a Product of Densities

Authors: Pietro Carlotti, Sifan Liu, Yuling Yao

Organizations: University of Texas at Austin, Department of Statistics and Data Sciences · Duke University, Department of Statistical Science

Abstract

Aggregating independently computed sub-sample inferences is a well-known challenge in large-scale Bayesian computation, federated learning, and opinion pooling. The core problem is to approximate a distribution whose density is proportional to a product of individual densities. We propose Product Monte Carlo (PMC), an algorithm for computing arbitrary expectations under a product of densities using only Monte Carlo draws from individual densities. We turn the intractable product expectation into a regression problem by expressing it through a family of computable conditional expectations. Besides the point estimate, our fully Bayesian regression formulation provides coherent uncertainty quantification and principled diagnostics. We establish theoretical guarantees on consistency and characterize convergence rates for the PMC estimator. Across simulated and real-data experiments varying the number, dimensionality, and heterogeneity of the component distributions, PMC achieves higher estimation accuracy, better-calibrated uncertainty coverage, and improved scalability at a lower computational cost than existing methods. In addition, PMC requires no density evaluation, making it applicable to simulation-based inference, as we demonstrate in an application with protein images.

Figures & tables

Panel of Figure 1 — Figure 1: PMC aggregates samples from individual components to estimate any expectation under the product distribution. No density evaluation of $\pi_{k}$ is needed.

Figure 2: For $\pi_{1}=\mathcal{N}(1,1)$ and $\pi_{2}=\mathcal{N}(2,1)$ , the product distribution $\pi_{\ast}=\mathcal{N}(1.5,0.5)\propto\pi_{1}\pi_{2}$ matches the conditional distribution on the diagonal slice.

Method	Distribution-free	Evaluation-free	UQ
Parametric approx.		$\checkmark$
Weighting	$\checkmark$
Fusion	$\checkmark$
PMC	$\checkmark$	$\checkmark$	$\checkmark$

Table 1: Comparison of available methods and PMC.

Figure 3: Using posterior predictive check in PMC. Left: true, well-specified, and misspecified regressions of

h(C)

D

. Right: posterior predictive distributions of

T

, a statistic of simulated data, under the two fitted models. Orange lines mark the observed statistics.

Figure 4: MSE and runtime versus

S

for the dimension-scaling experiment for the standardized sum of coordinates.

Method	Estimation error	Merge (min)	Total (min)
Full MCMC	0.0000	–	0.8152
Subsample	–	–	0.9796
GBF	0.0184	0.2047	1.1843
PMC-OLS	0.0110	0.0002	0.9799
PMC-MLP	0.0144	0.1069	1.0866
PMC-GP	0.0258	0.1009	1.0806

Table 2: Application in a privacy-restricted clinical trial.

Method	Estimation error	Merge (min)	Total (min)
Full MCMC	0.0000	–	308.1124
Subsample	–	–	35.2890
PMC-OLS	0.0020	0.0122	35.3012
PMC-MLP	0.0001	0.1092	35.3982
PMC-GP	0.0001	6.7192	42.0082

Table 3: Application in large-scale inference. PMC’s minimal overhead keeps total cost a fraction of full-data MCMC.

Figure 7: Application to SBI. Left: Image-level subposteriors and the combined posterior of the latent parameter

z

. PMC-MLP and PMC-GP estimates recover the true posterior mean. Right Table: PMC closely approximates the full-data posterior expectation with negligible merging cost relative to image-level subsampling SBI.

Appendix figures & tables58 assets

Supplementary material from the paper’s appendix.

Appendix

Figure 8: MSE and runtime versus

S

for

h(u)=u

with

K=2

univariate Gaussian components. Columns and rows vary

\mu_{2},\sigma_{2}\in\left\{1,2,3,4\right\}

, respectively. The first component is fixed with

\mu_{1}=0

and

\sigma_{1}=1

Figure 9: MSE and runtime versus

S

for

h(u)=u^{2}

with

K=2

univariate Gaussian components. Columns and rows vary

\mu_{2},\sigma_{2}\in\left\{1,2,3,4\right\}

, respectively. The first component is fixed with

\mu_{1}=0

and

\sigma_{1}=1

Figure 10: MSE and runtime versus

S

for

h(u)=u

with

K

homogeneous univariate Gaussian components. Columns vary

K\in\left\{2,3,5,10,20\right\}

, and rows show MSE and runtime, respectively.

Figure 11: MSE and runtime versus

S

for

h(u)=u^{2}

with

K

homogeneous univariate Gaussian components. Columns vary

K\in\left\{2,3,5,10,20\right\}

, and rows show MSE and runtime, respectively.

Figure 12: MSE and runtime versus

S

for

h(u)=u

with

K

heterogeneous univariate Gaussian components. Columns vary

K\in\left\{2,3,5,10,20\right\}

, and rows show MSE and runtime, respectively.

Figure 13: MSE and runtime versus

S

for

h(u)=u^{2}

with

K

heterogeneous univariate Gaussian components. Columns vary

K\in\left\{2,3,5,10,20\right\}

, and rows show MSE and runtime, respectively.

Figure 14: MSE and runtime versus

S

for

h(u)=u

with

K=2

bivariate Gaussian components. Columns vary

\mu\in\left\{1,2,3\right\}

. The first component is fixed with

\mu=0

Figure 15: MSE and runtime versus

S

for

h(u)=u\odot u

with

K=2

bivariate Gaussian components. Columns vary

\mu\in\left\{1,2,3\right\}

. The first component is fixed with

\mu=0

Figure 16: MSE and runtime versus

S

for

h(u)=u

with

K=2

bivariate Gaussian components. Columns vary

\sigma\in\left\{1,2,3,4\right\}

. The first component is fixed with

\sigma=1

Figure 17: MSE and runtime versus

S

for

h(u)=u\odot u

with

K=2

bivariate Gaussian components. Columns vary

\sigma\in\left\{1,2,3,4\right\}

. The first component is fixed with

\sigma=1

Figure 18: MSE and runtime versus

S

for

h(u)=u

with

K=2

bivariate Gaussian components. Columns vary

\sigma_{2}\in\left\{1,2,3,4\right\}

. The first component is fixed with

\sigma_{2}=1

Figure 19: MSE and runtime versus

S

for

h(u)=u\odot u

with

K=2

bivariate Gaussian components. Columns vary

\sigma_{2}\in\left\{1,2,3,4\right\}

. The first component is fixed with

\sigma_{2}=1

Figure 20: MSE and runtime versus

S

for

h(u)=u

with

K=2

bivariate Gaussian components. Columns vary

\rho\in\left\{-0.6,-0.3,0,0.3,0.6\right\}

. The first component is fixed with

\rho=0

Figure 21: MSE and runtime versus

S

for

h(u)=u\odot u

with

K=2

bivariate Gaussian components. Columns vary

\rho\in\left\{-0.6,-0.3,0,0.3,0.6\right\}

. The first component is fixed with

\rho=0

Figure 22: MSE and runtime versus

S

for

h(u)=u

with

K=2

multivariate Gaussian components. Columns vary

p\in\left\{2,3,5,10,20\right\}

, and rows show MSE and runtime, respectively.

Figure 23: MSE and runtime versus

S

for

h(u)=u\odot u

with

K=2

multivariate Gaussian components. Columns vary

p\in\left\{2,3,5,10,20\right\}

, and rows show MSE and runtime, respectively.

Figure 24: MSE and runtime versus

S

for

h(u)=u

with

K=2

symmetric bimodal Gaussian-mixture components.

Figure 25: MSE and runtime versus

S

for

h(u)=u^{2}

with

K=2

symmetric bimodal Gaussian-mixture components.

Figure 26: MSE and runtime versus

S

for

h(u)=u

with

K=2

Gaussian-mixture components with opposing mode weights.

Figure 27: MSE and runtime versus

S

for

h(u)=u^{2}

with

K=2

Gaussian-mixture components with opposing mode weights.

Figure 28: MSE and runtime versus

S

for

h(u)=u

with

K=2

unimodal Gaussian-mixture components with heterogeneous scales.

Figure 29: MSE and runtime versus

S

for

h(u)=u^{2}

with

K=2

unimodal Gaussian-mixture components with heterogeneous scales.

Figure 30: MSE and runtime versus

S

for

h(u)=u

with

K=3

Gaussian-mixture components.

Figure 31: MSE and runtime versus

S

for

h(u)=u^{2}

with

K=3

Gaussian-mixture components.

Figure 32: MSE and runtime versus

S

for

h(u)=u

with

K=2

heterogeneous Gaussian-mixture components.

Figure 33: MSE and runtime versus

S

for

h(u)=u^{2}

with

K=2

heterogeneous Gaussian-mixture components.

Figure 34: MSE and runtime versus

S

for

h(u)=u

with

K=2

homogeneous Beta components. Columns vary

a\in\left\{1,2,5\right\}

Figure 35: MSE and runtime versus

S

for

h(u)=u^{2}

with

K=2

homogeneous Beta components. Columns vary

a\in\left\{1,2,5\right\}

Figure 36: MSE and runtime versus

S

for

h(u)=u

with

K=2

heterogeneous Beta components with

\pi_{1}=\mathrm{Beta}(2,2)

and

\pi_{2}=\mathrm{Beta}(5,5)

Figure 37: MSE and runtime versus

S

for

h(u)=u^{2}

with

K=2

heterogeneous Beta components with

\pi_{1}=\mathrm{Beta}(2,2)

and

\pi_{2}=\mathrm{Beta}(5,5)

Figure 38: MSE and runtime versus

S

for

h(u)=u

with

K=2

Beta components with

(a_{1},b_{1})=(2,5)

and

(a_{2},b_{2})=(2,5)

Figure 39: MSE and runtime versus

S

for

h(u)=u^{2}

with

K=2

Beta components with

(a_{1},b_{1})=(2,5)

and

(a_{2},b_{2})=(2,5)

Figure 40: MSE and runtime versus

S

for

h(u)=u

with

K=2

Beta components with

(a_{1},b_{1})=(2,5)

and

(a_{2},b_{2})=(5,2)

Figure 41: MSE and runtime versus

S

for

h(u)=u^{2}

with

K=2

Beta components with

(a_{1},b_{1})=(2,5)

and

(a_{2},b_{2})=(5,2)

Figure 42: MSE and runtime versus

S

for

h(u)=u

with

K=2

Beta components with

(a_{1},b_{1})=(0.8,3)

and

(a_{2},b_{2})=(3,0.8)

Figure 43: MSE and runtime versus

S

for

h(u)=u^{2}

with

K=2

Beta components with

(a_{1},b_{1})=(0.8,3)

and

(a_{2},b_{2})=(3,0.8)

Method	Estimation error	Merge (min)	Total (min)
Oracle	0.0000	–	0.9936
Batch	–	–	0.8302
GBF	0.0475	27.2566	28.0868
PMC-OLS	0.0676	0.0050	0.8352
PMC-MLP	0.0101	0.1152	0.9455
PMC-GP	0.0091	0.5056	1.3359

Appendix

Table 4: Results on the MAGIC data.

Figure 44: MSE and runtime versus

S

for

h(u)=u

with

K=2

univariate Gaussian components. Columns and rows vary

\mu_{2},\sigma_{2}\in\left\{1,2,3\right\}

, respectively. The first component is fixed with

\mu_{1}=0

and

\sigma_{1}=1

Figure 45: MSE and runtime versus

S

for

h(u)=u^{2}

with

K=2

univariate Gaussian components. Columns and rows vary

\mu_{2},\sigma_{2}\in\left\{1,2,3\right\}

, respectively. The first component is fixed with

\mu_{1}=0

and

\sigma_{1}=1

Figure 46: MSE and runtime versus

S

for

h(u)=u

with

K=2

bivariate Gaussian components. Columns vary

\mu\in\left\{1,2,3\right\}

. The first component is fixed with

\mu=0

Figure 47: MSE and runtime versus

S

for

h(u)=u\odot u

with

K=2

bivariate Gaussian components. Columns vary

\mu\in\left\{1,2,3\right\}

. The first component is fixed with

\mu=0

Figure 48: MSE and runtime versus

S

for

h(u)=u

with

K=2

multivariate Gaussian components. Columns vary

p\in\left\{2,3,4\right\}

Figure 49: MSE and runtime versus

S

for

h(u)=u\odot u

with

K=2

multivariate Gaussian components. Columns vary

p\in\left\{2,3,4\right\}

Component	Specification
CPU	2 $\times$ AMD EPYC 7763 (128 physical cores, 256 threads, 255 online)
GPU	4 $\times$ NVIDIA RTX A4000 (16 GB memory each)
RAM	$\sim$ 881 GiB usable, shared across all users
OS	Ubuntu 20.04 LTS

Appendix

Table 5: Hardware of the shared machine used for most experiments.

Setting	Value
Methods	PMC-OLS, PMC-MLP, PMC-GP, GBF, and the Oracle
Targets	$h(u)=u$ and $h(u)=u\odot u$
Sample sizes $S$	$100$ , $500$ , $1000$ , $5000$ , $10^{4}$ , $5\times 10^{4}$ , $10^{5}$ , $5\times 10^{5}$ , $10^{6}$ , and, only for the univariate Gaussian location and scale experiment, $2\times 10^{6}$
Realized sample size	$\operatorname{Poisson}(S)$ draws per component
Replications	$20$ independent replications per experimental configuration
Number of components $K$	$2$ , except in the experiments on the number of components

Appendix

Table 6: Settings shared by the baseline and supplementary benchmark suite in Appendix B.1 .

PMC-OLS
Polynomial degree	1
PMC-MLP
Hidden layers	64, 64, 64
Activation	ReLU
Ensemble size	2
Optimizer	Adam

Appendix

Table 7: Common PMC settings used in the baseline and supplementary benchmark suite.

Experiments	Setting	Common value	Value used
Multivariate Gaussian, dimension ( $p=2,3,5,10,20$ )	Kernel	Matérn-3/2	RBF
	Fixed mean	$-0.0014$	$0.0000$
	Fixed lengthscale	6	3

Appendix

Table 8: Baseline and supplementary experiments whose PMC settings differ from the common settings.

Experiment	Method	Cell	Changed settings
Univariate Gaussian, location and scale	PMC-MLP	$S=2{,}000{,}000$ , $\mu_{2}=4$ , $\sigma_{2}=1$ , square	learning rate $6\times 10^{-4}$ , training steps 16,000
	PMC-MLP	$S=2{,}000{,}000$ , $\mu_{2}=4$ , $\sigma_{2}=2$ , square	hidden layers 256, 256, 256, learning rate $2\times 10^{-4}$ , batch size 4096
	PMC-GP	$S=2{,}000{,}000$ , $\mu_{2}=4$ , $\sigma_{2}=1$ , identity	fixed lengthscale 16
	PMC-GP	$S=2{,}000{,}000$ , $\mu_{2}=4$ , $\sigma_{2}=1$ , square	fixed lengthscale 20
Univariate Gaussian, number of components, homogeneous ( $K=20$ )	PMC-GP	$S=10^{6}$ , $K=20$ , identity	kernel RBF, fixed lengthscale 16
	PMC-GP	$S=10^{6}$ , $K=20$ , square	kernel RBF, fixed lengthscale 16

Appendix

Table 9: Per-cell hyperparameter overrides.

PMC-OLS
Polynomial degree	1
PMC-MLP
Hidden layers	64, 64, 64
Activation	ReLU
Discrepancy centering	mean
Optimizer	Adam

Appendix

Table 10: Method settings for the dimension-scaling experiment.

PMC-GP
Kernel	Matérn-3/2
Inducing points	64
Learn inducing locations	No
Learn kernel hyperparameters	Yes
Maximum updates	40
Minimum epochs	20

Appendix

Table 11: Method settings for the calibration experiment.

PMC-GP
Kernel	RBF
Inducing points	64
Learn inducing locations	No
Learn kernel hyperparameters	Yes
Maximum updates	1600
Minimum epochs	10

Appendix

Table 12: PMC-GP settings for the posterior-variance experiment.

Setting	Value
GBF
Output particles	$100{,}000$ (Gaussian and Beta components); $S$ (Gaussian mixtures)
Coalescence time and time mesh	chosen by the SSH guidance of the DCFusion package, with a minimum conditional effective sample size of $0.2$ at initialization
Preconditioning	sample variance (univariate Gaussians) or sample covariance (bivariate and multivariate Gaussians) of each component
Resampling threshold	resample when the effective sample size falls below half of the number of particles
Resampling method	multinomial

Appendix

Table 13: Settings of the fusion methods in the simulation experiments.

Components	Oracle
Gaussian	fit the mean and covariance of each component from its draws and combine them analytically into the Gaussian product
Beta	fit the two shape parameters of each component by maximum likelihood, with location $0$ and scale $1$ fixed, and combine them analytically into the Beta product
Gaussian mixture	fit each component by EM with the true number of modes, a full covariance, $10$ random restarts, at most $1000$ iterations, a tolerance of $10^{-8}$ , and a variance floor of $10^{-8}$ , and compute the product exactly

Appendix

Table 14: Computation of the Oracle in the simulation experiments.

Data
Setting	Respiratory	HIGGS	MAGIC
Observations $n$	$111$	$11{,}000{,}000$	$19{,}020$
Coefficients	$5$	$29$	$11$
Number of batches $K$	$2$	$10$	$10$
Oracle and batch sampler (NUTS)
Prior scale	$2$	$2$	$2$

Appendix

Table 15: Settings of the real-data experiments.

Setting	Respiratory	HIGGS	MAGIC
GBF
Run	Yes	No	Yes
Particles per level	$1000$	–	$1000,1000$
Merge factors	$2$	–	$5,2$
Time mesh	$0.5,1$	–	$0.5,1$
Preconditioning	Yes	–	Yes

Appendix

Table 16: Settings of the real-data experiments (continued).

Setting	Value
All PMC methods
Posterior draws per batch $S$	$40{,}000$
Sampling	independent pairing across batches
Consensus weights	equal ( $1/K$ )
Discrepancy construction	consecutive differences
Targets	$h(u)=u$ and $h(u)=u\odot u$

Appendix

Table 17: PMC settings in the real-data experiments.

Setting	Value
Experiment
Number of images $K$	$10$
Posterior draws per image	$1{,}000$
Parameter	$\theta=19\,\mathrm{sigmoid}(z)$ , with true value $\theta=12$
Image-level prior on $z$	$\mathcal{N}(0,10)$
Prior given all images	$\mathcal{N}(0,1)$

Appendix

Table 18: Settings of the HSP90 experiment.

Setting	Value
PMC-OLS
Polynomial degree	$1$
PMC-MLP
Hidden layers	$32,32$
Activation	ReLU
Optimizer	Adam

Appendix

Table 19: Settings of the HSP90 experiment (continued).

Explore similar work

SCOPE Before You Sample: Conformal Scores Guide Posterior Exploration for Fast and Calibrated Bayesian Predictive Inference
Oct 12, 2026Zahra Moslemi, Ziyi Liang, Babak ShahbabaBayesian Inference Posterior Sampling
University of California, Irvine · University of California, Santa Barbara
On two ways to use determinantal point processes for Monte Carlo integration
Apr 21, 2026Guillaume Gautier, Rémi Bardenet, Michal ValkoDeterminantal Point Process Variance Reduction
†Univ. Lille, CNRS, Centrale Lille, UMR 9189 – CRIStAL, 59651 Villeneuve d’Ascq, France · *Inria Lille-Nord Europe, 40 avenue Halley 59650 Villeneuve d’Ascq, France · ‡DeepMind Paris, 14 Rue de Londres, 75009 Paris, France +1
Picard Proximal Monte Carlo for Parallel Bayesian Imaging with Score-Based Generative Priors
Aug 18, 2026Deliang Wei, Evan Bell, Wenhan Guo +2Posterior Sampling Bayesian Inverse Problems
Johns Hopkins University · University of California, Los Angeles

Oct 12, 2026J/K move · Enter open · S save

SCOPE Before You Sample: Conformal Scores Guide Posterior Exploration for Fast and Calibrated Bayesian Predictive Inference

Zahra Moslemi, Ziyi Liang, Babak Shahbaba

University of California, Irvine · University of California, Santa Barbara

Bayesian Inference Posterior Sampling Conformal Prediction Markov Chain Monte Carlo

Open paper AskPDF arXiv

Bayesian predictive regions account for parameter uncertainty but can be poorly calibrated under model misspecification, even with expensive full MCMC. Cheaper approximations such as variational inference reduce the cost but often yield less accurate predictive distributions. We propose SCOPE (Score-guided Conformal Posterior Exploration), which retains the predictive quality of full MCMC while achieving calibrated coverage at a fraction of its cost. The key observation is that a prediction region can depend on the posterior mainly through a few parameter directions that move the region's boundary, and these are generally not the directions of largest posterior variance. SCOPE gives nonconformity scores a new role in guiding posterior computation: score gradients on a held-out guide set, weighted toward the target boundary, identify the relevant directions. MCMC explores this subspace, while an inexpensive projected variational approximation carries the remaining uncertainty. The same guide set calibrates the region's final density level. We establish asymptotic coverage for SCOPE and provide a split variant with finite-sample validity. Controlled simulations show that score-guided exploration yields more accurate calibrated regions than non-targeted direction choices at low active rank. Across simulated and real benchmarks, SCOPE has a median test-NLL gap of only 0.002 nats relative to full Hamiltonian Monte Carlo (HMC), and on the real benchmarks it runs 12-15 times faster.

Abstract

Explore similar work

SCOPE Before You Sample: Conformal Scores Guide Posterior Exploration for Fast and Calibrated Bayesian Predictive Inference

On two ways to use determinantal point processes for Monte Carlo integration

Picard Proximal Monte Carlo for Parallel Bayesian Imaging with Score-Based Generative Priors