Papers

The Bayesian Ideal Observer (IO) establishes the theoretical upper bound on task performance for binary detection tasks. However, analytical computation of the IO test statistic is generally intractable. Numerical approaches based on Markov-chain Monte Carlo (MCMC) methods, including their recent deep generative model-based extensions, typically require extensive posterior sampling for each test image. Supervised learning has also been investigated to approximate the IO performance. However, such methods are typically trained for a specific detection task and signal and may require retraining when the task or signal changes. The score function, defined as the gradient of the log probability density, encodes the local geometry of the data distribution and is a fundamental quantity in modern score-based generative modeling. This work reformulates the IO test statistic in terms of the score function and introduces a score-based ideal observer (SIO). The proposed SIO uses a denoising convolutional neural network trained exclusively on signal-absent images to estimate the signal-absent score function. Once trained, the resulting score model can be used to approximate the IO test statistic for detection tasks involving arbitrary additive signals, without per-image posterior sampling or signal-specific retraining. Numerical studies consider a signal-known-exactly (SKE) detection task with a stochastic lumpy-background model. The results demonstrate that the proposed SIO can closely approximate the IO performance.

Fast, accurate algorithms for the Faddeeva function w(z), and hence the Voigt profile K(x,a), have existed for four decades, and analytic first derivatives are available in some implementations. What the established libraries reviewed here have not provided is the full likelihood calculus of the normalized Voigt distribution: applications still commonly resort to pseudo-Voigt approximations, finite-difference derivatives, or numerical convolution for parameter inference. Because w'(z) = -2z w(z) + 2i/sqrt(pi), every derivative of the Voigt log-likelihood is an algebraic function of K and the dispersion part L(x,a) = Im w(z), from the single complex evaluation that delivers the profile. This yields the score and Hessian in closed form, and the expected Fisher information by one-dimensional quadrature of an analytic integrand; for fixed interior widths sigma, gamma > 0, the MLE of the center and both widths is consistent and asymptotically normal at rate sqrt(n), despite the distribution having no finite mean or variance, so conventional likelihood-based standard errors apply. The conditional mean of the Gaussian component given an observation is (y - mu) - gamma L/K: a redescending function that attributes moderate deviations to the Gaussian (Doppler/resolution) component and extreme ones to the Lorentzian tail. The package voigtinference (Python, NumPy/SciPy, with a cross-validated Julia companion) supplies the toolkit: score, full parameter Hessian, expected information, Newton-based unbinned maximum likelihood with boundary diagnostics, conditional component moments, and evaluation validated against high-precision references at extreme width ratios. It applies directly to unbinned non-relativistic, constant-width Breit-Wigner x Gaussian resonance fits and supplies analytic Jacobians for line-shape refinement. Companion paper: arXiv:2605.01665.

Data increasingly arrive as collections of curves - a voice recording, a growth trajectory, a day of sensor readings - where each observation is a whole function rather than a single number. The usual question asked of such data is how the average curve differs from one group to the next. But the average is only half the picture: two populations of curves can share almost the same mean and still differ profoundly in how they fluctuate around it, and it is often this variability - the pattern of covariation within a curve - that carries the scientific signal. Comparing populations at this level means comparing their covariance operators, and statistics on covariances is impaired by their non-linearity. The fdWasserstein package equips R users with functions to make such comparisons. It is centered on the geometry of optimal transport, under which covariance operators can be meaningfully averaged, contrasted, and interpolated. It provides the Procrustes-Wasserstein distance between covariance operators, their Frechet mean (barycenter), an ANOVA-type permutation test for the equality of several covariances, principal component analysis of covariance variation, and an entropy-regularized soft clustering of curves by their covariance structure. We outline the underlying ideas, discuss the implementation and demonstrate the complete workflow on the phoneme data shipped with the package.

Semi-supervised classifiers are commonly trained from samples in which all features are observed but some class labels are missing. When label missingness is independent of the observed data, unavailable class memberships reduce Fisher information relative to a completely classified sample. We study a different regime in which the probability of label missingness depends on posterior classification uncertainty, so that the observed missing-label indicators can themselves carry information about the Bayes decision boundary. Building on the conditionally weighted information decomposition of Ahfock and McLachlan, we develop this phenomenon for a two-component exponential mixture. Although the exponential model is non-Gaussian, asymmetric, and supported on the positive half-line, its log-posterior odds remain linear in the feature. We derive Bayes' rule and its exact error rate, formulate entropy-logistic and squared-discriminant missingness mechanisms, and obtain the full partially classified likelihood. We then derive a decomposition of the Fisher information into the complete-data information, the conditionally weighted loss due to missing labels, and the information contributed by the missing labels. Numerical quadrature identifies regions in which the full likelihood classifier has asymptotic relative efficiency above or below one. Monte Carlo experiments with finite training samples broadly support the population calculations, with the largest departures from the asymptotic predictions occurring near the transition at which the relative efficiency crosses one.

Theoretical guarantees for kernel regression are typically formulated in terms of smoothness, but practical accuracy depends critically on how the design resolution compares with the target and kernel lengthscales. We develop a finite-sample theory for Matérn regression on periodic domains with quasi-uniform designs, covering noiseless interpolation and noisy kernel ridge regression. Our minimax result shows that accurate recovery requires a design dense enough to resolve the target lengthscale and sufficient information at that scale to overcome noise. We prove that Matérn interpolation additionally requires the design to resolve the kernel lengthscale: if the kernel is too short relative to the point spacing, an additional error remains even when the target is well resolved. For noisy Matérn regression, we derive a three-term fixed-ridge risk characterization consisting of target-scale bias, kernel-scale bias, and variance, and show that optimizing over the ridge parameter yields four distinct contributions. When the target is no more than twice as smooth as the kernel, choosing a kernel lengthscale longer than the target does not worsen the oracle risk, whereas choosing one too short can. Thus these resolution conditions determine when accurate recovery becomes possible, while smoothness determines how rapidly the error decreases thereafter.

Two-factor ANOVA is widely used in experimental studies but relies on additivity, normality, independence, and homoscedasticity. These assumptions are often violated for nonnegative, positively skewed observations. Although Box–Cox-type transformations are commonly used, they may reduce interpretability and require a subjective choice of transformation. We propose an alternative framework in which nonnegative observations affected by two factors are modeled by gamma distributions with unknown shape and scale parameters that may depend on factor levels. We develop likelihood ratio tests (LRTs) for main and interaction effects. The asymptotic LRT (ALRT) uses the asymptotic chi-square distribution, which may be inaccurate for small to moderate samples. We therefore propose a parametric bootstrap LRT (PBLRT) that determines critical values by simulation. Extensive simulations show that the PBLRT maintains the nominal significance level well. Real-data examples demonstrate its applicability and show that its inferences can differ from those of traditional ANOVA.

This note revisits the convergence of the Parikh–Boyd two-way block-splitting algorithm for large-scale distributed learning through the He–Yuan prediction–correction framework. Simultaneous row–column partitioning is also relevant to hybrid federated learning, where data may be heterogeneous in both samples and features. We lift the reduced iteration to an equal-dimensional product space and reconstruct the primal and dual coordinates omitted by its implementation. The induced orthogonal-complement structure establishes exact iteration-by-iteration equivalence with the published updates. A mixed variational-inequality representation then yields a fundamental descent inequality, global convergence, an ergodic complexity bound, and current-iterate residual estimates under standard convexity, solvability, exact-subproblem, and invariant-initialization assumptions. The analysis also shows that the reduced state recursion is a metric proximal point iteration. No strong convexity, differentiability, or full-rank condition is imposed. The derivation clarifies which algebraic initialization conditions allow the reduced implementation to inherit the full-space convergence and complexity guarantees without modifying its local updates.

A common question about networked time series is whether outcomes changed because shocks transmit more strongly or because the pattern of connections changed. Standard practice inserts a recorded network into an outcome regression and reads movements of the fitted coefficient as changes in transmission strength. When the network is latent, time varying, and measured with error, this reading fails: changes in strength and changes in composition can produce the same outcome distribution at a single date, and the population coefficient moves under composition changes alone. The question becomes answerable when outcomes are analyzed jointly with repeated noisy measurements of the network, such as paired reports of bilateral trade flows. For the joint model we establish necessary and sufficient conditions for local identification, estimators of the strength and composition paths, a simultaneous confidence band for the strength path, confidence sets that remain exact under weak identification, breakdown bounds under common reporting bias, and an exactly sized test that detects changes on the observed path and attributes them to strength or to composition. Simulations assess each procedure at its stated boundary. On a mirror-reported trade panel of eighteen economies over 1995 to 2020, the diagnostics flag exactly the crisis years and the composition coordinate attached to European Union membership declines by roughly two thirds. The estimand is predictive dependence, not a causal effect.

The multivariate expectile-based distribution is a `spiked' perturbation of a multivariate Gaussian distribution, proposed by Arbel et al. in 2023. The rejection sampler proposed for this distribution in the initial work is valid, but its acceptance probability can degenerate badly both in high dimension and in the strongly asymmetric limit. We give a simple alternative. After affine whitening and a polar decomposition, the sampling problem reduces exactly to a univariate distribution, and an additional hyperbolic change of variables exhibits this distribution as log-concave, so that a universal and uniformly efficient construction of Devroye applies. The resulting exact sampler therefore enjoys an acceptance probability of at least $1 - \exp \left( -1 \right) = 0.632120 \ldots$, uniformly over the dimension and all admissible asymmetry parameters.

Multidimensional array data, or tensors, arise naturally in neuroimaging and other high-dimensional applications. We propose a parsimonious Bayesian tensor regression model for studies in which a brain image is the response and predictors are vector-valued covariates. The method extends Bayesian envelope dimension reduction to tensor responses, identifying material subspaces that contain regression information while removing variation that is immaterial to the predictors. This formulation leads naturally to a Tucker tensor decomposition and allows spatial dependence and multiple sources of uncertainty to be modeled jointly. We develop a computationally feasible Markov chain Monte Carlo algorithm based on Gibbs sampling and establish posterior consistency for the proposed model. Simulation studies demonstrate substantial gains in estimation accuracy and uncertainty quantification when meaningful dimension reduction is present. We apply the method to Human Connectome Project neuroimaging data to investigate associations between alcohol use and brain activity. The results illustrate the value of Bayesian tensor envelope regression for inference with high-dimensional, spatially dependent imaging responses.

The $K$-function is a fundamental summary statistic for assessing clustering or regularity of point processes in two or three dimensional Euclidean space. In practice, however, many planar point patterns arise from projecting locations of objects on a surface in three dimensional space to two dimensional space. For example, when events or objects occur in a landscape their elevation is often ignored. This can lead to erroneous conclusions regarding properties of the point process generating the point pattern. There is not a unique way to extend the classical $K$-function to point patterns on a complex surface. In this paper we propose, explore, and discuss several approaches in terms of their theoretical and computational properties. The best performing approach, coined the surface area $K$-function, can be viewed is an analogue of the classical $K$-function replacing counts of points in Euclidean balls with counts of points in surface geodesic balls. However, an important distinction is that the argument of our surface area $K$-function is area instead of radius of geodesic balls. The performances of the various surface $K$-functions are compared in applications to simulated and real data.

Complex-valued arrays arise in signal processing, where scientific interpretation depends on retaining amplitude and phase information. Existing covariance estimation methods either ignore the multiway organization of such data or rely on real-domain embeddings that do not directly exploit their complex structure. We develop sparse separable factor analysis (SSFA), a latent factor model for complex-valued arrays with a separable covariance structure across modes. Each mode-specific covariance matrix is modeled through a low-rank Hermitian factor structure and a diagonal residual covariance matrix. To obtain interpretable estimates, we impose elementwise lasso penalties on the complex loading matrices and estimate the SSFA parameters using a mode-wise parameter-expanded expectation-maximization procedure. The resulting loading updates admit closed-form complex soft-thresholding solutions, which shrink the modulus of each loading while preserving its phase. A separate balancing step resolves the scale nonidentifiability of the separable covariance structure. Simulation studies show that SSFA improves covariance estimation relative to vectorization-based methods, including complex principal component analysis. We apply SSFA to local field potential recordings from mice, where we compare separability structures induced by different groupings of brain region, frequency, and time and perform model-based imputation of recordings missing because of electrode misplacement.

We develop spectral algorithms for selecting state-space partitions that define averaging kernels for finite, ergodic and reversible Markov chains. For a partition $\mathcal O$, the Gibbs kernel $G_{\mathcal O}$ resamples within the current block from the stationary conditional distribution; when this update is tractable, composing or mixing it with a baseline kernel $P$ can accelerate convergence. We select $\mathcal O$ by rounding the bottom nonconstant eigenfunctions of $P^2$, or the algebraically smallest eigenfunctions of $P$ for additive mixtures, using weighted $k$-means. For $F(\mathcal O)=\|G_{\mathcal O}P-\Pi\|_{F,\pi}^2$, we derive exact trace and normalized-cut representations and show that $F$ equals the Pearson $\chi^2$-mutual information between the initial block label and the state after one transition, giving this matrix objective a natural probabilistic interpretation. In the two-block case, a threshold sweep exactly solves the associated one-dimensional weighted two-means rounding problem. For general $k \geq 2$, weighted $k$-means rounds the bottom $(k-1)$-dimensional embedding, after which candidates are rescored by $F$; the rounding distortion is a distance between subspaces that yields spectral approximation bounds. We extend the framework to additive mixtures, finite-horizon objectives, and discounted infinite-horizon objectives. In contrast to classical normalized spectral clustering, which uses top nonconstant modes to find low-flow persistent clusters, our method uses bottom modes to favor large normalized cross-block flow and rapid loss of block-label information. Experiments on a controlled-spectrum graph, a mean-field Ising model, and Bayesian variable selection show notable per-iteration improvements in convergence and statistical estimation.

Continuous positive data characterized by high skewness and heavy tails frequently arise in applied statistics. In other applications, these characteristics are accompanied by a point mass at zero, resulting in a non-negative response with a mixed discrete-continuous distribution. Standard regression models often fail to capture these complex features adequately, requiring more flexible approaches. In this paper, we introduce the BCSreg package for R, which provides a comprehensive and unified computational framework for fitting Box-Cox symmetric and log-symmetric regression models for positive continuous data and their zero-adjusted extensions for mixed non-negative data. These broad classes of models accommodate varying degrees of skewness and tail-heaviness while allowing the parameters to be interpreted directly on the original scale of the data. Through a user-friendly multi-part formula interface, the BCSreg package allows practitioners to simultaneously specify regression structures for the scale parameter (which is proportional to the quantiles of the response), the relative dispersion, and, when appropriate, the probability of zero occurrences. Furthermore, the package provides a complete suite of diagnostic tools specifically tailored to these classes of models, including randomized quantile residuals, simulated envelopes, and influence diagnostics. The package's features and capabilities are illustrated through applications to real data.

Binary segmentation is widely used for multiple change-point detection because it is fast, simple to describe, and simple to implement. Its validity rests on the requirement that, at each recursive stage, the procedure identifies one of the true change-points when several are present in the current interval. This holds for detecting changes in mean using the CUSUM statistic, but fails in some other settings, in particular in slope change detection for continuous piecewise-linear signals. We propose Look-Ahead Binary Segmentation (LABS), a modification in which the change-points returned by the two child recursions define a narrower interval on which the parent estimate is re-evaluated. LABS inherits the computational speed of standard binary segmentation, but achieves the near-optimal consistency rate of $O\{(n\log n)^{1/2}\}$ in the slope-change signal setting when the LABS model is chosen via either thresholding or a Schwarz-like information criterion. Simulations show that LABS is fast and achieves state-of-the-art performance.

The randomly pivoted Cholesky algorithm is one of the leading methods for computing a low-rank approximation to a large positive-semidefinite matrix. However, while it consistently achieves accuracy comparable to or better than competing methods of its type in experiments, its theoretical analysis lags somewhat behind other methods. This paper closes this gap, proving that randomly pivoted Cholesky produces an approximation with expected error within a $1+\varepsilon$ factor of the optimal rank-$r$ approximation in $\mathcal{O}(r/\varepsilon + r\sqrt{\log r})$ steps. This result nearly matches the optimal complexity $\Theta(r/\varepsilon)$ for any low-rank approximation method based on a partial Cholesky decomposition (also known as a column Nyström approximation). The paper also presents bounds on the randomly pivoted Cholesky trace and spectral-norm errors that hold with high probability. The mathematical argument is largely due to GPT 5.6-Sol (Pro), with some refinements by the author.

The standard approach to stochastic simulation from the von Mises-Fisher distribution is a rejection sampler proposed by Ulrich and Wood. This note provides a theoretical justification of its efficiency, lower-bounding its acceptance probability uniformly over the dimension and concentration parameter. A novel interpretation of the proposal is given which provides some intuition for this efficiency.

Many scientific problems are relational: the quantity of interest is a connection between two objects, while information about similarity is available for the objects themselves. We develop a covariance framework for unordered relationships that transfers object-level geometry to the relations they form while preserving endpoint identity and invariance to ordering. Building on symmetric pairwise-kernel representations, we develop theory for the loop-free domains used in undirected networks. We establish spectral interlacing and trace-loss results after self-pairs are removed, connect the relational spectrum to regularization and risk, derive an exact inferential error for spectral truncation, and quantify how perturbations of the underlying geometry propagate to pair covariance and estimation. Simulations show when structured borrowing improves estimation and how geometric misspecification can erode that benefit. We apply the framework to autism neuroimaging using resting-state functional-connectivity data from the Autism Brain Imaging Data Exchange (ABIDE). With a 116-region parcellation and 6,670 unique connections, the application shows that a large connectome can have a much smaller effective covariance dimension. It also demonstrates that high explained covariance alone is insufficient for choosing a low-rank representation when inferential accuracy is the goal. The framework provides a principled foundation for covariance and regularization when the statistical units are unordered relationships.

In many applications, from government to ecology, integrating data from diverse and noisy sources is critical for downstream inference. However, a unique identifier to link records cleanly from the same entity may not exist. Entity resolution (also referred to as de-duplication or record linkage) merges such databases to identify duplicates, allowing for more complete data for inference. A multitude of methods from statistics and computer science in recent years assume a single, immutable true value for each variable used in linkage. We argue this restrictive assumption is often violated in practice, leading researchers to discard useful data and biasing downstream inference. For example, a respondent's education status may truly change between two different surveys, yet existing methods would treat at least one observation as a distorted version of the truth, even though both are correct. In this paper, we propose a novel entity resolution model that comprehensively accommodates "multiple truths" by introducing a mixture of exponential family distributions that further handles multiple data types. We provide options to fit the model with Markov chain Monte Carlo routines and variational inference for massive datasets. We demonstrate the method's value via simulation and linking a longitudinal survey of Italian household wealth.

The Bingham distribution is a family of antipodally symmetric distributions on the unit sphere, characterised by an exponential-of-quadratic change of measure with respect to the uniform distribution. Kent, Ganeiber and Mardia proposed a rejection sampler for generating samples from Bingham distributions using proposals from an angular central Gaussian (ACG) distribution. Their empirical results suggest that the least efficient regime is the high-concentration limit, where the acceptance probability is of order $d^{-1/2}$ in dimension $d$, implying a polynomial complexity guarantee. In this note, we verify this dimension-dependent prediction, establishing the uniform guarantee $\inf\{\alpha_D:D=D^\top\in\mathbb{R}^{d\times d}\}\ge c_\star/\sqrt{d}$, where $c_\star=0.759\ldots$. A one-dimensional high-concentration limit demonstrates that the $d^{-1/2}$ rate is unimprovable and that even the constant $c_\star$ cannot be improved beyond $0.857\ldots$. The proof relies on a novel interpretation of the acceptance probability and a comparison principle for weighted sums of chi-squared random variables, which may be of independent interest.

In Markov chain Monte Carlo sampling, light-tailed target distributions present something of a poisoned chalice: their light tails offer good confinement, and tend to imply good mixing properties for natural continuous-time dynamics, but the steepness of their tail decay means that they often fall outside of the scope of modern quantitative convergence theory. For usual gradient-based samplers, this reflects a genuine instability issue, whereby Metropolis acceptance rates can degrade badly. In this work, we study the gradient-free random-walk Metropolis sampler, and show that for a wide range of light-tailed targets, the acceptance probability remains stable for reasonable choices of proposal variance, from which effective and favourable mixing time estimates can be deduced. The analysis relies on a simple relationship between the first and second derivatives of the log-density of the target distribution.

Autism spectrum disorder (ASD) is associated with heterogeneous alterations across distributed brain systems, creating challenges for whole-connectome inference. The difficulty arises not only from the large number of connections, but also from dependence among effects indexed by anatomically and functionally related region pairs. We introduce a Bayesian Edge-Space regression framework that treats each participant's connectome as a network-valued response and models the adjusted ASD effect over unordered brain-region pairs. The main methodological contribution is a positive-semidefinite covariance construction defined directly on connections. Anatomical and diagnosis-blind functional similarities are lifted from regions to edge space through a symmetrized endpoint-matching operation that preserves endpoint identity and is invariant to endpoint ordering. An additive Bayesian hierarchy estimates anatomical, functional, and interaction contributions together with multisite adjustments and connection-specific effects. Theoretical results establish covariance validity and continuous nesting of the structured components. Low-rank kernel representations and an exact sufficient-statistic reduction enable whole-connectome computation without preliminary edgewise estimation. Simulations show improved recovery of the effect surface, particularly under weak signals. In the Autism Brain Imaging Data Exchange, the framework identifies widespread reductions together with localized increases in ASD-associated connectivity. This pattern supports heterogeneous reorganization across distributed neural systems rather than uniform hyper- or hypoconnectivity. Under the fitted parameterization, the functional component has the largest structural scale, indicating organization beyond anatomical proximity alone.

Skchange is an open-source Python library for detecting structural changes in time series. It implements modern change detection algorithms within a unified and extensible framework. The algorithms are modular and composable, and they include changepoint search methods based on both cost minimisation and statistical tests. Key features include the detection of anomalous segments in addition to changepoints; theoretically well-founded fast and approximate search methods; theoretically well-founded algorithms for high-dimensional data, covering settings where either few or many features change simultaneously; utilities for automatic and data-driven penalty calibration, which balances false alarms against missed detections; and a large collection of built-in costs and statistical tests. The design follows established scikit-learn conventions to streamline both user and contributor experience, and Numba is used extensively to achieve high computational performance. Source code and documentation are available at https://github.com/NorskRegnesentral/skchange.

Coccidioidomycosis is a fungal lung infection endemic to parts of the Pacific Northwest and southwestern United States, Mexico, Central America, and South America. A 2017 study by Thompson et al. reported that many analyzed Coccidioides isolates had elevated minimum inhibitory concentration values for fluconazole but comparatively low minimum inhibitory concentration values for other triazole drugs. We constructed 172 synthetic joint minimum inhibitory concentration profiles from the published drug-specific marginal frequency tables. Values sampled from each marginal distribution were paired across drugs using a Gaussian copula, and potential cross-drug patterns were examined using k-means clustering, hierarchical clustering, and Gaussian mixture models. On the selected reconstruction, a forced four-cluster partition consistently identified an upper caspofungin minimum-inhibitory-concentration tail. Across the evaluated dependence scenarios, the gap statistic favored a single cluster in 96% to 100% of nested reconstructions, providing no evidence that the reconstructed multivariate data supported a broader multicluster structure. A separate simulation study examined recovery of an upper-severity-score category defined from prespecified quantiles of a composite minimum-inhibitory-concentration score. The study included 4,000 replications across 20 prevalence and measurement-noise conditions. The multiclass adjusted Rand index ranged from moderate to high depending on the method and prevalence, whereas binary partition agreement for the upper-severity-score category approached zero at 1% prevalence even though recall remained near 1.0. Clustering conclusions therefore depended on the assumed cross-drug dependence structure and on the realized number of observations in the target category.

We conduct a retrospective analysis of COVID-19 transmission dynamics across U.S. states using a modified population-based Susceptible-Exposed-Infectious-Asymptomatic-Recovered-Deceased (SEIARD) compartmental model. The proposed framework introduces time-varying transmission, reporting, and mortality rates to capture temporal variations in public behavior and policy interventions during the pandemic. In particular, the transmission rate is modeled as a function of population mobility (derived from Google Mobility Reports), with a residual time-decay term capturing the net effect of unobserved factors such as behavioral adaptation and control measures, while reporting is linked to nationwide testing strategies. We employ a Bayesian approach to integrate multiple data sources and quantify uncertainties in model parameters. The model explicitly distinguishes between symptomatic and asymptomatic infectious individuals and links the latent epidemic states to observable quantities, including reported cases and deaths, through a dynamic reporting function. This retrospective modeling framework provides insights into state-level epidemic trajectories and supports data-driven decision-making for optimal allocation of healthcare resources and evaluation of public health interventions during future pandemics. We further apply a clustering analysis to the posterior parameter estimates to identify groups of U.S. states exhibiting similar epidemiological characteristics, revealing substantial regional heterogeneity in transmission intensity, reproduction dynamics, and mortality burden.

Goal-oriented optimal experimental design (GO-OED) selects experiments according to the expected information gain (EIG) about a quantity of interest (QoI) rather than the full parameter vector. This work develops a variational GO-OED formulation for mixed discrete-continuous QoI laws arising in probabilistic mechanics when thresholding or event-based transformations map a positive-probability set of uncertain inputs to a common value while other inputs produce continuously varying responses. The motivating application is ship roll safety assessment in random waves, where the QoI is the temporal exceedance probability above a prescribed roll-angle threshold. This quantity is zero when no exceedance occurs and varies continuously over positive values otherwise. A purely continuous variational approximation does not dominate a posterior QoI law containing an atom, yielding an infinite Kullback-Leibler divergence and a trivial Barber-Agakov lower bound of $-\infty$. Scoring atom samples using continuous density values instead changes the objective and does not produce a valid lower-bound estimator. We introduce a mixed variational approximation that models the conditional atom probability and continuous component separately, with a normalizing flow used for the latter. An analytical example recovers the correct EIG landscape, while the ship roll application provides stable EIG lower-bound estimates and identifies informative wave conditions for temporal-exceedance-probability inference.

To address the computational and storage challenges often encountered in large-scale survival data analysis, we propose an efficient Poisson subsampling method for the partially linear additive Cox model. This model provides a flexible yet interpretable framework by incorporating linear covariate effects, additive nonparametric components for nonlinear covariates, and a nonparametric baseline hazard function. The proposed method adopts B-spline basis functions to approximate the nonparametric components and employs the decorrelated score technique to construct a Poisson subsampling-based estimation equation, based on which we establish the asymptotic normality of the resulting estimator and derive the optimal subsampling probabilities according to the L-optimality criterion. Furthermore, we design a two-step adaptive algorithm for practical implementation. The proposed approach enables computationally efficient statistical inference for large-scale survival analysis without processing the full dataset. We validate the performance of the proposed method through extensive simulation studies and a real-world application to a lymphoma cancer dataset, demonstrating its efficiency and accuracy in large-scale settings.

Bayesian inference via Markov Chain Monte Carlo (MCMC) provides effective parameter estimation, but its real-time application in complex physical systems is hindered by heavy computational bottlenecks and extreme sensitivity to statistical noise. We address this by proposing a neural-network-based probabilistic surrogate framework for rapid and robust MCMC inference. Using fusion plasma Thomson scattering diagnostics as a challenging, noise-dominated testbed, our approach employs a dual-head architecture to simultaneously estimate the expected physical emission spectrum and the channel-wise intrinsic measurement noise variance. By optimizing a Gaussian Negative Log-Likelihood (GNLL) objective, the learned aleatoric uncertainty dynamically buffers the sampler against pathological shot noise. Evaluations demonstrate that this surrogate framework achieves > 1500x acceleration over exact physical forward models, while simultaneously reducing inference error (RMSE) by >20% compared to standard homoscedastic neural baselines, offering a highly promising paradigm for real-time physical analysis.

Transportation provides access to employment opportunities and essential services such as healthcare services, while urban areas have various transportation options, the situation differs in rural areas. Rural residents often have longer commute distances, limited access to public transit, and extended waiting times for public transportation if they exist, which can significantly impact their access to vital services and job opportunities. This study used data from the 2017 NHTS survey to examine the commuting patterns in rural areas by utilizing multinomial logistic regression to determine how various factors impact the choice of mode of transport in rural areas. Findings from this study revealed a higher dependency, 92.1%, on using personal vehicles when making trips in rural areas. Multinomial logistic regression results showed that socio-demographics, household, and trip characteristics affect the mode of transport used in a trip. Older adults, females, and individuals with higher education levels than high school graduates are less likely to use public transit when making trips. For household characteristics, the availability of vehicles in a household and households with higher income levels have lower probabilities of making a trip using public transit. Longer trip distances reduce the likelihood of a trip using active commuting modes such as walking and biking. These findings provide insights into understanding the transportation behaviors in rural areas and provide knowledge to be used in the planning and developing of transportation projects to promote equitable and accessible transportation in rural areas.

Composite likelihood (CL) methods provide a computationally efficient alternative to full likelihood inference for complex multivariate models by replacing the joint likelihood with a product of lower-dimensional marginal or conditional components. Like the MLE, however, the maximum CL estimator (MCLE) is highly sensitive to data contamination. On the other hand, robust divergence-based procedures such as the minimum density power divergence (DPD) estimator require the full joint density and so scale poorly to complex multivariate models. We introduce the composite DPD (CDPD), a genuine statistical divergence built entirely from the low-dimensional component densities defining a CL, combining the computational scalability of CL with the robustness of the DPD. The resulting minimum CDPD estimator (MCDPDE) robustifies the MCLE without requiring integration over the full multivariate sample space. We establish consistency, asymptotic normality, and the influence function of the MCDPDE under regularity conditions on the component models alone, without requiring correct specification of the full joint distribution. We show that it is qualitatively robust for every positive value of its tuning parameter, unlike the MCLE recovered as the limit. Because its components can be chosen at the pairwise or cell level, the framework guards simultaneously against casewise and cellwise contamination. Operating directly on component densities rather than elliptical distance structures, it extends robust inference beyond the elliptical models to which most existing cellwise-robust procedures are confined. We develop computational algorithms implemented in the accompanying R package mvdpd. Simulation studies and real-data applications show that the MCDPDE achieves substantial robustness gains over the MCLE while retaining competitive efficiency under the assumed model.