Joint Distribution Estimator

• Joint distribution estimator is a statistical construct that estimates the joint probability law of multiple random variables using a range of parametric, nonparametric, and algorithmic methods.
• It leverages techniques such as tensor-product spline MLE, semi-supervised plug-in methods, and scalable kernel strategies to achieve consistent and optimal convergence rates.
• Its applications span biostatistics, causal inference, and privacy-preserving data analysis, though challenges remain for managing high dimensionality and strict regularity assumptions.

A joint distribution estimator is a statistical or algorithmic construct designed to estimate the joint probability law of two or more random variables based on observations, either fully or partially observed. Estimation methods for joint distributions span parametric, nonparametric, semiparametric, and algorithmic approaches, with particular attention to high-dimensional scaling, structural constraints (such as monotonicity or sparsity), privacy preservation, and the presence of missing, censored, or noisy data. Modern developments emphasize tractability, convergence properties, and real-world applicability across statistics, machine learning, privacy, and causal inference.

1. Spline-Based Sieve Maximum Likelihood for Joint CDFs with Censored Data

A theoretically and numerically robust approach for joint cumulative distribution function (CDF) estimation under bivariate current status data employs tensor spline-based sieve maximum likelihood estimation (MLE), as developed in (Wu et al., 2012). The method constructs an approximant for the joint CDF $F(s,t)$ and marginals $F_1(s)$, $F_2(t)$ via tensor–product B-spline bases:

$$F_n(s,t) = \sum_{i=1}^{p_n} \sum_{j=1}^{q_n} \alpha_{ij}\, N_i^{(1,l)}(s) N_j^{(2,l)}(t)$$

Coefficients are determined by maximizing the likelihood subject to monotonicity and agreement constraints. Crucially, the use of I-splines, which are cumulative sums of B-splines, embeds the monotonicity constraint at the basis level, converting the estimation problem into a well-conditioned numerical optimization under linear inequality constraints.

Estimation proceeds with a generalized gradient projection algorithm:

• Compute Newton-type step and project onto the feasible region defined by active constraints.
• Find a maximal feasible step size constrained by nonnegativity and linear sum.
• Conduct a line search to guarantee increase in log-likelihood.
• Update constraints and check for convergence.

The resulting estimator is consistent and achieves an $L_2$ convergence rate $O_p(n^{-\min\{pv,(1-2v)/2\}})$, potentially exceeding $n^{-1/3}$ for sufficiently smooth functions. Monte Carlo simulations using data generated from a Clayton copula with exponential marginals show uniformly reduced bias and mean squared error relative to traditional NPMLEs. This estimator is uniquely defined, computationally efficient, and naturally extends to interval-censored data, supporting practical applications such as bivariate event time modelling in biostatistics.

2. Joint Distribution Estimation under Semi-Supervised and f-Divergence Loss Frameworks

Recent work has formalized the semi-supervised joint distribution estimation problem, where $m$ samples of $(X,Y)$ pairs and $n$ samples with $X$ alone are observed. The minimax framework with $l_p^p$ and $f$-divergence losses has received particular attention (Erol et al., 15 May 2024, Erol et al., 2023).

The composite estimator is constructed as follows:

• Estimate the marginal $p_X$ using all available $X$ samples (both labeled and unlabeled), leveraging the larger (typically unlabeled) sample size for minimax-optimal marginal estimation.
• Estimate the conditional $p_{Y|X}$ by subdividing labeled pairs by value of $X$ and applying minimax univariate estimators (with optimal constants) on the $Y$-conditional subsamples.
• Form the estimated joint pmf as $\hat{p}_{X,Y} = \hat{p}_X \cdot \hat{p}_{Y|X}$.

It is rigorously established that this plug-in composition estimator achieves first-order minimax optimality for $1 \leq p \leq 2$, including for $l_1$ (total variation) loss, and for common $f$-divergences: Kullback-Leibler, $\chi^2$, Squared Hellinger, Le Cam, under mild regularity. The approach is robust for $m = o(n)$, and risk bounds explicitly quantify error contributions from both labeled and unlabeled samples. This result has broad implications for practical data analysis settings with partial labeling, showing that optimally leveraging both data sources yields joint estimators that attain the minimax rates with explicit first-order constants.

3. Algorithmic and Computational Considerations

Modern joint distribution estimation must contend with high dimensionality, privacy, and computational feasibility. Several methodologies address these challenges:

• Spline sieve estimators (Wu et al., 2012): The choice of tensor-product spline degrees and knot placement determines the approximation's bias-variance trade-off. Sieve dimension should scale as $n^v$ ($v<1$) for optimal convergence, and the I-spline reparameterization enables large-scale numerical optimization under shape constraints.
• Tensor decomposition with dictionaries (Haque et al., 2022): For multivariate joint densities, especially in hybrid discrete-continuous settings, Canonical Polyadic Decomposition (CPD) with per-axis dictionary mixtures allows flexible, sparse representations. Factor densities for each variable are expressed as mixtures from user-specified families (e.g., Gaussian, Laplacian, Uniform). Fitting is performed by minimizing Frobenius-norm errors to observed 2D marginals, subject to simplex constraints via projected mirror descent and SPA-based identifiability refinement.
• Privacy-preserving estimation: For differential privacy, joint distributions are reconstructed by imposing consistency with differentially private query answers while minimizing KL divergence to a public empirical prior (Tao et al., 2021). Optimization is iterative and alternates between enforcing workload consistency and minimal divergence updates. In local differential privacy, noisy multi-dimensional histograms are inverted with analytically derived means/covariances and matrix inversion, exploiting normal approximations to yield unbiased frequency (co-occurrence) estimates (Carette, 2022).
• Scalable kernel methods: Adaptive joint distribution learners express the Radon–Nikodym derivative of the joint with respect to the product of marginals in a tensor product RKHS, with efficient low-rank Cholesky decompositions for scalability (Filipovic et al., 2021).

4. Constrained, Robust, and Semi/Nonparametric Estimators

• Constrained joint estimators in regression: Joint asymptotic distributions for unrestricted and linearly restricted estimators in multivariate regression are characterized via exact central limit theorems. Asymptotic Distributional Risk (ADR) theory demonstrates that imposing a valid restriction reduces risk near the restriction, but if misspecified, restricted estimators can incur greater risk than unrestricted ones (Nkurunziza et al., 2017).
• Semiparametric and robust estimators: In complex elliptical distributions, joint estimation of location and scatter/shape matrices is accomplished by combining Tyler's robust $M$-estimator of location with the $R$-estimator of shape (which achieves semiparametric efficiency). The block-diagonal nature of the semiparametric Fisher information allows the two estimation problems to decouple asymptotically, yielding near-optimal performance for the location and strictly improved mean squared error for the scatter, especially in heavy-tailed or non-Gaussian settings (Fortunati et al., 2021).
• Infinite-dimensional joint estimation: For stochastic parameters in PDE-constrained mixed effects models, nonparametric MLEs over the space of probability measures (with existence and consistency established in the Prohorov metric) are computed via discrete measure approximations coupled to numerical solutions (e.g., Galerkin methods). Cross-validation and uncertainty quantification are performed via the propagation of parameter distributions through the model (Asserian et al., 2023).

5. Monte Carlo and Bootstrap Approaches for Joint Law Sampling

• Augmented estimator distributions: In high-dimensional regularized regression (e.g., Lasso), the paper (Zhou, 2014) demonstrates that the joint law of the estimator and its subgradient can be expressed in closed form, enabling efficient Markov chain Monte Carlo (MH Lasso sampler) and importance sampling procedures for sampling from the exact or estimated (plug-in) joint distribution. This is essential for accurate inference, uncertainty quantification, and hypothesis testing in sparse models.
• INAR (integer-valued autoregressive) time series: A joint semi-parametric estimator for both the model's coefficients and the discrete innovation distribution is constructed via Z-estimation. The corresponding bootstrap procedure mimics this joint estimation, yielding asymptotically valid inference for model parameters and innovation law even without any parametric assumptions on the innovations (Faymonville et al., 15 Jul 2025). Applications include joint model diagnostics and predictive interval construction for count time series.

6. Practical Applications

Joint distribution estimators are used in diverse settings:

• Causal inference: Identification and estimation of joint laws of potential outcomes or principal stratification subgroups, where joint laws are generally unidentified but can be recovered by combining data from multiple experiments under transportability constraints and solvability conditions (Wu et al., 29 Apr 2025).
• Multiple test calibration: Recursive computation of joint order statistics enables exact error rate calculation in procedures such as false discovery rate (FDR) control (Schroeder et al., 2018).
• Survival analysis and biostatistics: Joint spline-based CDF estimation for bivariate time-to-event data supports coherent modeling of complex dependencies under censoring.
• Machine learning and privacy: Kernel-based, tensor-based, and privacy-preserving estimators are at the core of robust generative, classification, synthetic data generation, and query answering pipelines, especially under privacy constraints or in high-dimensional, partially labeled scenarios.

7. Impact and Limitations

Joint distribution estimation is foundational for inference, prediction, and causal evaluation in modern statistics and machine learning. Spline-based and tensor decomposition methods address the curse of dimensionality via flexible, low-dimensional representations or structured regularization. Plug-in composition estimators are theoretically optimal for semi-supervised discrete estimation in a variety of distance/divergence losses.

Limitations arise from the modeling choices and regularity assumptions: smoothness constraints can hamper adaptability to distributions with sharp features; the curse of dimensionality persists in fully nonparametric settings unless strong structure is leveraged; certain consistency results require large sample sizes or regular design. In privacy-preserving contexts, approximation error is dictated by the noise mechanism and the invertibility of the underlying randomization matrix.

Emerging research continues to focus on the interplay between statistical efficiency, computational scalability, structural expressivity, and privacy or causal identification constraints—often necessitating intricate combinations of numerical optimization, modern empirical process theory, and interdisciplinary algorithmic design.