Papers

We consider sparsely observed functional data generated by a latent Gauss-Markov process modeled through a linear stochastic differential equation (SDE). Motivated by an application to tumor growth data, we allow for outcome-dependent dropout, where sampling terminates once the most recent observation exceeds a prescribed threshold. Such observation schemes arise naturally in longitudinal studies and violate the fundamental missing completely at random assumption commonly imposed in the analysis of partially observed functional data. We develop a likelihood-based estimation framework that, in contrast to existing moment-based methods, avoids the bias induced by outcome-dependent dropout. We establish convergence rates for the proposed estimators and show that they are mini-max optimal up to logarithmic factors, with the diffusion coefficient admitting a faster rate than the drift coefficients under a random initial condition for the SDE. An application to tumor growth data further demonstrates the advantages of a fully probabilistic functional data model for tasks beyond second-order inference, including prediction bands, first-passage times, and tumor-age estimation.

In high-dimensional tabular learning, feature screening provides a lightweight, model-agnostic way to remove irrelevant features before model fitting. However, scoring raw values directly can miss nonlinear or distributional structure. We introduce PACE (Plug-and-Play Contextual Embedding), which inserts a frozen tabular foundation model (TFM) column encoder before an existing feature-scoring rule, expanding each feature into a higher-dimensional contextual representation. Across controlled studies, PACE improves raw-space screening of complex nonlinear dependence with only modest additional encoding cost. These gains translate to downstream prediction on TALENT datasets: PACE-DC improves binary AUC by 0.077 and multiclass macro-AUC by 0.064, with a median normalized RMSE improvement of 0.063 across ten learners. Matched random-weight and random-feature controls show that PACE gains from pretrained structure beyond generic dimensional expansion. PACE further achieves favorable performance–time trade-offs against task-fitted selectors and attribution-based methods, positioning pretrained column geometry as a reusable upstream primitive for tabular learning.

We study second- and higher-order methods for solving smooth monotone variational inequalities (MVI). Monteiro and Svaiter (SIAM J. Optim., 2012) showed that a second-order method, NPE, converges at a rate of $\mathcal{O}(T^{-1.5})$. For convex-concave minimax optimization, a subclass of MVI problems, Chen, Liu, Luo, and Zhang (COLT 2025) recently improved this rate to $\tilde{\mathcal{O}}( T^{-1.75})$. However, the result has a substantial gap compared to the lower bound of $\Omega(T^{-2.5})$ established by Chen et al. (2026). In this paper, we propose a novel second-order method that achieves the optimal rate of $\mathcal{O}(T^{-2.5})$. Our algorithm also extends to higher-order methods: for any integer $ p \ge 1$, we obtain a $p$th-order method with a convergence rate of $\mathcal{O}(T^{-(3p-1)/2})$, matching the known lower bounds and therefore establishing optimal rates across all orders.

The classical semiparametric convolution theorem characterizes the limiting distributions of regular estimators of pathwise differentiable parameters. We extend this analysis to directionally pathwise differentiable functionals. Regularity along a subspace or convex cone forces the directional derivative to agree there with a bounded linear functional, and every such estimator has a Gaussian convolution factor determined by that functional. We characterize the relevant subspaces, compare their efficiency bounds, and identify the influence functions of efficient estimators. These comparisons show when a stronger regularity requirement increases the variance bound. Under an additivity condition, the Gaussian limit experiment separates into a linear estimation problem and a nonlinear one. Their minimax risks add under squared error; for more general losses, the Gaussian component enters through convolution of the loss. Applications to treatment values, instrumental-variable bounds, and calibration give explicit efficient influence functions, compare existing procedures, and establish feasible attainment in finite-stratum models.

Causal discovery from observational data is a fundamental yet challenging task in scientific research. While existing approaches are primarily based on conditional independence tests, structure scores, or restrictive functional assumptions, we propose Decoupled Causal Discovery (DCD), a novel decoupling-based perspective that does not rely on these methodologies. DCD directly identifies the Markov boundary (MB) by decoupling non-target variables via weighting functions, such that only variables within the MB preserve dependence with the target under the decoupled distribution. Building on this, DCD iteratively constructs the Completed Partially Directed Acyclic Graph (CPDAG) by exploiting structural asymmetries within the MBs. We establish the theoretical identifiability, soundness, and completeness of DCD. Empirical evaluations demonstrate that DCD achieves strong performance, particularly excelling in challenging noise regimes.

In clinical trials, randomization of treatment and control groups is typically used to ensure the groups have similar covariate distributions on average, resulting in unbiased causal effect estimation. Such balanced covariate distributions are hard to achieve in practice, however, due to recruitment costs, patient dropouts, and more. One possible solution to this problem is to include control patients from external, real world databases. In this paper, we propose DBML, a method that selects control units from a real world database by matching the distribution between the treatment group and the potential control group, instead of matching units between two groups. DBML has three steps. First, we detect anomalous database units (with respect to the treatment distribution) using a variational autoencoder. Second, we re-weight the remaining database units to match the distribution of the treatments. Finally, we use these weights to sample units for our control group. The proposed method is compared to alternative matching-based algorithms and weighted algorithms, achieving superior performance in covariate balance.

Evaluating treatment effects in clinical trials with time-to-event outcomes is complicated by intercurrent events (ICEs), e.g., treatment discontinuation or competing events. While the ICH E9 (R1) addendum outlines five strategies to address these challenges, accessible software implementing these methods is limited. This article introduces the \pkg{tteICE} package for \proglang{R}, which implements five strategies to facilitate valid causal inference in the presence of ICEs. \pkg{tteICE} enables researchers to either nonparametrically or semiparametrically efficiently estimate potential cumulative incidence functions under both treatment and control conditions. The treatment effect is assessed by contrasting these functions, with statistical uncertainty quantified using pointwise confidence intervals and $p$-values. The package supports data from randomized controlled trials and observational studies, including settings with competing risks and semi-competing risks. To reduce barriers for applied researchers, \pkg{tteICE} features integrated plotting functions and an accompanying interactive Shiny application that provides a user-friendly graphical interface with step-by-step instructions.

We study exact sign recovery for the PCA connectivity construction of Ding and He (2004), based on the first K-1 principal components. For any fixed number of components with affinely independent means, we derive the oracle signal projection explicitly. Under anisotropic sub-Gaussian noise and random component sizes, a rowwise eigenspace bound yields a sufficient condition, governed by the weakest signal direction, for recovering all connectivity signs without requiring a bounded condition number of the nonzero signal spectrum. Under isotropic Gaussian noise with a fixed mean shape, we obtain an asymptotically sharp cutoff, with an explicit formula when the mixing proportions are equal and the means form a regular simplex.

Bayesian filtering provides a principled framework for online state estimation under uncertainty, yet its application to systems with high-dimensional states and complicated posterior distributions remains challenging. Recent generative models, such as flow matching, have shown potential in Bayesian filtering. However, they still rely on particle-based representations of the posterior, which lose the rich information of the full distribution, or tackle a trajectory-level inverse problem that conflicts with the recursive structure of Bayesian filtering. To address this, we propose a new perspective of directly encoding the evolving distribution into flow matching model weights, namely, the Belief Flow Filter (BFF). It is a generative filtering framework that updates model weights via gradient descent at test time to track the posterior evolution. Thereby, BFF bypasses the scalability issue of particle representations or the flexibility limitation of Gaussian assumptions in conventional filters. We theoretically justify that the BFF design is structurally aligned with Bayesian filtering, and its training objective targets the recursive filtering operator. BFF is empirically verified across 5 different physical systems, including ones with chaotic dynamics and highly sparse, non-linear observations. The results show that BFF attains the best score in 8 of 9 metric-benchmark cells across the three standard 1D and 2D PDE benchmarks, and similarly leads on the extreme single-moving-sensor setting and on a real-world-grounded tokamak plasma estimation task, demonstrating its potential to accurately approximate the Bayesian filtering operator in high-dimensional probability space.

Conformal prediction intervals are used to construct marginally well-calibrated regions around point predictions. Recent work has focused on extending these ideas to conformal predictive distributions, whose marginal probability integral transform (PIT) follows a Unif[0, 1] distribution. In this work, we investigate how to recover predictive densities from conformal predictive distributions. We extend the theory of conformal predictive distributions by proving asymptotic marginal validity for a tail-corrected version of conformal predictive distributions. We propose a method called quantile matching, which preserves an upper bound on the deviation of the marginal PIT from uniformity. Furthermore, we show that the distribution induced by quantile matching is equivalent to the crisp version of conformal predictive distributions when the number of quantiles equals the size of the calibration set. For the recovery of conformal densities, we construct a fidelity-constrained bandwidth optimization for kernel smoothing that preserves asymptotic marginal validity and has a closed-form solution. We apply and compare our proposed methodology on a large simulated real estate transactions dataset based on the Hierarchical Trend Model.

\citet{barczy2014asymptotic} showed that the ordinary least-squares (OLS) estimator of the innovation mean is inconsistent for a unit-root INAR(2) process. We construct a consistent intercept estimator using inverse-time weighted least squares (WLS) and derive its mixed-rate asymptotics: the intercept estimator is asymptotically normal at rate $\sqrt{\log n}$, with a Gaussian limit independent of the autoregressive limits, while the autoregressive estimators retain their OLS rates. We develop an OLS unit-root test calibrated with WLS nuisance estimates. Under the maintained unit root, residual-score Gaussian intervals for the innovation mean and long-run drift remain asymptotically valid after test nonrejection. Simulations show that WLS reduces root mean squared error for both quantities relative to OLS, with further gains from imposing the unit root. An application to two Canadian flood inventories shows that conclusions about persistence depend on the inventory and weighting offset.

We study principal component analysis (PCA) under memory constraints, a setting that is increasingly important in large-scale data analysis. Our focus is on Oja's algorithm, which is a one-pass, memory-efficient algorithm requiring only $O(p)$ storage in the rank-one case and $O(pk)$ storage for $k$- PCA. The main goal is to develop a procedure based on Oja's algorithm that is optimal without prior knowledge of the eigengap, which is otherwise needed to tune the learning rates effectively. We do this by first introducing a geometric perspective on the convergence of $k$-PCA: by embedding the Oja iterates into the exterior space, we establish an exact equivalence between $k$-PCA in the ambient space and $1$-PCA in the exterior space. This geometric viewpoint enables us to establish several convergence guarantees, including one based only on a single learning-rate schedule. Building on this geometric perspective, we develop an averaging theory for $k$-PCA and show that the resulting averaged estimator is adaptive: it achieves the nearly optimal convergence rate without prior knowledge of the eigengap. As an application, we construct memory-efficient estimators for elliptical component analysis (ECA) \cite{han2014scale,han2018eca}. Simulation studies and real data analysis are conducted to demonstrate the benefits of our proposed algorithms.

Online controlled experiments, or A/B tests, are widely used to estimate causal effects on digital platforms. A central challenge is to improve experimental sensitivity, or statistical power, without increasing the experimental sample size. Stratified sampling is a classical variance reduction technique; however, its effectiveness depends critically on how the strata are constructed. We thus propose an optimization-based stratification framework for stratified sampling using optimal multi-way decision trees. Our method, called Optimal Multi-way Stratification Trees (OMST), formulates stratification as a path-selection problem over a feature graph. The selected paths define interpretable stratification rules and are optimized using an exact variance-minimizing binary optimization formulation under continuous proportional allocation and a Neyman-type optimal allocation. We incorporate supervised optimal binning to generate outcome-relevant candidate splits for numerical features. Furthermore, we introduce reduction procedures for redundant candidate paths and assignment constraints, substantially reducing the optimization problem size. Experiments on both a real-world and a simulated dataset demonstrate that OMST achieves comparable or superior variance reduction to existing methods while maintaining shallow and interpretable stratification trees.

Insurance pricing combines pooling with differentiation: a tariff may leave benchmark differences in expected loss partly mutualized or translate them into policy-level premium differences. We propose a benchmark-relative pooling–profiling scale with two complementary coordinates. The coupled $L^p$ coordinate measures policy-level alignment between an evaluated tariff and a stated benchmark pure premium, whereas the marginal Wasserstein coordinate compares their exposure-weighted premium distributions. The difference between their residual $p$-costs defines an allocation mismatch. Under portfolio balance, the coupled $L^1$ coordinate has an exact actuarial interpretation: it is the fraction of the transfer volume induced by full pooling that the tariff removes. Synthetic and motor-insurance applications show that broad classes, proxies, shrinkage and tail caps can affect marginal differentiation, policy-level allocation and transfers differently. A barycentric group-parity intervention further shows that conditional premium disparities can fall mainly through reallocation and restored benchmark-relative transfers, with little change in marginal differentiation.

High-dimensional count data are common in scientific applications, but most diffusion and flow models are designed for continuous or categorical data, and generation often requires many sequential model evaluations. We propose Count Flow Map, a generative model that learns finite-time transitions directly in count space for one- or few-step generation. Our model directly learns stochastic transitions over finite time intervals, using Poisson births and Binomial deaths to preserve nonnegative integer counts without a predefined maximum. These transition models are trained to match the underlying local birth–death dynamics and to maintain consistency across step sizes. We characterize the connection between local dynamics and finite-time transition consistency and derive a bound on the generation error. After validating Count Flow Map in several simulations, including a high-dimensional, high-count setting, we apply it to single-cell drug perturbation prediction and neural population forecasting, where it captures perturbation responses and supports forecasts of high-activity events with only one or a few model evaluations. Together, these experiments demonstrate that Count Flow Map enables high-quality generation directly in count space across inference budgets, from one-step to few-step generation, using a single trained model.

Without interference, cross-fitting enables flexible covariate adjustment while preserving finite-sample unbiasedness under independent unit-level randomization. Under network interference, out-of-sample prediction alone no longer guarantees unbiasedness: assignments entering evaluation-fold Horvitz–Thompson weights may also affect outcomes in the training sample, inducing dependence between fitted predictions and those weights. We develop neighborhood-excluded cross-fitting, which constructs estimand- and design-specific training samples to restore the conditional independence needed for finite-sample unbiasedness without a correctly specified outcome model. We establish asymptotically valid design-based Wald inference for direct and indirect effects under Bernoulli randomization and for the global average treatment effect under Bernoulli cluster randomization. Neighborhood exclusion creates a trade-off in choosing the number of folds: unit-level splitting may require the number of folds to grow with average exclusion-neighborhood size, while cluster-level splitting can substantially relaxes this requirement, permitting a fixed number of folds under partial interference. For linear adjustment under Bernoulli randomization, we derive variance-optimal and confidence-interval-length-optimal procedures, establish explicit rate conditions allowing the covariate dimension to diverge, and show that the variance-optimal procedure is asymptotically no-harm. Simulations illustrate the bias from omitting neighborhood exclusion and the precision gains from adjustment. An application to a social network experiment yields confidence intervals substantially shorter than those from the unadjusted estimator.

Data assimilation aims to infer the state of complex dynamical systems based on observational data. However, accurate inference of the multimodal posteriors induced by nonlinear or non-injective observation operators remains a key challenge under high-dimensional and sparse observation conditions. Ensemble filters scale to high dimensions but are confined by restrictive distributional assumptions, while training-free generative filters (e.g., EnSF, EnFF) alleviate this limitation but may introduce structural errors and hinder information propagation under sparse observations. To address these issues, we propose a training-free, asymptotically exact posterior transport method. Firstly, a covariance-aware posterior transport scheme is designed, which embeds the forecast cross-covariance into flow-based transport and accurately recovers unobserved states while preserving the non-Gaussian posterior structure. Furthermore, the method combines a tractable observation-adaptive proposal with posterior correction, ensuring accurate approximation of the nonlinear posterior distribution. Finally, we establish the corresponding posterior flow theory, from which the asymptotic exactness of the proposed method relative to finite-ensemble surrogates and the Wasserstein error bound are derived. Experiments on Double-Well, Lorenz-96, and Kolmogorov flow show that the proposed method captures complex posterior structure and remains accurate under sparse, nonlinear, and non-injective observations. In the sparse non-injective setting, it reduces RMSE by 56% relative to the best baseline.

Self-calibrating monitors adapt their threshold online to guarantee a prescribed long-run false-alarm rate under arbitrary drift. We compute the price of that guarantee, stating every law with its exact domain of validity. First, the guarantee is an accounting identity, insensitive to what the monitor is meant to detect. Two evidence identities make the cost exact for the online quantile tracker: a persistent step of height $\delta$ yields excess alarm mass within one alarm of $\delta/\eta$, and exactly $\delta/\eta$ pathwise when $\delta$ is a lattice multiple of the gain $\eta$; a ramp of slope $c$ yields a stationary excess rate of exactly $c/\eta$, independent of accumulated size, up to a boundary $c=\eta(1-\alpha)$ coinciding with the alarm-rate cap. Second, the certificate's own fluctuation obeys an exact law: the windowed alarm rate has standard deviation of order $1/L$, not the binomial $1/\sqrt{L}$, since the windowed mass telescopes to a difference of a tight internal state; the closed-form constant is validated with no fitted parameter. Detectors calibrated on the binomial scale are miscalibrated by $\sqrt{\eta\varphi(q_0)L}$, and correct calibration turns detection windows from quadratic to linear in the inverse fault speed. Third, any monitor required to tolerate a drift class $\mathcal{D}$ is blind, at any horizon and for any rule, to every fault in $\mathcal{D}-\mathcal{D}$; the proof is a deliberately elementary two-point argument and the contribution is the object it identifies: for speed-bounded classes the blind set is exactly the doubled-speed class, and the tracker absorbs a speed class fixed by its own gain, so that under a certification regime declaring absorbed drift normal, the monitor manufactures $\mathcal{D}$. An exact Gaussian projection bound, sharper than Pinsker and never vacuous, quantifies power outside it.

Cryptocurrency markets exhibit extreme price swings and sentiment-driven regime shifts, which traditional volatility models often fail to capture. To address this, we develop a partially observed Markov processes (POMP) model augmented with sentiment and heavy-tailed distributions. Specifically, we extend Breto's framework by including the Fear and Greed Index (FGI) as an exogenous regressor in the latent volatility dynamics and replacing Gaussian measurement noise with a Student's t distribution. We fit the model via simulation-based inference using daily Bitcoin returns and FGI data from January 2020 to April 2025. Our approach significantly outperforms three benchmark models in log-likelihood and filter stability. Moreover, our model yields interpretable parameters consistent with known features of financial volatility. These results demonstrate that incorporating sentiment signals and heavy-tailed noise improves the modeling of volatility and regime shifts in cryptocurrency markets.

Bayesian graph alignment estimates correspondence probabilities, but convergence of an alignment-score trace need not imply accurate correspondence marginals. We audit this gap on 240 new exact graph pairs from four source families, 240 larger pairs with 20-100 vertices, and a separate 60-case exact implementation check. Under an explicit edge-flip likelihood, we compare three samplers and score, marginal, indicator, categorical, and classifier-based diagnostics. Marginal disagreement improves error discrimination over score R-hat for the exact informed sampler, but its improvement for vanilla local sampling is uncertain. Assignment-based R* and short indicator panels are competitive; no diagnostic dominates across samplers and endpoints. At larger sizes, diagnostics predict subsequent marginal changes, not posterior error, and classification performance depends on the drift threshold. Disjoint-window and held-out-chain checks attenuate but preserve positive associations. Only 22 of 240 original reference sets pass an agreement screen. On forty failure-selected cases, eightfold SMC particle escalation does not resolve disagreement, whereas additional rejuvenation helps. Longer informed runs remain unstable. An elementary feasible-alignment bound demonstrates severely unrepresentative SMC and informed-chain scores in concentrated 100-vertex cases, independently of approximate reference consensus. We also exhibit common-start chains with near-zero disagreement despite exact marginal error near .967. These results support assignment-sensitive auditing while identifying limits of finite budgets, diagnostic rankings, and reference agreement as evidence of accuracy.

We provide sufficient conditions for the consistency of penalized least squares procedures that select the order (dimension) of a regression model from a sequence of nested classes, allowing for dependent, martingale-difference errors. The main contribution is to relax the classical identifiability requirement: parameters indexing classes larger than the true order need not be identified, provided the additional, excess directions admit a linear approximation to the truth in a neighbourhood of the true parameter. This relaxation lets the number of candidate models grow with the sample size, removing the usual need for a fixed upper bound. We verify the resulting high-level conditions for two classes of nonlinear regression models used in applied work: multiple-regime smooth transition regression and mixture-of-experts models with generic generalized linear model experts. Under a BIC-type penalty, the resulting order selection rule is consistent in either model class whenever the number of candidate models grows slower than the logarithm of the sample size.

Unobserved confounding is a fundamental challenge in causal inference from observational data. This article develops a mixture-learning perspective, viewing latent confounders as sources of heterogeneity that induce mixture structure in observed data. Under suitable structural and identifiability assumptions, recovering the mixing distribution and component mechanisms enables estimation of interventional distributions and causal estimands. Using variants of Bernoulli mixtures as a running example, we contextualize mixture-learning techniques and their structural assumptions, and connect them to causal inference in panel-data settings, including latent factor models and synthetic interventions.We then consider high-dimensional exponential-family mixtures with dependent outcome trajectories, moving beyond counterfactual means to model counterfactual distributions. We situate this perspective relative to complementary approaches for unobserved confounding. Together, these ideas provide a bridge between mixture learning and causal inference, connecting recent advances in high-dimensional mixture learning to scalable identification and estimation of causal effects while raising new challenges for mixture learning.

The Deck-of-cards-based Ordinal Regression (DOR) infers a value function from a ranking of reference alternatives in which the Decision Maker (DM) inserts blank cards between consecutive levels to express preference intensity. DOR, and its stochastic extension (SMAA-DOR), treat these answers as hard constraints defining a set of compatible value functions. We propose B-DOR, a probabilistic reformulation of DOR in which each pair of adjacent levels yields an ordinal observation, the declared direction and the number of cards, modelled through a cumulative-link likelihood that relates the number of blank cards to the latent value difference between alternatives. Two Bayesian inference algorithms are proposed: BAYES-DOR samples the whole posterior distribution by Hamiltonian Monte Carlo; FTRL-DOR tracks the maximum a posteriori estimate by constrained convex optimization. Moreover, through a multi-step elicitation process, elicitation can be spread over several short sessions reducing the cognitive burden on the DM. Both algorithms enjoy logarithmic regret bounds for prediction that hold for any sequence of DM responses and that guide the choice of the prior hyperparameters. A Monte Carlo study over 768 configurations shows that accuracy grows with the number of sessions, that blank cards add significant information over preference directions alone, that both algorithms maintain good performance under inconsistent answers, and that both outperform DOR and SMAA-DOR. An illustrative application to Italian regional healthcare performance demonstrates the practical applicability of the approach for building composite indicators.

Clinical trials in critical care often evaluate outcomes that are truncated by death, such as time to discharge alive, for which treatment effects are not well defined among patients who would not survive under one or both treatment strategies. Principal stratification provides a natural framework for defining survivor causal effects, but existing approaches often rely on monotonicity assumptions that may be implausible in settings where treatment can affect survival in competing directions. This concern is illustrated by the ARDS Network trial of lower versus higher positive end-expiratory pressure (PEEP), in which higher PEEP may benefit some patients by improving alveolar recruitment while harming others through over-distention or hemodynamic compromise. Moreover, the largely null average findings of the original trial do not rule out the possibility of clinically meaningful subgroups that may benefit from or be harmed by higher PEEP. We propose a monotonicity-free framework for identifying conditional survivor average causal effects (CSACE) using an interpretable sensitivity parameter that characterizes latent principal-stratum membership. We then develop a flexible Bayesian estimation strategy based on BART, a posterior-mean-based variable importance measure, and a distribution-aware conditional inference tree that uses the full posterior draws of individualized CSACEs. In simulation studies, the proposed BART-based approach improves individualized CSACE estimation and variable-importance recovery compared with linear modeling, especially under nonlinear treatment-effect heterogeneity. Applied to the ARDS PEEP trial, our method reveals clinically interpretable heterogeneity among estimated always-survivors, identifying subgroups with posterior evidence of benefit or harm that are obscured by the overall null trial results.

Modern prediction-driven decision systems often rely on black-box predictors, but a point forecast alone does not provide the uncertainty scale required for robust downstream decision-making. We build a score-calibrated robustness framework that converts any fixed point predictor into a decision-relevant uncertainty representation through distribution-free conformal calibration. We use the conformal score, rather than a particular uncertainty set, as the primitive unit of robustness. The same score determines coverage-calibrated uncertainty sets for reliability-based robust optimization and normalizes target violations in a target-oriented formulation, Conformal Robust Satisficing. This formulation induces a conformal fragility measure that quantifies how rapidly performance deteriorates as the realized parameter departs from the forecast on the conformal score scale. For objective-uncertainty problems under standard convexity and duality conditions, we show that the reliability-based and target-oriented formulations parameterize the same score-calibrated robust decision frontier. This equivalence yields a data-driven mapping between reliability levels and acceptable targets and characterizes the marginal cost of robustness. Synthetic experiments validate the theoretical guarantees and illustrate the reliability-target correspondence. A real-data online-grocery case study demonstrates how the interface combines deep-learning demand forecasts with tractable inventory optimization, thereby improving reliability and reducing operational costs. Overall, our work shows that conformal scores endow fixed black-box predictors with an interpretable uncertainty scale for downstream decision-making while enabling reliability guarantees, acceptable-target selection, and fragility analysis within a unified framework.

We present an interactive classroom activity designed to address a central challenge in teaching introductory Bayesian statistics: how to formalize subjective knowledge and available information into prior distributions and then update them with empirical data. Role-playing as data analysts for a racing team, students evaluate candidate tires by converting qualitative engineering reports into prior distributions, collecting primary data via a virtual racing game, and using a Beta-Binomial model to inform team strategy. This discovery-based exercise allows small groups to observe directly how different prior choices and sample data jointly shape posterior inference. Student feedback ($n=32$) highlights high enjoyment, engagement and improved conceptual clarity. Open-access materials to implement the activity are provided, alongside recommendations for adapting it to other teaching contexts.

Comparing probability measures in machine learning trades transport geometry against computational cost: Wasserstein distances encode the geometry of $\mathbb R^d$ but require solving a transport problem, while kernel discrepancies are cheap to evaluate yet depend delicately on their test class. We study the Toscani–Fourier family $\mathrm T_{s,p}$, the weighted $L^p$ norm of the difference of two characteristic functions, as a continuous Fourier-side discrepancy on $\mathbb R^d$. For $1\le p<\infty$ we show that $d/p<s<1+d/p$ is exactly the window in which $\mathrm T_{s,p}$ is finite on $\mathcal P_p(\mathbb R^d)$, both endpoints already failing for a pair of Dirac measures, and we establish the metric, embedding, and compactness structure of the resulting space, which we prove to be complete. Duality identifies $\mathrm T_{s,p}$ as an integral probability metric over a homogeneous Fourier–Lebesgue ball, with an explicit extremizer when $1<p<\infty$. We prove the global bound $\mathrm T_{s,p}\lesssim W_p^{\,s-d/p}$, whose exponent is sharp, show that no global converse of any form can hold, and recover topological equivalence with $W_p$ on bounded-support and uniform-tail classes, together with explicit reverse moduli on bounded-support classes that improve the imported energy-kernel exponent at $p=2$. There, $\mathrm T_{s,2}$ is a constant multiple of the classical energy distance, which yields an exact finite-sample identity for the mean of the empirical discrepancy; the numerical experiments are otherwise diagnostic.

Health and well-being indices are widely used to assess population health outcomes and inform policy decisions. Individual-level assessment of well-being can be used to develop community-level indices that measure wellness for different geographical units. While many existing indices operate at coarse geographic levels such as counties or states, finer spatial resolution can offer more actionable insights. We present a novel Bayesian double generalized beta regression framework to model a bounded individual-level well-being index (WBI) using annual survey data collected from 2021 to 2023 in Massachusetts. Although survey respondents may differ across years, responses are geotagged to ZIP Code Tabulation Areas (ZCTAs), enabling the integration of both spatial and temporal information. Our framework incorporates spatial dependencies via a graph Laplacian matrix that encodes driving time-based ZCTA neighborhood structure, and leverages temporal borrowing by using posterior spatial effect estimates from one year to inform priors in the next. This dual-borrowing strategy within a Bayesian double generalized beta regression framework enhances estimation precision, particularly in areas with sparse data, and improves inference for smaller geographic units. We demonstrate the utility of our method through a realistic simulation study that highlights improved estimation when borrowing spatial and temporal information. In the real data analysis, we model individual-level well-being for residents in Massachusetts and find that income, education status, and marital status are most associated with WBI. Additionally, we observe that ZCTAs in Western Massachusetts, Cape Cod, and those near Boston perform best with the highest spatial effects.

Mechanistic epidemic calibration can be sensitive to short-lived reporting anomalies that give a few observations high leverage under least squares (LSQ). We develop direct forward-solver least absolute deviations (LAD) calibration for time-varying susceptible-exposed-infectious-removed models fitted to reporting-interval incidence, with finite-horizon regularity and large-sample results. A phase-aware Monte Carlo study compares LAD and LSQ across three transmission drivers, four epidemic phases, five reporting mechanisms, and four forecast horizons. All 120,000 primary fits converged. LAD was favored in 56.7% of mean-absolute-error comparisons and 57.9% of weighted-interval-score comparisons overall; under isolated spikes, backlog release, or temporary underreporting, these proportions rose to 73.6% and 72.2%, with strong gains near the epidemic peak. Rolling-origin validation on four synthetic RAPIDD Ebola scenarios also identifies settings where LSQ is preferred, demonstrating the framework's ability to distinguish robustness gains from sustained-trend behavior. Multiple uncertainty and identifiability diagnostics, complete code and data, and an interactive application provide a reproducible computational workflow.

In a spiked covariance model, let $\mathcal{H}$ and $\mathcal{B}$ be selected sample and population eigenvector matrices. Their alignment matrix $\Xi=\mathcal{H}^\top\mathcal{B}$ is a latent oracle diagnostic of direction-specific correspondence beyond subspace overlap. We study its uniform pathwise recovery in a nested high-dimensional, low-sample-size experiment with fixed sample size and growing dimension. At every fixed model, an observable spectral contrast consistently estimates the row Gram matrix $\Xi\Xi^\top$. For two or more spikes, however, no measurable estimator sequence achieves uniformly reliable almost-sure path recovery of $\Xi$ over a family having common spike strengths and a common signal subspace but different ordered eigenframes. For every $0<\epsilon<1/\sqrt{2}$, a finite orbit forces each estimator to have some orbit point at which the approximate path-recovery probability is at most a variational bound $u(\epsilon)$, where $u(\epsilon)$ tends to zero as $\epsilon$ tends to zero. The worst-model exact-recovery probability over a countable orbit is zero. Thus the Gram summary is recoverable model by model, while the full alignment matrix is not uniformly recoverable along the nested path.