Papers

This paper investigates fractional Riesz-Bessel equations with random initial conditions. The spectra of these random initial conditions exhibit singularities both at zero frequency and at non-zero frequencies, which correspond to the cases of classical long-range dependence and cyclic long-range dependence, respectively. Using spectral methods and asymptotic theory, it is shown that the rescaled solutions of the equations converge to spatio-temporal Gaussian random fields. The limit fields are stationary in space and non-stationary in time. The covariance and spectral structures of the resulting asymptotic random fields are provided. The paper further establishes multiscaling limit theorems for the case of regularly varying asymptotics. A numerical example illustrating the theoretical results is also presented.

The scale mixture of the exponential distribution provides a flexible framework for modelling lifetime and reliability data. This model is widely used in survival analysis, biomedical studies, statistical finance, and other related disciplines. In this work, we investigate the estimation of the ordered scale parameter of two scale mixture of exponential distributions under Stein loss and symmetric loss functions. Under certain conditions, we prove the inadmissibility of the affine equivariant estimator and exhibit several improved estimators. Consequently, we propose a class of estimators that uniformly dominate the best affine equivariant estimators (BAEE). Furthermore, we have proved that the boundary estimator of this class is a generalized Bayes estimator. As an application, we have proposed improved estimators for the ordered scale parameters of the multivariate Lomax and exponential inverse Gaussian distributions. For each case, we have conducted a simulation study to compare the risk performance of the improved estimators. Finally, we have given two real-life data analysis for implementation purposes.

Consider the additive Gaussian model $Y = X + \sigma Z$, where $X \sim P$ is an unknown signal, $Z \sim N(0,1)$ is independent of $X$, and $\sigma > 0$ is known. Let $Q$ denote the law of $Y$. We construct a hierarchy of denoisers $T_0, T_1, \ldots, T_\infty \colon \mathbb{R} \to \mathbb{R}$ that depend only on higher-order score functions $q^{(m)}/q$, $m \geq 1$, of $Q$ and require no knowledge of the law $P$. The $K$-th order denoiser $T_K$ involves scores up to order $2K{-}1$ and satisfies $W_r(T_K \sharp Q, P) = O(\sigma^{2(K+1)})$ for every $r \geq 1$; in the limit, $T_\infty$ recovers the monotone optimal transport map (Brenier map) pushing $Q$ onto $P$. We provide a complete characterization of the combinatorial structure governing this hierarchy through partial Bell polynomial recursions, making precise how higher-order score functions encode the Brenier map. We further establish rates of convergence for estimating these scores from $n$ i.i.d.\ draws from $Q$ under two complementary strategies: (i) plug-in kernel density estimation, and (ii) higher-order score matching. The construction reveals a precise interplay among higher-order Fisher-type information, optimal transport, and the combinatorics of integer partitions.

Positional encoding (PE) is a core architectural component of Transformers, yet its impact on the Transformer's generalization and robustness remains unclear. In this work, we provide the first generalization analysis for a single-layer Transformer under in-context regression that explicitly accounts for a completely trainable PE module. Our result shows that PE systematically enlarges the generalization gap. Extending to the adversarial setting, we derive the adversarial Rademacher generalization bound. We find that the gap between models with and without PE is magnified under attack, demonstrating that PE amplifies the vulnerability of models. Our bounds are empirically validated by a simulation study. Together, this work establishes a new framework for understanding the clean and adversarial generalization in ICL with PE.

This paper addresses a long-standing open problem in the analysis of linear mixed models with crossed random effects under unbalanced designs: how to find an analytic expression for the inverse of $\mathbf{V}$, the covariance matrix of the observed response. The inverse matrix $\mathbf{V}^{-1}$ is required for likelihood-based estimation and inference. However, for unbalanced crossed designs, $\mathbf{V}$ is dense and the lack of a closed-form representation for $\mathbf{V}^{-1}$, until now, has made using likelihood-based methods computationally challenging and difficult to analyse mathematically. We use the Khatri–Rao product to represent $\mathbf{V}$ and then to construct a modified covariance matrix whose inverse admits an exact spectral decomposition. Building on this construction, we obtain an elegant and simple approximation to $\mathbf{V}^{-1}$ for asymptotic unbalanced designs. For non-asymptotic settings, we derive an accurate and interpretable approximation under mildly unbalanced data and establish an exact inverse representation as a low-rank correction to this approximation, applicable to arbitrary degrees of unbalance. Simulation studies demonstrate the accuracy, stability, and computational tractability of the proposed framework.

We examine the non-asymptotic properties of robust density ratio estimation (DRE) in contaminated settings. Weighted DRE is the most promising among existing methods, exhibiting doubly strong robustness from an asymptotic perspective. This study demonstrates that Weighted DRE achieves sparse consistency even under heavy contamination within a non-asymptotic framework. This method addresses two significant challenges in density ratio estimation and robust estimation. For density ratio estimation, we provide the non-asymptotic properties of estimating unbounded density ratios under the assumption that the weighted density ratio function is bounded. For robust estimation, we introduce a non-asymptotic framework for doubly strong robustness under heavy contamination, assuming that at least one of the following conditions holds: (i) contamination ratios are small, and (ii) outliers have small weighted values. This work provides the first non-asymptotic analysis of strong robustness under heavy contamination.

Understanding how vaccines perform against different pathogen genotypes is crucial for developing effective prevention strategies, particularly for highly genetically diverse pathogens like HIV. Sieve analysis is a statistical framework used to determine whether a vaccine selectively prevents acquisition of certain genotypes while allowing breakthrough of other genotypes that evade immune responses. Traditionally, these analyses are conducted with a single sequence available per individual acquiring the pathogen. However, modern sequencing technology can provide detailed characterization of intra-individual viral diversity by capturing up to hundreds of pathogen sequences per person. In this work, we introduce methodology that extends sieve analysis to account for intra-individual viral diversity. Our approach estimates vaccine efficacy against viral populations with varying true (unobservable) frequencies of vaccine-mismatched mutations. To account for differential resolution of information from differing sequence counts per person, we use competing risks Cox regression with modeled causes of failure and propose an empirical Bayes approach for the classification model. Simulation studies demonstrate that our approach reduces bias, provides nominal confidence interval coverage, and improves statistical power compared to conventional methods. We apply our method to the HVTN 705 Imbokodo trial, which assessed the efficacy of a heterologous vaccine regimen in preventing HIV-1 acquisition.

Batch effects pose a significant challenge in the analysis of single-cell omics data, introducing technical artifacts that confound biological signals. While various computational methods have achieved empirical success in correcting these effects, they lack the formal theoretical guarantees required to assess their reliability and generalization. To bridge this gap, we introduce Mixture-Model-based Data Harmonization (MoDaH), a principled batch correction algorithm grounded in a rigorous statistical framework. Under a new Gaussian-mixture-model with explicit parametrization of batch effects, we establish the minimax optimal error rates for batch correction and prove that MoDaH achieves this rate by leveraging the recent theoretical advances in clustering data from anisotropic Gaussian mixtures. This constitutes, to the best of our knowledge, the first theoretical guarantee for batch correction. Extensive experiments on diverse single-cell RNA-seq and spatial proteomics datasets demonstrate that MoDaH not only attains theoretical optimality but also achieves empirical performance comparable to or even surpassing those of state-of-the-art heuristics (e.g., Harmony, Seurat-V5, and LIGER), effectively balancing the removal of technical noise with the conservation of biological signal.

There has been significant progress in Bayesian inference based on sparsity-inducing (e.g., spike-and-slab and horseshoe-type) priors for high-dimensional regression models. The resulting posteriors, however, in general do not possess desirable frequentist properties, and the credible sets thus cannot serve as valid confidence sets even asymptotically. We introduce a novel debiasing approach that corrects the bias for the entire Bayesian posterior distribution. We establish a new Bernstein-von Mises theorem that guarantees the frequentist validity of the debiased posterior. We demonstrate the practical performance of our proposal through Monte Carlo simulations and two empirical applications in economics.

Semi-continuous exposure variables, featuring a point mass at zero and a continuous dose in exposed individuals, arise frequently in epidemiology. In such settings, investigators may wish to distinguish the causal effect of exposure at a fixed dose with no exposure, from the dose-response effect among exposed individuals. We develop a two-stage framework for causal inference with semi-continuous exposures and multiple outcomes. The proposed approach yields outcome-specific estimators for these two causal components, establishes their large sample properties, and provides Wald-type tests for homogeneity across outcomes. We further extend the framework to the multi-cohort setting through hierarchical homogeneity testing and evidence synthesis, enabling principled pooling across cohorts and outcomes. We apply the proposed methods to several cohort studies designed to evaluate the effects of prenatal alcohol exposure on executive function, mathematics achievement, and learning and memory. Simulation studies show that the proposed estimators and testing procedures perform well. The framework provides a principled basis for estimation, comparison, and synthesis of causal effects with multiple outcomes and semi-continuous exposures.

Revisiting the continuous-time Mean-Variance (MV) Portfolio Optimization problem, we model the market dynamics with a jump-diffusion process and apply Reinforcement Learning (RL) techniques to facilitate informed exploration within the control space. We recognize the time-inconsistency of the MV problem and adopt the time-inconsistent control (TIC) approach to analytically solve for an exploratory equilibrium investment policy, which is a Gaussian distribution centered on the equilibrium control of the classical MV problem. Our approach accounts for time-inconsistent preferences and actions, and our equilibrium policy is the best option an investor can take at any given time during the investment period. Moreover, we leverage the martingale properties of the equilibrium policy, design a RL model, and propose an Actor-Critic RL algorithm. All of our RL model parameters converge to the corresponding true values in a simulation study. Our numerical study on 24 years of real market data shows that the proposed RL model is profitable in 13 out of 14 tests, demonstrating its practical applicability in real world investment.

Background: Alzheimer's Disease and Related Dementias (ADRD) affects millions worldwide. Significant disparities exist in ADRD diagnosis and care, disproportionately impacting minority and socioeconomically vulnerable populations Objective: In this study, we investigate the relationship between ADRD density and accessibility to healthcare. We identify underserved and overserved areas in Maryland based on diagnosed cases and mortality due to ADRD, focusing on geographic disparities in care. Methods: 2023 Maryland ADRD patients were identified using ICD-10 codes from. Accessibility was measured using the Kernel Density Two-Step Floating Catchment Area (KD2SFCA) method. The Gini index and t-tests were used to analyze disparities between urban and rural areas. Hot Spot Analysis Getis-Ord Gi* and local bivariate relationships analysis were applied to assess spatial correlations. Principal component analysis (PCA) was applied to calculate the health risk index. Results: Hospital accessibility was unevenly distributed. Mortality rates from ADRD were higher in underserved areas with fewer hospitals. Hot spot analysis shows eastern and southern Maryland have zones with high mortality per population and per ADRD patient, surrounded by similarly high-rate zones. Central Maryland shows lower death rates per patient but more hospital facilities. In eastern Maryland, higher poverty areas are surrounded by zones with lower accessibility and higher health risk indices. Conclusion: Hospital accessibility is unevenly distributed, creating major rural disparities. Underserved regions in terms of access to healthcare facilities, particularly in eastern and southern Maryland, exhibit high ADRD mortality rates despite low diagnosis rates. This suggests that many ADRD cases remain undiagnosed, underdiagnosed, or subject to delayed treatment.

Machine learning (ML) is transforming healthcare, but safe clinical decisions demand reliable uncertainty estimates that standard ML models fail to provide. Conformal prediction (CP) is a popular tool that allows users to turn heuristic uncertainty estimates into uncertainty estimates with statistical guarantees. CP works by converting predictions of a ML model, together with a calibration sample, into prediction sets that are guaranteed to contain the true label with any desired probability. An often cited advantage is that CP theory holds for calibration samples of arbitrary size, suggesting that uncertainty estimates with practically meaningful statistical guarantees can be achieved even if only small calibration sets are available. We question this promise by showing that, although the statistical guarantees hold for calibration sets of arbitrary size, the practical utility of these guarantees does highly depend on the size of the calibration set. This observation is relevant in medical domains because data is often scarce and obtaining large calibration sets is therefore infeasible. We corroborate our critique in an empirical demonstration on a medical image classification task.

Do Ethereum's Layer-2 (L2) rollups actually decongest the Layer-1 (L1) mainnet once protocol upgrades and demand are held constant? Using a 1245-day daily panel from August 5, 2021 to December 31, 2024 that spans the London, Merge, and Dencun upgrades, we link Ethereum fee and congestion metrics to L2 user activity, macro-demand proxies, and targeted event indicators. We estimate a regime-aware error-correction model that treats posting-clean L2 user share as a continuous treatment. Over the pre-Dencun (London+Merge) window, a 10 percentage point increase in L2 adoption lowers median base fees by about 13% – roughly 5 Gwei at pre-Dencun levels – and deviations from the long-run relation decay with an 11-day half-life. Block utilization and a scarcity index show similar congestion relief. After Dencun, L2 adoption is already high and treatment support narrows, so blob-era estimates are statistically imprecise and we treat them as exploratory. The pre-Dencun window therefore delivers the first cross-regime causal estimate of how aggregate L2 adoption decongests Ethereum, together with a reusable template for monitoring rollup-centric scaling strategies.

We read with interest the above article by Zavorsky (2025, Respiratory Medicine, doi:10.1016/j.rmed.2024.107836) concerning reference equations for pulmonary function testing. The author compares a Generalized Additive Model for Location, Scale, and Shape (GAMLSS), which is the standard adopted by the Global Lung Function Initiative (GLI), with a segmented linear regression (SLR) model, for pulmonary function variables. The author presents an interesting comparison; however there are some fundamental issues with the approach. We welcome this opportunity for discussion of the issues that it raises. The author's contention is that (1) SLR provides "prediction accuracies on par with GAMLSS"; and (2) the GAMLSS model equations are "complicated and require supplementary spline tables", whereas the SLR is "more straightforward, parsimonious, and accessible to a broader audience". We respectfully disagree with both of these points.

The Weibull distribution is a commonly adopted choice for modeling the survival of systems subject to maintenance over time. When only proxy indicators and censored observations are available, it becomes necessary to express the distribution's parameters as functions of time-dependent covariates. Deep neural networks provide the flexibility needed to learn complex relationships between these covariates and operational lifetime, thereby extending the capabilities of traditional regression-based models. Motivated by the analysis of a fleet of military vehicles operating in highly variable and demanding environments, as well as by the limitations observed in existing methodologies, this paper introduces WTNN, a new neural network-based modeling framework specifically designed for Weibull survival studies. The proposed architecture is specifically designed to incorporate qualitative prior knowledge regarding the most influential covariates, in a manner consistent with the shape and structure of the Weibull distribution. Through numerical experiments, we show that this approach can be reliably trained on proxy and right-censored data, and is capable of producing robust and interpretable survival predictions that can improve existing approaches.

This article presents an Integral Gaussian Process (IntegralGP) framework for volumetric estimation of subterranean properties in mineral deposits. It provides a unified representation for data with different spatial supports, which enables blasthole geochemical assays to be properly modelled as interval observations rather than points. This approach is shown to improve regression performance and boundary delineation. A core contribution is a description of the mathematical changes to the covariance expressions which allow these benefits to be realised. The gradient and anti-derivatives are obtained to facilitate learning of the kernel hyperparameters. Numerical stability issues are also discussed. To illustrate its application, an IntegralGP data fusion algorithm is described. The objective is to assimilate line-based blasthole assays and update a block model that provides long-range prediction of Fe concentration beneath the drilled bench. Heteroscedastic GP is used to fuse chemically compatible but spatially incongruous data with different resolutions and sample spacings. Domain knowledge embodied in the structure and empirical distribution of the block model must be generally preserved while local inaccuracies are corrected. Using validation measurements within the predicted bench, our experiments demonstrate an improvement in bench-below grade prediction performance. For material classification, IntegralGP fusion reduces the absolute error and model bias in categorical prediction, especially instances where waste blocks are mistakenly classified as high-grade.

This article contributes to the discussion on the relationship between the Neyman-Rubin and the graphical frameworks for causal inference. We present specific examples of data-generating mechanisms - such as those involving undirected or deterministic relationships and cycles - where analyses using a directed acyclic graph are challenging, but where the tools from the Neyman-Rubin causal framework are readily applicable. We also provide examples of data-generating mechanisms with M-bias, trapdoor variables, and complex front-door structures, where the application of the Neyman-Rubin approach is complicated, but the graphical approach is directly usable. The examples offer insights into commonly used causal inference frameworks and aim to improve comprehension of the languages for causal reasoning among a broad audience.

Deep learning models for medical image segmentation suffer significant performance drops due to distribution shifts, but the causal mechanisms behind these drops remain poorly understood. We extend causal attribution frameworks to high-dimensional segmentation tasks, quantifying how acquisition protocols and annotation variability independently contribute to performance degradation. We model the data-generating process through a causal graph and employ Shapley values to fairly attribute performance changes to individual mechanisms. Our framework addresses unique challenges in medical imaging: high-dimensional outputs, limited samples, and complex mechanism interactions. Validation on multiple sclerosis (MS) lesion segmentation across 4 centers and 7 annotators reveals context-dependent failure modes: annotation protocol shifts dominate when crossing annotators (7.4% $\pm$ 8.9% DSC attribution), while acquisition shifts dominate when crossing imaging centers (6.5% $\pm$ 9.1%). This mechanism-specific quantification enables practitioners to prioritize targeted interventions based on deployment context.

Accurate forecasting of urban air pollution is essential for protecting public health and guiding mitigation policies. While Deep Learning (DL) and hybrid pipelines dominate recent research, their complexity and limited interpretability hinder operational use. This study investigates whether lightweight additive models – Facebook Prophet (FBP) and NeuralProphet (NP) – can deliver competitive forecasts for particulate matter (PM$_{2.5}$, PM$_{10}$) in Beijing, China. Using multi-year pollutant and meteorological data, we applied systematic feature selection (correlation, mutual information, mRMR), leakage-safe scaling, and chronological data splits. Both models were trained with pollutant and precursor regressors, with NP additionally leveraging lagged dependencies. For context, two machine learning baselines (LSTM, LightGBM) and one traditional statistical model (SARIMAX) were also implemented. Performance was evaluated on a 7-day holdout using MAE, RMSE, and $R^2$. Results show that FBP consistently outperformed NP, SARIMAX, and the learning-based baselines, achieving test $R^2$ above 0.94 for both pollutants. These findings demonstrate that interpretable additive models remain competitive with both traditional and complex approaches, offering a practical balance of accuracy, transparency, and ease of deployment.

Accurate and efficient surrogate modeling is essential for modern computational science, and there are a staggering number of emulation methods to choose from. With new methods being developed all the time, comparing the relative strengths and weaknesses of different methods remains a challenge due to inconsistent benchmarking practices and (sometimes) limited reproducibility and transparency. In this work, we present a large-scale, fully reproducible comparison of $29$ distinct emulators across $60$ canonical test functions and $40$ real emulation datasets. To facilitate rigorous, apples-to-apples comparisons, we introduce the R package duqling, which streamlines reproducible simulation studies using a consistent, simple syntax, and automatic internal scaling of inputs. This framework allows researchers to compare emulators in a unified environment and makes it possible to replicate or extend previous studies with minimal effort, even across different publications. Our results provide detailed empirical insight into the strengths and weaknesses of state-of-the-art emulators and offer guidance for both method developers and practitioners selecting a surrogate for new data. We discuss best practices for emulator comparison and highlight how duqling can accelerate research in emulator design and application.

Stochastic optimization of engineering systems is often infeasible due to repeated evaluations of a computationally expensive, high-fidelity simulation. Bi-fidelity methods mitigate this challenge by leveraging a cheaper, approximate model to accelerate convergence. Most existing bi-fidelity approaches, however, exploit either design-space curvature or random-space correlation, not both. We present BISTRO - a BI-fidelity Stochastic Trust-Region Optimizer for unconstrained optimization under uncertainty through a stochastic approximation procedure. This approach exploits the curvature information of a low-fidelity objective function to converge within a basin of a local minimum of the high-fidelity model where low-fidelity curvature information is no longer valuable. The method then switches to a variance-reduced stochastic gradient descent procedure. We provide convergence guarantees in expectation under certain regularity assumptions and ensure the best-case $\mathcal{O}(1/n)$ convergence rate for stochastic optimization. On benchmark problems and a 20-dimensional space shuttle reentry case, BISTRO converges faster than adaptive sampling and variance reduction procedures and cuts computational expense by up to 29x.

Kernel density estimation is a key component of a wide variety of algorithms in machine learning, Bayesian inference, stochastic dynamics and signal processing. However, the unsupervised density estimation technique requires tuning a crucial hyperparameter: the kernel bandwidth. The choice of bandwidth is critical as it controls the bias-variance trade-off by over- or under-smoothing the topological features. Topological data analysis provides methods to mathematically quantify topological characteristics, such as connected components, loops, voids et cetera, even in high dimensions where visualization of density estimates is impossible. In this paper, we propose an unsupervised learning approach using a topology-based loss function for the automated and unsupervised selection of the optimal bandwidth and benchmark it against classical techniques – demonstrating its potential across different dimensions.

Climate change has driven the market to seek new ways of raising funds to mitigate its effects. One such innovation is the emergence of Green Bonds financial assets specifically designed to support sustainable projects. This study explores the fractal behavior of daily price changes in thirty-five Green Bond funds using the Multifractal Detrended Fluctuation Analysis (MFDFA) method. Our results indicate that price changes exhibit persistent behavior and high multifractality, characterized by large fluctuations. Only one of the thirty-five time series analyzed showed an outlier result, suggesting that the funds display very similar behavior. By shuffling the series, we were able to reduce multifractality significantly. These findings suggest that Green Bond funds exhibit multifractal behavior typical of other financial assets.

We seek to conduct statistical inference for a large collection of primary parameters, each with its own nuisance parameters. Our approach is partially Bayesian, in that we treat the primary parameters as fixed while we model the nuisance parameters as random and drawn from an unknown distribution which we endow with a nonparametric prior. We compute partially Bayes p-values by conditioning on nuisance parameter statistics, that is, statistics that are ancillary for the primary parameters and informative about the nuisance parameters. The proposed p-values have a Bayesian interpretation as tail areas computed with respect to the posterior distribution of the nuisance parameters. Similarly to the conditional predictive p-values of Bayarri and Berger, the partially Bayes p-values avoid double use of the data (unlike posterior predictive p-values). A key ingredient of our approach is that we model nuisance parameters hierarchically across problems; the sharing of information across problems leads to improved calibration. We illustrate the proposed partially Bayes p-values in two applications: the normal means problem with unknown variances and a location-scale model with unknown distribution shape. We model the scales via Dirichlet processes in both examples and the distribution shape via Pólya trees in the second. Our proposed partially Bayes p-values increase power and calibration compared to purely frequentist alternatives.

Accurately quantifying uncertainty of individual treatment effects (ITEs) across multiple decision points is crucial for personalized decision-making in fields such as healthcare, finance, education, and online marketplaces. Previous work has focused on predicting non-causal longitudinal estimands or constructing prediction bands for ITEs using cross-sectional data based on exchangeability assumptions. We propose a novel method for constructing prediction intervals using conformal inference techniques for time-varying ITEs with weaker assumptions than prior literature. We guarantee a lower bound for coverage, which is dependent on the degree of non-exchangeability in the data. Although our method is broadly applicable across decision-making contexts, we support our theoretical claims with simulations emulating micro-randomized trials (MRTs) – a sequential experimental design for mobile health (mHealth) studies. We demonstrate the practical utility of our method by applying it to a real-world MRT - the Intern Health Study (IHS).

With the NBA's adoption of in-game limb tracking in 2023, Sony's Hawk-Eye system now captures high-resolution, 3D poses of players and the ball 60 times per second. Linking these data to key events opens a new era in NBA analytics. Here, we leverage a large dataset of 21,964 shot attempts from 72 NBA players to introduce a novel ensemble of metrics for evaluating free-throw shooting. Inspired by baseball analytics, we introduce command, which quantifies the quality of a free throw by measuring a shooter's accuracy and precision near the basket's bullseye. This metric recognizes that some makes (or misses) are better than others and captures a player's ability to execute quality attempts consistently. We demonstrate that command captures underlying skill more effectively than traditional make-or-miss statistics; early-season command predicts late-season success more reliably than traditional shooting percentage. To identify what drives command, we define launch-based metrics assessing consistency in release velocity, angle, and 3D position. Players with greater touch, i.e., more consistent launch dynamics, exhibit stronger command as they can reliably control their shot trajectory. Finally, we develop a physics model to identify the range of launch conditions that result in a make and to determine which launch conditions are most robust to small perturbations. This framework reveals ”safe” launch regions and explains why certain players excel at free throws, providing actionable insights for player development.

For differences between means of continuous data from independent groups, the customary scale-free measure of effect is the standardized mean difference (SMD). To justify use of SMD, one should be reasonably confident that the group-level variances are equal. Empirical evidence often contradicts this assumption. Thus, we have investigated an alternate approach, based on stochastic ordering of the treatment and control distributions, that takes into account means and variances. For applying stochastic ordering, our development yields a key quantity, $\mathsf{A}$, the outcome value at which the direction of the ordering of the treatment and control distributions changes. Using an extensive simulation, we studied relative bias of point estimators of $\mathsf{A}$ and coverage and relative width of bootstrap confidence intervals.

Non-invasive measurements of the human brain using magnetic resonance imaging (MRI) have significantly improved our understanding the brain's network organization by enabling measurement of anatomical connections between brain regions (structural connectivity) and their coactivation (functional connectivity). Heritability analyses have established that genetics account for considerable intersubject variability in structural and functional connectivity. However, characterizing how genetics shape the relationship between structural and functional connectomes remains challenging, since this association is obscured by unique environmental exposures in observed data. To address this, we develop a regression analysis framework that enables characterization of the relationship between latent genetic contributions to structural and functional connectivity. Implementing the proposed framework requires estimating genetic covariance matrices in multivariate random effects models, which is computationally intractable for high-dimensional connectome data using existing methods. We introduce a constrained method-of-moments estimator that is several orders of magnitude faster than existing methods without sacrificing estimation accuracy. For the genetic regression analysis, we develop regularized estimation approaches, including ridge, lasso, and tensor regression. Applying our method to Human Connectome Project data, we find that functional connectivity is moderately predictable from structure at the genetic level (max R^2 = 0.34), though it is not directly predictable in the observed data (max R^2 = 0.03). This stark contrast suggests that unique environmental factors mask strong genetically-encoded structure-function relationships.

Stationary points or derivative zero crossings of a regression function correspond to points where a trend reverses, making their estimation scientifically important. Existing approaches to uncertainty quantification for stationary points cannot deliver valid joint inference when multiple extrema are present, an essential capability in applications where the relative locations of peaks and troughs carry scientific significance. We develop a principled framework for functions with multiple regions of monotonicity by constraining the number of stationary points. We represent each function in the diffeomorphic formulation as the composition of a simple template and a smooth bijective transformation, and show that this parameterization enables coherent joint inference on the extrema. This construction guarantees a prespecified number of stationary points and provides a direct, interpretable parameterization of their locations. We derive non-asymptotic confidence bounds and establish approximate normality for the maximum likelihood estimators, with parallel results in the Bayesian setting. Simulations and an application to brain signal estimation demonstrate the method's accuracy and interpretability.