Papers

Modern learning problems increasingly leverage data from heterogeneous sources to improve robustness and generalizability under distribution shift, often with high-dimensional parameters admitting low-rank structure. Existing group distributionally robust optimization (group-DRO) approaches optimize against adverse source mixtures, but extending this principle to low-rank matrix estimation with a broad class of losses leads to a challenging nonconvex-concave problem for which methods and theory remain limited. We propose Prior-Informed generalizable Low-rank representation Learning via Adversarial Reweighting (PILLAR), which smooths the group-DRO objective by regularizing adversarial source weights toward a user-specified reference distribution while learning a shared low-rank representation. By tempering extreme adversarial reweighting, PILLAR improves estimation stability and reduces the potential conservativeness of group-DRO. Using Burer-Monteiro factorization, we optimize the nonconvex-concave objective by bifactored gradient descent with mirror ascent. We establish linear convergence to a statistical neighborhood of the population saddle point and derive an error bound that separates low-rank estimation from adversarial reweighting. We further show that PILLAR improves on vanilla group DRO when the latter's excess conservativeness dominates source heterogeneity and prior misspecification. Simulations and an electronic health record study of suicide-attempt risk demonstrate improved worst-source and transported predictive accuracy.

What does it mean to measure the replicability of scientific results? Is it merely determining whether results are consistent across studies? Much of the existing literature on replicability draws conclusions exclusively from the comparison of study results, in terms of observed estimates, uses a variety of different metrics, and ultimately collapses the comparison into a binary verdict of success or failure. This fails to differentiate between repeated errors and corroboration of the shared generative effect. We develop a unified Bayesian framework in which distinct but related notions of replicability are formally defined and jointly embedded within a single inferential structure. A single model answers four complementary questions: agreement among independent studies; consistency between the individual studies and the shared generative effect; agreement between an existing body of evidence and a new study from the same super-population; inference on the shared generative effect. In practice, we suggest quantifying replicability through posterior probabilities relative to regions of practical relevance. We illustrate these ideas in an eQTL (expression quantitative trait locus) analysis assessing the replicability of a regulatory effect across human ancestral groups.

Quasi-periodic oscillations are characterized by fitting Lorentzian functions to a nonparametric estimate of the power spectral density, a two-step procedure that separates estimation from modeling and does not accommodate measurement error. We develop a parametric alternative based on autoregressive moving-average (ARMA) processes observed with measurement error. We derive an exact decomposition of the ARMA spectral density into real-valued components, one for each class of autoregressive roots under complex conjugation, plus a polynomial term when the moving-average order is at least the autoregressive order. Each component is a Wrapped Cauchy kernel, and for complex roots it is accompanied by a conjugate-Poisson term that makes the peak asymmetric and displaces its maximum from the frequency given by the argument of the root. We obtain closed-form expressions for the location, height, width, and asymmetry of each peak, and a matching variance decomposition. Parameters are estimated by maximum likelihood through the Kalman filter, and all descriptors follow by plug-in. Simulations show that the closed-form peak location is recovered without appreciable bias. An X-ray binary light curve illustrates the methodology.

Quantile regression has received considerable attention in high dimensional methodology due to its robustness and ability to accommodate heteroskedasticity. Despite this attention, current inferential methods for high dimensional quantile regression have significant limitations—only testing significance, testing at a single quantile at a time, or assuming homoskedasticity. To generalize the inferences that can be drawn in high dimensional quantile regression, we introduce partial penalized Wald and rank score tests for testing linear hypotheses uniformly over a set of quantiles $\mathcal T$ under a heteroskedastic model. We propose globally adaptive local linear approximation (LLA) algorithms for folded concave penalized estimation and show that our global LLA algorithms recover “testing oracles” for the full and reduced models with high probability. We then show that our global LLA estimators achieve the $\sqrt{s/n}$ oracle estimation rate uniformly over $\mathcal T$, establish Bahadur representations, and derive asymptotic distributions for the test statistics under the null and local alternatives. We develop a dual formulation incorporating adaptive penalty weights for the rank scores and use the strong oracle property to establish their oracle equivalence. Using empirical process techniques, we derive reference distributions based on suprema of squared Bessel processes, extending partial penalized inference beyond pointwise chi-square calibration.

Treatment decisions are often sequential and depend on information absent from the recorded history. Latent sequential ignorability models such assignment by assuming that baseline covariates and principal strata, defined by joint potential intermediate responses, account for confounding of later treatment. For two treatment periods with a binary intermediate response, we derive sharp population and principal effect bounds under a specified association between the potential responses and show that the two targets can require different optimizations to obtain sharp bounds. For population effects, optimization over latent treatment probabilities reduces to their feasible endpoints. Stronger selection between response types can be offset by making their outcome distributions more similar, preserving both the observed distribution and population marginal outcome distributions. This replacement can change effects within principal strata, whose sharp bounds may instead require interior treatment probabilities. Restrictions on differences in conditional effects can tighten the bounds without equating outcome levels. For binary outcomes, we give an exact finite algorithm that includes these interior assignments. We also develop inference under smoothness and estimation conditions. We illustrate the bounds in a study of gym attendance after an incentive program.

Sharing time-to-event data can leak sensitive information about individual participants in clinical studies. To protect the privacy of data, differential privacy (DP) framework offers rigorous protection for data sharing. The local variant of differential privacy is particularly appealing because each record is perturbed before it leaves the data owner or contributor. However, the injected perturbation distorts the true distribution of the source data, yielding low utility and invalid statistical inference with the released data. In this work, we propose a nonparametric estimator for the true survivor function via the deconvolution kernel function with locally differential private (LDP) right-censored survival data in which both the event times and the event indicators are privatized. Our estimator simultaneously de-biases the perturbation applied to the times and the errors induced by flipping the event indicators due to LDP. Theoretically, we establish the uniform consistency and point-wise asymptotic normality of the de-biased distribution estimator, and derive its uniform convergence rate together with the optimal bandwidth attaining that rate. Simulation studies and the analysis of a semi-synthetic dataset generated from the METABRIC breast cancer cohort demonstrate the performance of the proposed estimator across privacy budgets, censoring levels, and sample sizes. By the post-processing property of DP, the proposed de-biased estimator inherits the privacy guarantee while restoring valid inference on the survivor function.

The cvxoptdes R-package provides a unified framework to perform convex optimal experiment design and online response optimization for linear, generalized linear and additive, nonlinear, and random intercept mixed models, using exclusively open-source solvers. It supports optimization and augmentation of weighted, penalized and convex constrained continuous designs across a discrete grid of candidate points and rounding to integer-valued designs using point-exchange algorithms or mixed-integer programming for orthogonally constrained designs, including the class of orthogonal minimally aliased response surface (OMARS) designs. Models and design constraints are defined and updated in an object-oriented context with the aim to facilitate the process of sequential (penalized) design augmentation analogous to online design procedures known from active learning and Bayesian optimization.

Rating procedures are crucial in psychological and medical diagnosis, educational assessment, peer review, and grant evaluation. A major methodological challenge is disentangling multiple sources of variability in ratings, such as subject ability, rater behavior, and item characteristics. Item Response Theory-based models, most notably many-facet Rasch models and hierarchical rater models, have been widely adopted for jointly estimating subject, rater, and item effects, while retaining high interpretability. We introduce the Dynamic Extended Generalized Many-Facet Rasch model, a general class of Bayesian semiparametric models that integrate and extend existing approaches. Its measurement model allows location and discrimination parameters to depend jointly on item and rater characteristics and incorporates subject-item and rater-subject variation. The dynamic extension introduces a rater-specific cross-lagged effect capturing dependence on the previously rated subject, and a Bayesian P-spline capturing systematic changes in discrimination over the rating sequence, accounting for learning, fatigue, or adaptation effects. Subject heterogeneity is addressed through Bayesian nonparametric mixtures, flexibly learning the latent ability distribution without strict parametric assumptions and enabling model-based subject clustering. We also derive a closed-form, model-based inter-rater consensus coefficient and study its characteristics. We discuss theoretical properties, develop posterior computation strategies, and evaluate the model in simulation studies and two real-data applications.

Many scientific models can be simulated but have intractable likelihoods, making frequentist uncertainty quantification difficult. We construct exact anytime-valid confidence sequences for the parameters of such models. At each observation, fresh simulations from each candidate null conformalize a positive predictable score into a conditional e-value; their product is a test martingale, which we invert over the parameter. Validity holds for any such score, whatever its accuracy. We show that the oracle likelihood-ratio score maximizes expected log growth, which motivates approximating it with an amortized neural ratio estimator. On synthetic benchmarks, an audit of model routing in an LLM API, and microlensing data, the resulting confidence sequences stay valid under continuous monitoring, where fixed-sample methods do not, and neural scores remain competitive with exact likelihood ratios and application-specific scores.

Latent field models are essential for spatial statistics, but the posterior distribution of the latent field, even with a Gaussian process (GP) prior, is analytically intractable whenever the likelihood is non-Gaussian. Existing approximations are either too slow or too restrictive for large non-Gaussian datasets. We present a scalable variational inference framework that mirrors the sparse nearest-neighbor autoregressive structure of the prior in the variational family, uses normalizing-flow conditionals for expressiveness, and amortizes all variational parameters via a single shared neural conditioner. Doubly stochastic optimization with reduced ancestor sets achieves polylogarithmic complexity per iteration.

Circular data recorded over time, such as wind direction, are increasingly observed in intervals rather than as single point-valued observations, yet no existing framework accommodates data that are simultaneously circular, interval-valued, and temporally dependent. Time series models for point-valued circular data, including recently proposed score-driven approaches, assume that each observation is a single angle. This paper addresses this gap by proposing a dynamic conditional score framework for interval-valued circular time series. The lower and upper boundary sequences of the observed intervals are treated as circular time series governed by score-driven dynamics, which preserves the periodicity of the circular domain and the dependency of boundaries, in which the boundaries are coupled through a toroidal distribution and a jointly derived score. Maximum likelihood estimation is discussed, and their finite-sample performance is evaluated through simulation. An application to interval-valued wind direction data illustrates the practical value of retaining the interval structure within a fully rotationally invariant circular framework.

The external validity of a study is the extent to which its findings generalize to unseen experimental conditions, such as datasets. We formalize the results of a study together with its research question, and define external validity as the probability that two independent and same-sized replications of the study produce results whose maximum mean discrepancy (MMD) is below a given tolerance. This approach, centered on the choice of a space of results and of an appropriate reproducing kernel for the MMD, is flexible enough to be deployed on multiple different types of results and research questions. Additionally, we establish an asymptotic power-law relation between the size of a study and the quantiles of the MMD, which we use to estimate the desired size of an externally valid study from a small-scale preliminary study. We deploy our framework to rankings with ties and discuss kernels suited to different research questions. Experiments on existing case studies and on a simulated newly designed study illustrate the applicability of the framework. We release our code as part of the genexpy Python module.

In spatial data analysis, global measures of spatial correlation provide useful insights into the general spatial patterns, but fail to capture local variations. Local Indicators of Spatial Association (LISA) address this limitation by identifying local clusters and spatial outliers, offering valuable applications in many fields of research such as public health, urban planning and environmental monitoring, to name only a few. However, the application of LISA tools in hypothesis testing is undermined by the presence of the multiple testing problem which dramatically increases the risk of committing Type I errors especially when dealing with very large datasets. Standard corrections procedure (like, e. g., Bonferroni and Benjamini-Hochberg procedures) fail to consider the spatial dependence among the various LISA tests and, furthermore, they are overly conservative. This paper proposes an adjustment that, by incorporating the notion of the "effective sample size", improves statistical inference balancing false discovery rate control and statistical power.

In this paper we propose and study arctan-skewing, a general approach to skew any symmetric distribution on the unit hypertorus. Arctan-skewing retains the main advantages of the popular sine-skewing mechanism, but improves on it by overcoming its two major limitations. These are a condition on the skewness parameter and a singular Fisher Information matrix in the vicinity of symmetry when combined with well-known symmetric models from the literature such as the circular von Mises or bivariate Cosine models. Arctan-skewing thus represents a compelling alternative to model skewed data on the circle and torus, notably from the field of protein bioinformatics. After a thorough investigation of its stochastic properties, we present a maximum likelihood estimation procedure and devise optimal tests for symmetry against our new class of skewed alternatives. Finally we show that arctan-skewed models better fit protein data compared to their sine-skewed counterparts.

Threshold and change-plane models provide a useful framework for capturing regime-dependent dynamics, but most existing approaches typically rely on scalar thresholds or linear splitting boundaries. We propose a scalable semiparametric change-plane regression framework that allows the regime boundary to depend nonlinearly on multivariate state variables. Random Fourier features provide a flexible and computationally efficient approximation to the nonlinear boundary function. A two-stage procedure first estimates the induced regime classification using a smooth surrogate and then refits the structural parameters by hard-threshold OLS. We establish convergence of the first-stage regime classification and show that its estimation error is asymptotically negligible for the second-stage refit, yielding oracle equivalence. We further develop a dependent-multiplier bootstrap Sup-Wald test for detecting the existence of regime effects when the boundary is unidentified under the null, and establish its asymptotic validity. Simulations demonstrate accurate regime recovery, near-oracle performance of the proposed estimator, and satisfactory finite-sample behavior of the Sup-Wald test. An application to postwar U.S. unemployment data provides strong evidence of regime dependence, while the estimated nonlinear change plane identifies a distinct stress regime characterized by significant shifts in short-run unemployment dynamics.

Simulation-based inference (SBI) enables Bayesian inference when a simulator is available but the likelihood is intractable. Modern SBI methods train neural networks on simulated data to obtain an amortized approximation of the posterior, making inference for new observations fast. A central challenge in SBI is to extract informative, low-dimensional summary statistics; without them, performance can degrade rapidly as the observation dimension grows. We therefore highlight summary learning as a critical component of SBI and propose a general framework that jointly learns a neural-network-parameterized summary and an amortized posterior estimator. We further introduce a local refinement strategy that adapts the amortized estimator to a specific observation, thereby reducing the amortization gap. Both amortized training and local refinement accommodate implicit posterior estimators without tractable densities. On the theoretical side, we analyze the finite-sample behavior of the joint learning procedure through the lens of nonparametric conditional distribution estimation, and show on representative examples how low-dimensional summaries can mitigate the classical curse of dimensionality. Empirically, our results support this motivation and demonstrate favorable performance relative to existing methods.

Effect sizes for single-case experimental designs usually rely on models for the baseline trend and the serial dependence, which are challenging to check with few sessions. We develop design-based effect sizes for multiple-baseline experiments with a randomized intervention onset. We define case-level and pooled standardized effects with potential outcomes and propose weighted estimators that are unbiased or asymptotically unbiased without an outcome model. The usual difference between phase means is generally biased when the untreated outcome changes over sessions. For research synthesis, we provide variance estimators and asymptotic confidence intervals that require no model for serial dependence, and variance formulas for planning. Simulations and three randomized studies of reading and writing interventions illustrate the methods.

Spatial prediction is a critically important task in data analyses that involve spatially or spatio-temporally-referenced data. It is commonly carried out using hierarchical models, which comprise a process model describing the latent process of interest, and a data model describing how the process is observed. The target quantity is the predictive distribution, which is the conditional distribution of the latent spatial process given the observed data. Unfortunately, the predictive distribution is rarely available in closed form and can be difficult to sample from. Here, we show how recent developments in guided diffusion modelling from machine learning can address a wide class of problems in spatial prediction when the likelihood function derived from the data model is known and differentiable, possibly after relaxation. To show the versatility of the approach, we apply it to both Gaussian and non-Gaussian process models and a variety of data models that describe direct partial observations, linear and areal aggregates, counts, right censoring, and area-exceedance functionals. Where possible, we validate the framework against exact and Markov-chain-Monte-Carlo (MCMC)-derived predictive distributions. We demonstrate our approach on several simulation studies and in a controlled experiment of methane-flux inversion over the UK and Ireland. We show that guided diffusion can sample from the correct predictive distribution in a variety of settings, that it can be applied to models where MCMC is intractable, and that it scales well with the dimensionality of the latent process.

Comparing distributions under severe sample imbalance requires effectively using abundant majority observations while retaining sensitivity to differences involving the smaller sample. This challenge is particularly relevant for complex, non-Euclidean data, such as graphs, whose structural characteristics call for flexible testing procedures. Classifier two-sample tests (C2STs) provide such flexibility, but existing procedures are primarily designed for balanced or comparably sized samples. In this paper, we study C2STs in a highly imbalanced regime, and show that severe sample imbalance may substantially degrade testing performance. Different from conventional imbalanced learning approaches that treat sample imbalance merely as a nuisance for classifier training, we turn it into a source of repeated testing evidence. Building on this insight, we propose an Enhanced Imbalance-Aware Deep Classifier Two-Sample Test (EIAD-C2ST), a classifier-agnostic framework designed for sample imbalance that accommodates diverse data modalities by tailoring classifiers to data structure. For computationally expensive deep classifiers, we further develop a scalable representation-amortized implementation strategy. Theoretically, we establish rigorous guarantees on asymptotic size and characterize conditions for power consistency, allowing the imbalance ratio to diverge. We demonstrate the empirical gains using simulations and an application to blockchain forensics. The results reveal that the majority-sample surplus, when properly leveraged, can be organized into useful statistical evidence, rather than an obstacle, for distributional comparisons with complex data.

Conformal prediction provides finite-sample, distribution-free coverage guarantees that can support decision-making in specific settings. In multistage problems, earlier decisions determine the populations relevant to later calibration. When later scores are refit using earlier calibration results, reusing those cases for ranking need not preserve exchangeability with a surviving test case. We study sequential conformal calibration with a fresh fold at each stage: conditional on earlier folds, the selected fold scores and the surviving test score are exchangeable. This yields finite-sample coverage guarantees for every prefix of the decision sequence. For atomless survivor score distributions, we derive a universal Binomial-thinned Beta chain for the realized coverage process. Its law depends only on fold sizes and target levels, even under predictable refitting, and yields exact coverage slack, abstention probabilities, and calibration-conditional guarantees. We further characterize the tradeoff between freshness and data reuse. For fixed scores, strict-survivor reuse is exactly valid. For adaptive scores refit using past calibration data, reusing survivors can incur worst-case rank loss, while blocked designs provide an intermediate option between full reuse and stagewise fresh folds. As a case study, we apply the framework to serial multi-echelon inventory problem and propose a sequential conditional-quantile formulation. Each stock level is learned as a context-dependent quantile of cumulative demand among scenarios satisfying earlier no-backlog constraints. We connect these targets to common-shift stationarity and derive decision-level implications for the resulting calibrated policies. In a three-echelon simulation, fresh-fold refitting achieves 81.98% final-prefix coverage against an 81.82% target while reducing mean cost by 5.5% relative to the Bonferroni baseline.

We study the asymptotic properties of the nonlinear stochastic GMM algorithm proposed by Chen, Lee, Liao, Seo, Shin and Song for recursively estimating parameters defined by possibly nonlinear and overidentified moment restrictions $\E[g(X,\beta^*)]=0$. Each observation updates estimates of the moment Jacobian, inverse covariance matrix and Gauss–Newton preconditioner at a cost independent of sample size. For linear and nonlinear instrumental-variable moments, we derive an exact rank-two Sherman–Morrison–Woodbury update. The algorithm combines preconditioned stochastic approximation with iterate averaging. Under a preconditioned monotonicity condition, satisfied globally in linear models and locally around a consistent pilot estimate in general, we establish almost sure convergence and rates, as well as a functional central limit theorem with the efficient GMM variance. These results require four moments of the moment function and its Jacobian when the step exponent exceeds $2/3$. The inverse Gauss–Newton matrix consistently estimates this variance, and a projected linearised statistic enables online testing of overidentifying restrictions. We extend the results to $L$-dependent data with long-run variance weights under stronger moment assumptions. Simulations show that the averaged estimator approaches offline two-step GMM in accuracy and coverage at $10^6$ observations, while requiring a fraction of its computing time and fixed memory at up to $10^7$ observations. The online test achieves correct size and comparable power. In an application estimating the effect of a third child on weeks worked using $254{,}654$ women from the 1980 US Census, single-pass estimates agree with offline GMM within a fraction of a standard error on average across data orderings, and both methods yield the same overidentification-test conclusion in most orderings.

Multivariate longitudinal count data arise in biological studies, where responses may be over-dispersed, correlated within subjects and over time, measured at irregular times, and exhibit unobserved subgroup structure. Our goal was to develop a model that can accommodate over-dispersion, longitudinal dependence, cross-outcome dependence, covariate effects, irregular observation times, and covariate-dependent latent-class membership. We propose a latent-class Poisson-lognormal mixed model for clustering heterogeneous multivariate count trajectories. Estimation is performed using a variational EM algorithm designed for scalable inference. Simulation studies show accurate parameter estimation, strong latent-class recovery, and reliable selection of the number of classes, as compared to state-of-the-art approaches. We also applied the proposed approach to longitudinal neutrophil RNA-seq data from hospitalized COVID-19 patients as an illustration. The proposed framework provides a flexible approach for clustering complex longitudinal biological count data. Implementation code is available at https://github.com/PengqiLiu-Mc/latent-class-mpln-mixed-models.

Bayesian simulation-based inference (SBI) can struggle to recover increasingly concentrated posteriors as observation size grows. We consider learning sufficient summary statistics jointly with the posterior as a general principle for scalable inference in various SBI methods. A common measure taken in SBI is dimension reduction through predetermined summaries, which may partially discard posterior information. Our approach instead imposes structure on how the posterior learner depends on a learned summary (e.g., affinity), and learns the summary simultaneously with the posterior. We identify an important mechanism named two-way guidance: under a structural restriction on the summary, the posterior sufficiency of the learned summary simplifies posterior learning, whereas a suitable posterior objective guides the learned summary towards posterior sufficiency. This mechanism explains a counterintuitive finding that the structural restriction on the summary statistics improves posterior learning, rather than limit its performance. We establish theoretical guarantees on both posterior estimation error and the learned summary's departure from sufficiency, thus identifying conditions on the simulation budget for recovering increasingly concentrated posteriors as observation size grows.

Estimating (heterogeneous) treatment effects is still challenging when the outcome is semicontinuous (i.e., exhibiting both a discrete point mass at zero and a skewed continuous distribution for positive values), despite the strong presence of this data type in real world settings. Existing methods either ignore the zero mass, impose fragile parametric assumptions, or suffer from extreme computational costs and data augmentation noise. This paper proposes Zero-Inflated Continuous Bayesian Causal Forests with Duan's Smearing (ZIC-BCF-Smear), a decoupled two-part nonparametric estimator that fits a Probit BCF on the full sample to model participation and a Gaussian BCF strictly on the active subpopulation to model positive intensity. Subpopulation Propensity Adjustment (SPA) mitigates regularization-induced confounding within the active subpopulation, and Duan's non-parametric smearing estimator re-transforms predictions to the response scale, bypassing the fragile scale multipliers of parametric models. A comprehensive simulation benchmark shows that ZIC-BCF-Smear matches well-specified parametric models without their assumptions and remains stable and calibrated where competitors either diverge or pathologically over-cover for nonlinear data. The method is then applied to individual-participant data from the SPECTRUM randomized phase III head-and-neck cancer trial, estimating the causal effect of assignment to panitumumab plus chemotherapy on the number of days patients spend with severe toxicity, an outcome for which a quarter of patients contribute an exact zero. Panitumumab assignment raises that burden by about ten days per patient under the trial's own differential treatment duration, and the decoupled structure localizes the effect mainly in the positive-intensity part rather than in whether the reconstructed severe-AE burden is positive.

We develop a class of semiparametric zero-altered models (ZAMs) for count data, estimated through the generalized method of moments (GMM). The proposed framework extends the utility of ZAMs by allowing for both under- and over-dispersion while preserving support on the entire set of natural numbers. This offers an alternative to other models that can experience computational challenges (e.g. Conway-Maxwell-Poisson) or whose supports can be truncated as functions of the parameter space (e.g. generalized Poisson). By formulating estimation through moment restrictions rather than full likelihood specification, the GMM approach provides substantial robustness against model misspecification. In particular, inference based on likelihood methods for parametric generalized linear models (GLMs), such as those built on Poisson, negative binomial, or related count distributions, can become unreliable when the assumed distributional form is violated, often leading to biased estimates and invalid inference. In contrast, our methodology requires only the first two moments of the distribution, ensuring consistent estimation even when higher-order distributional features are misspecified. We establish the asymptotic properties of the estimators, including consistency and asymptotic normality, under general regularity conditions. The flexibility of the approach is illustrated through both simulated experiments and real-data case studies, demonstrating superior performance relative to popular parametric GLMs, particularly in settings characterized by distributional misspecification. To facilitate practical adoption, we provide accompanying R functions for easy implementation, enabling practitioners to apply the proposed methods in applied settings without requiring bespoke programming.

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.

Leave-one-out (LOO) cross-validation, a standard tool for Bayesian model assessment, evaluates each observation using a LOO predictive distribution derived from the model itself. In this work, we introduce an alternative construction in which posterior predictive simulation connects the Bayesian model to an external pretrained predictive learner. Our method, Predictive-Context LOO (PC-LOO), fits a Bayesian model once, generates posterior-predictive synthetic data, and uses them as context for the learner to approximate each LOO predictive density. Under an oracle construction, averaging contexts generated from the LOO posterior recovers the exact LOO predictive density. Practical PC-LOO replaces the oracle predictor by the pretrained learner and the LOO posterior by the full posterior. We quantify these errors through learner compatibility $C_i$ and deletion influence $V_i$, obtaining $H^2 (p^{LOO}_{i} , q^{PC}_{i}) \leq C_{i} + V_{i}$. The identity also holds for any fixed subset of $M$ context points, allowing $M \ll n - 1$ as the learner's context budget, and scalability to large datasets. Experiments with TabICL support this decomposition, deletion influence decays approximately as $n^{-1}$, while learner compatibility can remain limiting at large $n$. PC-LOO remains accurate with reduced context sizes, while model comparison accuracy is governed by differences in approximation error across models rather than their absolute errors. More broadly, PC-LOO establishes posterior predictive simulation as a principled interface for integrating pretrained predictive learners into Bayesian workflows.

We develop standard errors for estimates of causal effects in settings where the distribution of treatment exposures is known but potentially exhibits complicated dependence patterns across units. Such dependence arises, for example, in settings with spillovers/interference across units or when the treatment assignment is spatially correlated. Our variance estimator generalizes several properties of the Neyman variance estimator for completely randomized experiments to general designs: it is unbiased for an upper bound on the variance from the design-based perspective; the conservativeness of this upper bound depends on the heterogeneity of treatment effects; and it is also valid under certain forms of super-population sampling of the potential outcomes. Our estimator reduces to standard approaches in simple designs, such as the Neyman variance estimator in completely randomized experiments and the cluster-robust variance estimator in cluster-randomized experiments. We show that the proposed estimator is admissible among the class of quadratic-form variance estimators, and derive asymptotic conditions under which it is consistent for an upper bound on the asymptotic variance of the treatment effect estimator. In simulations and applications based on two empirical studies, one with peer effects and another with a spatially correlated treatment, we find that the proposed variance estimator performs well, and in some cases yields standard errors that are 2-3 times shorter than those of existing alternatives.

Generating data with desired properties for a logistic regression model under controlled conditions is challenging. We propose an integrated iterative bisection algorithm that determines values of the intercept and regression coefficients capable of generating data with a specified prevalence and predictive accuracy. The performance of the proposed method is demonstrated through simulation studies. We further show that the choice of parameter values, or, in a Bayesian setting, the choice of priors for these parameters, induces a particular level of predictive accuracy. This relationship should therefore be taken into account when generating data for model evaluation under controlled conditions.

Real decisions are made under incomplete information. If we observe only some of the random variables we need, we can predict the others. The conditional marginals over the missing variables are the key ingredient for computing Bayes risk and Value of Information (VOI), the expected gain from acquiring one more observation before deciding. We present the Marformer, a Transformer trained to directly predict conditional marginals given any set of observed values. Like BERT, which is trained to predict missing words from context, the Marformer constructs a hidden-vector representation for each distribution $p(X_i)$ and iteratively refines it through attention to other distributions $p(X_j)$. Unlike generative approaches, the Marformer does not model the full joint distribution, requires no domain knowledge of the data-generating process, and makes all predictions in a single forward pass. We evaluate across three synthetic domains with missing data—Bayesian networks, discretized multivariate Gaussians, and structured annotation data. The Marformer can match or outperform classical missing-data methods, even when those methods are given the true model family and prior that generated the synthetic data. We also evaluate on a real annotation dataset, where the Marformer outperforms the evaluated baselines at the largest training size. In both cases, the Marformer is substantially faster than the evaluated generative baselines.