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

Oct 12, 2026

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

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.

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

Ask PDF

Apr 21, 2026

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

The standard Monte Carlo estimator

\widehat{I}_N^{\mathrm{MC}}

\int fdω

relies on independent samples from

ω

and has variance of order

1/N

. Replacing the samples with a determinantal point process (DPP), a repulsive distribution, makes the estimator consistent, with variance rates that depend on how the DPP is adapted to

f

and

ω

. We examine two existing DPP-based estimators: one by Bardenet & Hardy (2020) with a rate of

\mathcal{O}(N^{-(1+1/d)})

for smooth

f

, but relying on a fixed DPP. The other, by Ermakov & Zolotukhin (1960), is unbiased with rate of order

1/N

, like Monte Carlo, but its DPP is tailored to

f

. We revisit these estimators, generalize them to continuous settings, and provide sampling algorithms.

Guillaume Gautier, Rémi Bardenet, Michal Valko

†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

Determinantal Point Process Variance Reduction

Ask PDF

Aug 18, 2026

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

Bayesian imaging inverse problems often require sampling from high-dimensional posterior distributions. While recent score-based and diffusion models provide expressive Bayesian priors, their sampling procedures remain inherently sequential and computationally expensive for large-scale imaging applications. We propose PiX-MC, a time-parallel posterior sampling framework based on proximal Langevin dynamics and Picard iteration. The proximal-likelihood formulation exploits the fact that many imaging likelihoods admit efficient, problem-specific proximal operators, while Picard refinement exposes parallelism across discretization nodes and naturally supports multi-GPU implementation. To further improve practical scalability and sampling performance, we develop multi-block and annealed variants of the proposed framework. We establish convergence guarantees under transparent assumptions, accommodating non-log-concave posteriors, imperfect learned score models, multi-block implementations, and annealing schedules. Experiments on a diverse collection of imaging inverse problems demonstrate that PiX-MC substantially reduces wall-clock time while preserving reconstruction quality. On a

512\times512\times80

sparse-view computed tomography (CT) problem, annealed multi-block PiX-MC achieves up to a

50\times

runtime speedup over the standard Langevin sampler using eight GPUs.

Deliang Wei, Evan Bell, Wenhan Guo +2

Johns Hopkins University · University of California, Los Angeles

Posterior Sampling Bayesian Inverse Problems Diffusion-Based Inverse Problems Langevin Dynamics

Ask PDF