Papers

Archetypal analysis serves as an exploratory tool that interprets a collection of observations as convex combinations of pure (extreme) patterns. When these patterns correspond to actual observations within the sample, they are termed archetypoids. For the first time, we propose applying archetypoid analysis to nominal observations, specifically for identifying archetypal cases from questionnaires featuring nominal multiple-choice questions with a single possible answer. This approach can enhance our understanding of a nominal data set, similar to its application in multivariate contexts. We compare this methodology with the use of archetype analysis and probabilistic archetypal analysis and demonstrate the benefits of this methodology using a real-world example: the German credit dataset.

As ancient stellar systems, globular clusters (GCs) offer valuable insights into the dynamical histories of large galaxies. Previous studies of GC populations in the inner and outer regions of the Andromeda Galaxy (M31) have revealed intriguing subpopulations with distinct kinematic properties. Here, we build upon earlier studies by employing Bayesian modelling to investigate the kinematics of the combined inner and outer GC populations of M31. Given the heterogeneous nature of the data, we examine subpopulations defined by GCs' metallicity and by associations with substructure, in order to characterise possible relationships between the inner and outer GC populations. We find that lower-metallicity GCs and those linked to substructures exhibit a common, more rapid rotation, whose alignment is distinct from that of higher-metallicity and non-substructure GCs. Furthermore, the higher-metallicity GCs rotate in alignment with Andromeda's stellar disk. These pronounced kinematic differences reinforce the idea that different subgroups of GCs were accreted to M31 at distinct epochs, shedding light on the complex assembly history of the galaxy.

Empirical models of multi-product demand rely on low-dimensional product representations to capture substitution patterns, increasingly using proxies built from unstructured data. When proxies are imperfect, standard workflows yield biased counterfactuals and invalid inference. We develop a practical toolkit to address these issues. Our methods apply to market-level and/or individual data, require minimal additional computation, provide simple standard-error formulas, and accommodate proxies from fine-tuned models. Further, we propose diagnostics to assess proxy quality. Our methods yield meaningful improvements in predicting substitution in empirically calibrated simulations and in an application where we assess counterfactual prediction performance against a ground truth.

Gaussian Process (GP) regression is a powerful nonparametric Bayesian framework, but its performance depends critically on the choice of covariance kernel. Selecting an appropriate kernel is therefore central to model quality, yet remains one of the most challenging and computationally expensive steps in probabilistic modeling. We present a Bayesian optimization framework built on kernel-of-kernels geometry, using expected divergence-based distances between GP priors to explore kernel space efficiently. A multidimensional scaling (MDS) embedding of this distance matrix maps a discrete kernel library into a continuous Euclidean manifold, enabling smooth BO. In this formulation, the input space comprises kernel compositions, the objective is the log marginal likelihood, and featurization is given by the MDS coordinates. When the divergence yields a valid metric, the embedding preserves geometry and produces a stable BO landscape. We demonstrate the approach on synthetic benchmarks, real-world time-series datasets, and an additive manufacturing case study predicting melt-pool geometry, achieving superior predictive accuracy and uncertainty calibration relative to baselines including Large Language Model (LLM)-guided search. This framework establishes a reusable probabilistic geometry for kernel search, with direct relevance to GP modeling and deep kernel learning.

Modern data analysis increasingly requires flexible conditional inference P(X_B | X_A) where (X_A, X_B) is an arbitrary partition of observed variable X. Existing approaches are either restricted to a fixed conditioning structure or depend strongly on the distribution of conditioning masks during training. To address these limitations, we introduce Bayesian generative modeling (BGM), a unified framework for arbitrary conditional inference. BGM learns a generative model of X via a stochastic iterative Bayesian updating algorithm in which model parameters and latent variables are updated until convergence. Once trained, any conditional distribution can be obtained without retraining. Empirically, BGM achieves superior predictive performance with posterior predictive intervals, demonstrating that a single learned model can serve as a universal engine for conditional prediction with principled uncertainty quantification. We provide theoretical guarantees for convergence of the stochastic iterative algorithm, statistical consistency, and conditional risk bounds. The proposed BGM framework leverages modern AI to capture complex relationships among variables while adhering to Bayesian principles, offering a promising approach for a wide range of applications in modern data science. Code for BGM is available at https://github.com/liuq-lab/bayesgm. Document of BGM is available at https://bayesgm.readthedocs.io.

Regression models, where the response variable is circular, are common in areas such as biology, geology and meteorology. A typical model assumes that the conditional distribution of the response follows a von-Mises distribution. However, this assumption is inadequate when the response variable is multimodal. For this reason, in this paper, a finite mixture of regressions model is proposed for the case of a circular response variable and a set of circular and/or linear covariates. Mixture models are very useful when the underlying population is multimodal. Despite the prevalence of multimodality in regression modelling of circular data, the use of mixtures of regressions has received no attention in the literature. This paper aims to close this knowledge gap. To estimate the proposed model, we develop a maximum likelihood estimation procedure via the Expectation-Maximization algorithm. An extensive simulation study is used to demonstrate the practical use and performance of the proposed model and estimation procedure. In addition, the model is shown to be useful as a model-based clustering tool. Lastly, the model is applied to a real dataset from a wind farm in South Africa.

Canonical Polyadic (CP) tensor decomposition is a workhorse algorithm for discovering underlying low-dimensional structure in tensor data. This is accomplished in conventional CP decomposition by fitting a low-rank tensor to data with respect to the least-squares loss. Generalized CP (GCP) decompositions generalize this approach by allowing general loss functions that can be more appropriate, e.g., to model binary and count data or to improve robustness to outliers. However, GCP decompositions do not explicitly account for any symmetry in the tensors, which commonly arises in modern applications. For example, a tensor formed by stacking the adjacency matrices of a dynamic graph over time will naturally exhibit symmetry along the two modes corresponding to the graph nodes. In this paper, we develop a symmetric GCP (SymGCP) decomposition that allows for general forms of symmetry, i.e., symmetry along any subset of the modes. SymGCP accounts for symmetry by enforcing the corresponding symmetry in the decomposition. We derive gradients for SymGCP that enable its efficient computation via all-at-once optimization with existing tensor kernels. The form of the gradients also leads to various stochastic approximations that enable us to develop stochastic SymGCP algorithms that can scale to large tensors. We demonstrate the utility of the proposed SymGCP algorithms with a variety of experiments on both synthetic and real data.

Neuro-symbolic reasoning systems face fundamental challenges in maintaining semantic coherence while satisfying physical and logical constraints. Building upon our previous work on Ontology Neural Networks, we present an enhanced framework that integrates topological conditioning with gradient stabilization mechanisms. The approach employs Forman-Ricci curvature to capture graph topology, Deep Delta Learning for stable rank-one perturbations during constraint projection, and Covariance Matrix Adaptation Evolution Strategy for parameter optimization. Experimental evaluation across multiple problem sizes demonstrates that the method achieves mean energy reduction to 1.15 compared to baseline values of 11.68, with 95 percent success rate in constraint satisfaction tasks. The framework exhibits seed-independent convergence and graceful scaling behavior up to twenty-node problems, suggesting that topological structure can inform gradient-based optimization without sacrificing interpretability or computational efficiency.

Accurate prediction of structural dynamics is imperative for preserving digital twin fidelity throughout operational lifetimes. Parametric models with fixed nominal parameters often omit critical physical effects due to simplifications in geometry, material behavior, damping, or boundary conditions, resulting in model form errors (MFEs) that impair predictive accuracy. This work introduces a comprehensive framework for MFE estimation and correction in high-dimensional finite element (FE) based structural dynamical systems. The Gaussian Process Latent Force Model (GPLFM) represents discrepancies non-parametrically in the reduced modal domain, allowing a flexible data-driven characterization of unmodeled dynamics. A linear Bayesian filtering approach jointly estimates system states and discrepancies, incorporating epistemic and aleatoric uncertainties. To ensure computational tractability, the FE system is projected onto a reduced modal basis, and a mesh-invariant neural network maps modal states to discrepancy estimates, permitting model rectification across different FE discretizations without retraining. Validation is undertaken across five MFE scenarios-including incorrect beam theory, damping misspecification, misspecified boundary condition, unmodeled material nonlinearity, and local damage demonstrating the surrogate model's substantial reduction of displacement and rotation prediction errors under unseen excitations. The proposed methodology offers a potential means to uphold digital twin accuracy amid inherent modeling uncertainties.

We prove tight lower bounds for online multicalibration, establishing an information-theoretic separation from marginal calibration. In the general setting where group functions can depend on both context and the learner's predictions, we prove an $\Omega(T^{2/3})$ lower bound on expected multicalibration error using just three disjoint binary groups. This matches the upper bounds of Noarov et al. (2025) up to logarithmic factors and exceeds the $O(T^{2/3-\varepsilon})$ upper bound for marginal calibration (Dagan et al., 2025), thereby separating the two problems. We then turn to lower bounds for the more difficult case of group functions that may depend on context but not on the learner's predictions. In this case, we establish an $\widetilde{\Omega}(T^{2/3})$ lower bound for online multicalibration via an $O(\log^3 T)$-sized group family constructed from an orthonormal basis, again matching upper bounds up to logarithmic factors.

I propose a novel framework that integrates stochastic differential equations (SDEs) with deep generative models to improve uncertainty quantification in machine learning applications involving structured and temporal data. This approach, termed Stochastic Latent Differential Inference (SLDI), embeds an Itô SDE in the latent space of a variational autoencoder, allowing for flexible, continuous-time modeling of uncertainty while preserving a principled mathematical foundation. The drift and diffusion terms of the SDE are parameterized by neural networks, enabling data-driven inference and generalizing classical time series models to handle irregular sampling and complex dynamic structure. A central theoretical contribution is the co-parameterization of the adjoint state with a dedicated neural network, forming a coupled forward-backward system that captures not only latent evolution but also gradient dynamics. I introduce a pathwise-regularized adjoint loss and analyze variance-reduced gradient flows through the lens of stochastic calculus, offering new tools for improving training stability in deep latent SDEs. My paper unifies and extends variational inference, continuous-time generative modeling, and control-theoretic optimization, providing a rigorous foundation for future developments in stochastic probabilistic machine learning.

One-shot prediction enables rapid adaptation of pretrained foundation models to new tasks using only one labeled example, but lacks principled uncertainty quantification. While conformal prediction provides finite-sample coverage guarantees, standard split conformal methods are inefficient in the one-shot setting due to data splitting and reliance on a single predictor. We propose Conformal Aggregation of One-Shot Predictors (CAOS), a conformal framework that adaptively aggregates multiple one-shot predictors and uses a leave-one-out calibration scheme to fully exploit scarce labeled data. Despite violating classical exchangeability assumptions, we prove that CAOS achieves valid marginal coverage using a monotonicity-based argument. Experiments on one-shot facial landmarking and RAFT text classification tasks show that CAOS produces substantially smaller prediction sets than split conformal baselines while maintaining reliable coverage.

We revisit a fundamental question in hypothesis testing: given two sets of probability measures $\mathcal{P}$ and $\mathcal{Q}$, when does a nontrivial (i.e.\ strictly unbiased) test for $\mathcal{P}$ against $\mathcal{Q}$ exist? Le~Cam showed that, when $\mathcal{P}$ and $\mathcal{Q}$ have a common dominating measure, a test that has power exceeding its level by more than $\varepsilon$ exists if and only if the convex hulls of $\mathcal{P}$ and $\mathcal{Q}$ are separated in total-variation distance by more than $\varepsilon$. The requirement of a dominating measure is frequently violated in nonparametric statistics. In a passing remark, Le~Cam described an approach to address more general scenarios, but he stopped short of stating a formal theorem. This work completes Le~Cam's program, by presenting a matching necessary and sufficient condition for testability: for the aforementioned theorem to hold without assumptions, one must take the closures of the convex hulls of $\mathcal{P}$ and $\mathcal{Q}$ in the space of bounded finitely additive measures. We provide simple elucidating examples, and elaborate on various subtle measure-theoretic and topological points regarding compactness and achievability.

In the advent of big data and machine learning, researchers now have a wealth of congressional candidate ideal point estimates at their disposal for theory testing. Weak relationships raise questions about the extent to which they capture a shared quantity – rather than idiosyncratic, domain-specific factors – yet different measures are used interchangeably in most substantive analyses. Moreover, questions central to the study of American politics implicate relationships between candidate ideal points and other variables derived from the same data sources, introducing endogeneity. We propose a method, consensus multidimensional scaling (CoMDS), which better aligns with how applied scholars use ideal points in practice. CoMDS captures the shared, stable associations of a set of underlying ideal point estimates and can be interpreted as their common spatial representation. We illustrate the utility of our approach for assessing relationships within domains of existing measures and provide a suite of diagnostic tools to aid in practical usage.

In this work, we give a ${\rm poly}(d,k)$ time and sample algorithm for efficiently learning the parameters of a mixture of $k$ spherical distributions in $d$ dimensions. Unlike all previous methods, our techniques apply to heavy-tailed distributions and include examples that do not even have finite covariances. Our method succeeds whenever the cluster distributions have a characteristic function with sufficiently heavy tails. Such distributions include the Laplace distribution but crucially exclude Gaussians. All previous methods for learning mixture models relied implicitly or explicitly on the low-degree moments. Even for the case of Laplace distributions, we prove that any such algorithm must use super-polynomially many samples. Our method thus adds to the short list of techniques that bypass the limitations of the method of moments. Somewhat surprisingly, our algorithm does not require any minimum separation between the cluster means. This is in stark contrast to spherical Gaussian mixtures where a minimum $\ell_2$-separation is provably necessary even information-theoretically [Regev and Vijayaraghavan '17]. Our methods compose well with existing techniques and allow obtaining ”best of both worlds" guarantees for mixtures where every component either has a heavy-tailed characteristic function or has a sub-Gaussian tail with a light-tailed characteristic function. Our algorithm is based on a new approach to learning mixture models via efficient high-dimensional sparse Fourier transforms. We believe that this method will find more applications to statistical estimation. As an example, we give an algorithm for consistent robust mean estimation against noise-oblivious adversaries, a model practically motivated by the literature on multiple hypothesis testing. It was formally proposed in a recent Master's thesis by one of the authors, and has already inspired follow-up works.

Feature selection (FS) is essential for biomarker discovery and clinical predictive modeling. Over the past decades, methodological literature on FS has become rich and mature, offering a wide spectrum of algorithmic approaches. However, much of this methodological progress has not fully translated into applied biomedical research. Moreover, challenges inherent in biomedical data, such as high-dimensional feature space, low sample size, multicollinearity, and missing values, make FS non-trivial. To help bridge this gap between methodological development and practical application, we propose ROOFS (RObust biOmarker Feature Selection), a Python package available at https://gitlab.inria.fr/compo/roofs, designed to help researchers in the choice of FS method adapted to their problem. ROOFS benchmarks multiple FS methods on the user's data and generates reports summarizing a comprehensive set of evaluation metrics, including downstream predictive performance estimated using optimism correction, stability, robustness of individual features, and true positive and false positive rates assessed on semi-synthetic data with a simulated outcome. We demonstrate the utility of ROOFS on data from the PIONeeR clinical trial, aimed at identifying predictors of resistance to anti-PD-(L)1 immunotherapy in lung cancer. Of the 34 FS methods gathered in ROOFS, we evaluated 23 in combination with 11 classifiers (253 models) and identified a filter based on the union of Benjamini-Hochberg false discovery rate-adjusted p-values from t-test and logistic regression as the optimal approach, outperforming other methods including widely used LASSO. We conclude that comprehensive benchmarking with ROOFS has the potential to improve the reproducibility of FS discoveries and increase the translational value of clinical models.

Simulation studies are used to understand the properties of statistical methods. A key luxury in many simulation studies is knowledge of the true value (i.e. the estimand) being targeted. With this oracle knowledge in-hand, the researcher conducting the simulation study can assess across repeated realizations of the data how well a given method recovers the truth. In causal inference simulation studies, the truth is rarely a simple parameter of the statistical model chosen to generate the data. Instead, the estimand is often an average treatment effect, marginalized over the distribution of confounders and/or mediators. Luckily, these variables are often generated from common distributions such as the normal, uniform, exponential, or gamma. For all these distributions, Gaussian quadratures provide efficient and accurate calculation for integrands with integral kernels that stem from known probability density functions. We demonstrate through four applications how to use Gaussian quadrature to accurately and efficiently compute the true causal estimand. We also compare the pros and cons of Gauss-Hermite quadrature to Monte Carlo integration approaches, which we use as benchmarks. Overall, we demonstrate that the Gaussian quadrature is an accurate tool with negligible computation time, yet is underused for calculating the true causal estimands in simulation studies.

A key challenge in network science is the detection of communities, which are sets of nodes in a network that are densely connected internally but sparsely connected to the rest of the network. A fundamental result in community detection is the existence of a nontrivial threshold for community detectability on sparse graphs that are generated by the planted partition model (PPM). Below this so-called “detectability limit”, no community-detection method can perform better than random chance. Spectral methods for community detection fail before this detectability limit because the eigenvalues corresponding to the eigenvectors that are relevant for community detection can be absorbed by the bulk of the spectrum. One can bypass the detectability problem by using special matrices, like the non-backtracking matrix, but this requires one to consider higher-dimensional matrices. In this paper, we show that the difference in graph energy between a PPM and an Erdős–Rényi (ER) network has a distinct transition at the detectability threshold even for the adjacency matrices of the underlying networks. The graph energy is based on the full spectrum of an adjacency matrix, so our result suggests that standard graph matrices still allow one to separate the parameter regions with detectable and undetectable communities.

Building efficient and effective generative models for neural network weights has been a research focus of significant interest that faces challenges posed by the high-dimensional weight spaces of modern neural networks and their symmetries. Several prior generative models are limited to generating partial neural network weights, particularly for larger models, such as ResNet and ViT. Those that do generate complete weights struggle with generation speed or require finetuning of the generated models. In this work, we present DeepWeightFlow, a Flow Matching model that operates directly in weight space to generate diverse and high-accuracy neural network weights for a variety of architectures, neural network sizes, and data modalities. The neural networks generated by DeepWeightFlow do not require fine-tuning to perform well and can scale to large networks. We apply Git Re-Basin and TransFusion for neural network canonicalization in the context of generative weight models to account for the impact of neural network permutation symmetries and to improve generation efficiency for larger model sizes. The generated networks excel at transfer learning, and ensembles of hundreds of neural networks can be generated in minutes, far exceeding the efficiency of diffusion-based methods. DeepWeightFlow models pave the way for more efficient and scalable generation of diverse sets of neural networks.

Drug overdose deaths, including from opioids, remain a significant public health threat to the United States (US). To abate the harms of opioid misuse, understanding its prevalence at the local level is crucial for stakeholders in communities to develop response strategies that effectively use limited resources. Although there exist several state-specific studies that provide county-level prevalence estimates, such estimates are not widely available across the country, as the datasets used in these studies are not always readily available in other states, which, therefore, has limited the wider applications of existing models. To fill this gap, we propose a Bayesian multi-state data integration approach that fully utilizes publicly available data sources to estimate county-level opioid misuse prevalence for all counties in the US. The hierarchical structure jointly models opioid misuse prevalence and overdose death outcomes, leverages existing county-level prevalence estimates in limited states and state-level estimates from national surveys, and accounts for heterogeneity across counties and states with counties' covariates and mixed effects. Furthermore, our parsimonious and generalizable modeling framework employs horseshoe+ prior to flexibly shrink coefficients and prevent overfitting, ensuring adaptability as new county-level prevalence data in additional states become available. Using real-world data, our model shows high estimation accuracy through cross-validation and provides nationwide county-level estimates of opioid misuse for the first time.

Quantitative understanding of stochastic dynamics in limit order price changes is essential for execution strategy design. We analyze intraday transition dynamics of ask and bid orders across market capitalization tiers using high-frequency NASDAQ100 tick data. Employing a discrete-time Markov chain framework, we categorize consecutive price changes into nine states and estimate transition probability matrices (TPMs) for six intraday intervals across High ($\mathtt{HMC}$), Medium ($\mathtt{MMC}$), and Low ($\mathtt{LMC}$) market cap stocks. Element-wise TPM comparison reveals systematic patterns: price inertia peaks during opening and closing hours, stabilizing midday. A capitalization gradient is observed: $\mathtt{HMC}$ stocks exhibit the strongest inertia, while $\mathtt{LMC}$ stocks show lower stability and wider spreads. Markov metrics, including spectral gap, entropy rate, and mean recurrence times, quantify these dynamics. Clustering analysis identifies three distinct temporal phases on the bid side – Opening, Midday, and Closing, and four phases on the ask side by distinguishing Opening, Midday, Pre-Close, and Close. This indicates that sellers initiate end-of-day positioning earlier than buyers. Stationary distributions show limit order dynamics are dominated by neutral and mild price changes. Jensen-Shannon divergence confirms the closing hour as the most distinct phase, with capitalization modulating temporal contrasts and bid-ask asymmetry. These findings support capitalization-aware and time-adaptive execution algorithms.

We show how to construct the implied copula process of response values from a Bayesian additive regression tree (BART) model with prior on the leaf node variances. This copula process, defined on the covariate space, can be combined with any marginal distribution for the dependent variable to construct a flexible distributional BART model. Bayesian inference is performed via Markov chain Monte Carlo on an augmented posterior, where we show that key sampling steps can be realized as those of Chipman et al. (2010), preserving scalability and computational efficiency despite the high dimensionality of the copula process. The posterior predictive distribution from the copula process model is derived in closed form as the push-forward of the posterior predictive distribution of the underlying BART model with an optimal transport map. Under suitable conditions, and allowing for model misspecification, we show that the model-implied conditional density closest in Kullback-Leibler divergence to the true data-generating process exists, is unique and has a closed form expression. We further establish that fractional posteriors for the conditional densities at the observed covariate values concentrate around this projection at a rate that adapts to the smoothness of the corresponding regression function and to the number of relevant covariates. Simulation studies demonstrate improved accuracy of distributional predictions compared to the original BART model and leading benchmarks. Applications to five real datasets with 506 to 515,345 observations and 8 to 90 covariates further highlight the efficacy and scalability of our proposed BART copula process model.

We propose an estimator of a concave cumulative distribution function under the measurement error model, where the non-negative variables of interest are perturbed by additive independent random noise. The estimator is defined as the least concave majorant on the positive half-line of the deconvolution estimator of the distribution function. We show its uniform consistency and its square root convergence in law in $\ell_\infty(\mathbb R)$. To assess the validity of the concavity assumption, we construct a test for the nonparametric null hypothesis that the distribution function is concave on the positive half-line, against the alternative that it is not. We calibrate the test using bootstrap methods. The theoretical justification for calibration led us to establish a bootstrap version of Theorem 1 in Söhl and Trabs (2012), a Donsker-type result from which we obtain, as a special case, the limiting behavior of the deconvolution estimator of the distribution function in a bootstrap setting with measurement error. Combining this Donsker-type theorem with the functional delta method, we show that the test statistic and its bootstrap version have the same limiting distribution under the null hypothesis, whereas under the alternative, the bootstrap statistic is stochastically smaller. Consequently, the power of the test tends to one, for any fixed alternative, as the sample size tends to infinity. In addition to the theoretical results for the estimator and the test, we investigate their finite-sample performance in simulation studies.

The aim of this paper is to use Maximum Likelihood (ML) Classification on multispectral data by means of qualitative and quantitative approaches. Maximum Likelihood is a supervised classification algorithm which is based on the Classical Bayes theorem. It makes use of a discriminant function to assign pixel to the class with the highest likelihood. Class means vector and covariance matrix are the key inputs to the function and can be estimated from training pixels of a particular class. As Maximum Likelihood need some assumptions before it has to be applied on the data. In this paper we will compare the results of Maximum Likelihood Classification (ML) before apply the Weierstrass Transform and apply Weierstrass Transform and will see the difference between the accuracy on training pixels of high resolution Quickbird satellite image. Principle Component analysis (PCA) is also used for dimension reduction and also used to check the variation in bands. The results shows that the separation between mean of the classes in the decision space is to be the main factor that leads to the high classification accuracy of Maximum Likelihood (ML) after using Weierstrass Transform than without using it.

Estimating causal effects from observational data remains a fundamental challenge in causal inference, especially in the presence of latent confounders. This paper focuses on estimating causal effects in Gaussian Linear Structural Causal Models (GL-SCMs), which are widely used due to their analytical tractability. However, parameter estimation in GL-SCMs is often infeasible with finite data, primarily due to overparameterization. To address this, we introduce the class of Centralized Gaussian Linear SCMs (CGL-SCMs), a simplified yet expressive subclass where exogenous variables follow standardized distributions. We show that CGL-SCMs are equally expressive in terms of causal effect identifiability from observational distributions and present a novel EM-based estimation algorithm that can learn CGL-SCM parameters and estimate identifiable causal effects from finite observational samples. Our theoretical analysis is validated through experiments on synthetic data and benchmark causal graphs, demonstrating that the learned models accurately recover causal distributions.

This paper considers estimation and model selection of quantile vector autoregression (QVAR). Conventional quantile regression often yields undesirable crossing quantile curves, violating the monotonicity of quantiles. To address this issue, we propose a simplex quantile vector autoregression (SQVAR) framework, which transforms the autoregressive (AR) structure of the original QVAR model into a simplex, ensuring that the estimated quantile curves remain monotonic across all quantile levels. In addition, we impose the smoothly clipped absolute deviation (SCAD) penalty on the SQVAR model to mitigate the explosive nature of the parameter space. We further develop a Bayesian information criterion (BIC)-based procedure for selecting the optimal penalty parameter and introduce new frameworks for impulse response analysis of QVAR models. Finally, we establish asymptotic properties of the proposed method, including the convergence rate and asymptotic normality of the estimator, the consistency of AR order selection, and the validity of the BIC-based penalty selection. For illustration, we apply the proposed method to U.S. stock market data, highlighting the usefulness of our SQVAR method.

This study presents a cluster-based Bayesian SIRD model to analyze the epidemiology of chickenpox (varicella) in India, utilizing data from 1990 to 2021. We employed an age-structured approach, dividing the population into juvenile, adult, and elderly groups, to capture the disease's transmission dynamics across diverse demographic groups. The model incorporates a Holling-type incidence function, which accounts for the saturation effect of transmission at high prevalence levels, and applies Bayesian inference to estimate key epidemiological parameters, including transmission rates, recovery rates, and mortality rates. The study further explores cluster analysis to identify regional clusters within India based on the similarities in chickenpox transmission dynamics, using criteria like incidence, prevalence, and mortality rates. We perform K-means clustering to uncover three distinct epidemiological regimes, which vary in terms of outbreak potential and age-specific dynamics. The findings highlight juveniles as the primary drivers of transmission, while the elderly face a disproportionately high mortality burden. Our results underscore the importance of age-targeted interventions and suggest that regional heterogeneity should be considered in public health strategies for disease control. The model offers a transparent, reproducible framework for understanding long-term transmission dynamics and supports evidence-based planning for chickenpox control in India. The practical utility of the model is further validated through a simulation study.

Fine suspended particulates (FSP), commonly known as PM2.5, are among the most harmful air pollutants, posing serious risks to population health and environmental integrity. As such, accurately identifying clusters of air-quality monitoring stations exhibiting similar FSP levels is essential for guiding targeted interventions and localized public-health response measures. This task, however, is notably nontrivial as FSP levels may depend on various regional and temporal factors, which should be incorporated in the modeling process. Thus, we capitalize on Bayesian nonparametric dynamic clustering ideas, in which clustering structures may be influenced by complex dependencies. Existing implementations of dynamic clustering, however, rely on copula-based dependent Dirichlet processes (DPs), presenting considerable computational challenges for real-world deployment. With this in mind, we propose a more efficient alternative for dynamic clustering by incorporating the novel ideas of logistic-beta dependent DPs. We also adopt a Stirling-gamma prior–a novel distribution family–on the concentration parameter of our underlying DP, easing the process of incorporating prior knowledge into the model. Efficient computational strategies for posterior inference are also presented. We apply our proposed method to identify dynamic clusters of air-quality monitoring stations across Chile and demonstrate its superior performance over existing approaches.

A forecasting model compresses its predictors into an estimate of a conditional mean, and the systematic structure that estimate omits survives in the second moment of its forecast errors. Accuracy comparisons do not measure this structure, and variance-based extraction does not recover the part of it that a given decision bears. In this paper, we propose a distributionally robust framework that recovers the omitted structure from the residuals of a fixed forecaster: a decision is made robust over a two-layer moment ambiguity set on the standardized residual cross-section, and the discovery statistic is the covariance forcing, the component of the decision's residual risk transverse to its exposure. We demonstrate that the forcing is invariant to shrinkage and to isotropic inflation of the covariance, so the recovered direction is a property of the residuals rather than of the regularization, the sense in which the recovery is ground truth. This is confirmed on monthly U.S. Treasury zero-coupon yields from 2006 to 2025, where the recovered factor is named by factor-adjusted robust selection against a panel of 111 macroeconomic, Treasury supply-and-demand, and financial indicators. From the residuals of the linear factor-augmented dynamic Nelson-Siegel benchmark the factor names as a leading business-cycle factor, anchored on the Conference Board leading index and certified by a block-permutation test; from those of the more accurate nonlinear random forest benchmark the same procedure selects the same real-activity family without certification. A neutralization test completes the evidence: removing the recovered factor from the deployed duration position leaves volatility essentially unchanged and worsens the tail, so the factor is a material systematic risk the position bears.

We study eigenvalue fluctuations of kernel and adjacency matrices in dense graphon-based random graph models, for modes associated with simple, isolated nonzero eigenvalues of the graphon integral operator. Our results require no smoothness or Lipschitz continuity. They include degenerate non-leading modes whose normalized eigenfunctions have constant modulus and change sign between $-1$ and $+1$ on measurable sets. Continuous graphons exclude this possibility: their nonzero eigenfunctions have continuous representatives, so constant modulus forces a constant eigenfunction, which corresponds to the largest eigenvalue of a nonnegative graphon. Nevertheless, non-leading degeneracy already occurs at the second eigenvalue of a two-block model with equal blocks. The limiting behavior is governed by the variance of the squared population eigenfunction. When this variance is positive, the empirical eigenvalue, divided by $n-1$ and centered at its population value, has a Gaussian limit under $\sqrt{n}$ scaling. When the variance vanishes, the first-order fluctuation disappears, and the centered unnormalized kernel eigenvalue converges to a weighted series of centered chi-square variables. For adjacency matrices, Bernoulli edge noise is negligible at the $\sqrt{n}$ scale of the unnormalized eigenvalues. In the degenerate regime, after the same deterministic centering used for the kernel eigenvalue, the adjacency eigenvalue converges to the same weighted chi-square series plus a deterministic shift and a Gaussian term independent of that series. The proofs combine spectral localization, exact Schur identities, finite-mode approximations, and Hoeffding decompositions.