Papers

Building on a reparameterization for multivariate linear regression that yields a jointly convex penalized likelihood in the reparameterized regression coefficient matrix and the precision matrix, we show that the resulting scaled Gaussian loss is standard self-concordant. This places the joint estimation problem within composite self-concordant optimization and leads to two algorithms: a proximal gradient method and a damped proximal Newton method. In simulations, we evaluate algorithmic robustness, iterations to convergence, and elapsed time. In a protein expression application, compared with the classical-parameterization formulation, the proposed convex formulation attains similar mean squared prediction error and can be substantially faster when the fitted precision matrix is dense.

I develop a reproducible framework for evaluating individual batting and bowling performances in Twenty20 (T20) cricket on one interpretable scale of runs above expectation, built from two ball-level primitives. The first, Runs Above Expected (RAE), is the residual between the runs scored on a delivery and a contextual expectation of what an average performer would produce in the same situation. That expectation is a multiplicative log-linear model of the cohort scoring rate, whose per-cell estimator is shown to be a conditional Poisson maximum likelihood multiplier. It is fitted by iterated backfitting and stabilised by empirical Bayes shrinkage, so that thinly sampled contexts are pooled towards the population. A single opposition symmetry places run-scoring and run-prevention on the same footing. The second primitive, Dismissal Adjusted Runs (DAR), prices a dismissal in that currency as the runs it forgoes, read off a batting side value function solved by dynamic programming through a Bellman expectation recursion, under observed play. Because dismissal is the expected end of every innings, the realised wicket cost (realDAR) is centered against its expectation (xDAR) under a league dismissal hazard rate. The centered quantity is a run-weighted mean zero martingale residual of the dismissal process, so a player is charged only for departing from average behaviour. The two primitives sum to a symmetric Impact, which makes batting and bowling comparable in centre as well as in unit. Estimated on over 2.7 million legal deliveries of men's T20 cricket, the framework recovers known contextual structure, agrees with the conventional rates it refines while correcting their context-blindness, and yields face-valid player, innings and season leaderboards for the Indian Premier League.

Spectral imaging provides chemically specific and spatially resolved analysis of microplastics, but its routine application is hindered by large data volumes, acquisition artefacts, spectral variability, and misidentification of polymers due to alike spectra. This study proposes a multi-level preprocessing and modelling framework for FT-IR spectral imaging of microplastics that integrates image-level, tile-level, and spectral-level corrections with scalable identification strategies. Image-level variation associated with changing acquisition conditions was done with latent variable selection, while a background-based tile correction reduced illumination-related artefacts. Spectral preprocessing combined baseline correction, smoothing, derivative calculation, normalization, and wavelength selection, and only particle spectra were retained for further analysis to improve computational efficiency. For scalable identification, clustering was applied to particle spectra and spectral library matching was performed on cluster centroids instead of individual pixels. Among twelve evaluated matching strategies, a sign-invariant derivative-based cosine similarity method achieved perfect classification accuracy for polystyrene (PS), polyethylene terephthalate (PET), polyethylene (PE), and polypropylene (PP). The clustering-based workflow also produced more spatially coherent particle maps than direct software-based matching while substantially reducing processing time. The framework was evaluated for supervised classification-based MP indentification. These results show that multi-level correction combined with cluster-centroid spectral matching improves the robustness, efficiency, and interpretability of spectral-imaging-based microplastic identification.

We study uniform sampling of binary matrices with fixed row and column sums, a recurring problem in ecological null models, Rasch-model testing, network analysis, and combinatorics. We propose the Snake algorithm, a rejection-free Markov chain Monte Carlo sampler that grows an alternating path until its first self-intersection and flips the resulting loop. The chain is reversible and irreducible on the fixed-margin state space, hence has the uniform stationary distribution. We prove that one step flips on the order of $\sqrt{n}$ entries in sparse and balanced square regimes, give upper bounds on the per-step path length, and show that the resulting work per flipped entry is rate optimal in sparse and balanced regimes and near-optimal up to a polylogarithmic factor under a one-sided half-balanced condition. A Markov-chain comparison, combined with the recently established universal spectral-gap bound for the swap chain, proves that the lazy Snake chain is rapidly mixing for every feasible pair of margins; in the permutation-matrix case, the raw chain has the sharp total-variation mixing time $\Theta(n \log n)$. We also describe a directed-graph extension and an equal-margin label-shuffling variant. Numerical experiments against Swap, Rectangle Loop, Curveball, sequential importance sampling, and a directed edge-swap algorithm show consistent gains in move size, wall-clock convergence, and sampling efficiency.

This article investigates several physics-informed and hybrid machine learning strategies that incorporate physics knowledge in experimental data-driven deep-learning models for predicting the bond quality and porosity of fused filament fabrication (FFF) parts. Three types of strategies are explored to incorporate physics constraints and multi-physics FFF simulation results into a deep neural network (DNN), thus ensuring consistency with physical laws: (1) incorporate physics constraints within the loss function of the DNN, (2) use physics model outputs as additional inputs to the DNN model, and (3) pre-train a DNN model with physics model input-output and then update it with experimental data. These strategies help to enforce a physically consistent relationship between bond quality and tensile strength, thus making porosity predictions physically meaningful. Eight different combinations of the above strategies are investigated. The results show how the combination of multiple strategies produces accurate machine learning models even with limited experimental data.

This paper addresses the limitations of Gaussian distribution assumptions in statistical sparse learning, particularly in modeling correlated and heterogeneous data. Conventional Gaussian models often lack robustness towards outliers and underlying distribution assumptions. To overcome these limitations, we propose the use of the $q$Gaussian distribution, derived from Tsallis entropy maximization, as a robust alternative. This is notably relevant in biostatistics, where the presence of correlated observations and heterogeneity, such as in genetic and longitudinal studies, is prevalent. Our contributions include modeling of correlated data through the re-derived multivariate probability density function from Tsallis entropy maximization, thereby addressing the limitations inherent in conventional Gaussian models. Furthermore, we introduce a novel framework that adapts numerical methods designed to find equilibria in flows to tackle composite optimization problems prevalent in statistical sparse learning. Applying this framework to the Hager-Zhang conjugate gradient algorithm \cite{Hager2005}, we develop a numerically stable and efficient algorithm for sparse statistical learning. The $q$Gaussian distribution, informed by the principle of maximizing Tsallis entropy, presents a viable and flexible alternative to Gaussian-based methods. This paper not only contributes to the theoretical understanding of statistical distributions and optimization techniques, but also paves the way for practical data analysis.

Quasi-Monte Carlo sampling is a numerical integration method that uses points with a space-filling property in $[0,1]^s$ designed to give better estimates than plain Monte Carlo methods do. For integrands of bounded variation in the sense of Hardy and Krause, errors of $O(n^{-1+\epsilon})$ for any $\epsilon>0$ are obtained from $n$ sample points. Randomized quasi-Monte Carlo (RQMC) points are individually uniformly distributed but collectively space-filling and then independent replications provide variance estimates. For smooth enough integrands the randomization can give a root mean squared error of $O(n^{-3/2+\epsilon})$. This article explains RQMC for a statistical readership recounting some history and presenting some current directions.

Generalised Bayesian inference uses weights of the form $\exp\{-\beta_n\ell_n(\theta)\}$ even when the loss is only estimated. Exponentiating an unbiased loss estimate changes the target, and when $\beta_n\asymp n$ an ordinary Monte Carlo loss estimate with variance of order $M^{-1}$ requires a per-proposal budget $M$ of order $n^2$ to keep the leading log-weight variance bounded. We introduce Sign-Corrected Bessel Debiasing (SCBD), a signed pseudo-marginal method based on independent block estimates of the loss, and study its ordinary-MC and independently randomised quasi-Monte Carlo (RQMC) implementations. Under an i.i.d. Gaussian block model, a Bessel factor constructed from the block sample variance exactly removes the Gaussian exponential inflation despite the variance being unknown. For general non-Gaussian finite blocks, the method targets a posterior differing from the intended posterior by a parameter-dependent multiplicative factor. Under regularity conditions, the uncorrected and corrected ordinary-MC targets have total-variation errors of orders $\beta_n^2/M_n$ and $\beta_n^3/M_n^2$. If an RQMC block estimator has variance $\mathcal O\{B^{-\alpha}(\log B)^{d-1}\}$, the corresponding errors are of orders $\delta^{\mathrm{RQ}}_{n,M_n}$ and $(\delta^{\mathrm{RQ}}_{n,M_n})^{3/2}$, where $\delta^{\mathrm{RQ}}_{n,M_n}=\beta_n^2M_n^{-\alpha}{\log(2+M_n)}^{d-1}$. The same variance rate gives a sufficient budget of order $n^{2/\alpha}$, up to logarithmic factors, for bounded leading log-weight variance when $\beta_n\asymp n$. The numerical examples show that favourable RQMC representations can inherit this budget scaling and that variance reduction and Bessel correction are complementary. Compared to existing exact corrections, Bessel debiasing is essentially "for free". It is generic, easy to code and supported by theory.

Ornstein-Uhlenbeck processes with flexible multivariate drift matrices are powerful models for capturing coupled, asymmetric, and damped-oscillatory mean reversion. However, likelihood-based inference is challenging because the drift matrix must remain Hurwitz stable, while likelihood and gradient evaluations require repeated, costly computation and differentiation of the drift matrix exponential and of solutions to the associated Lyapunov equation. We introduce the Hurwitz smooth spectral block parametrization (H-SSBP), which represents the drift as a change of basis applied to independent one- and two-dimensional stable blocks. Simple scalar constraints enforce Hurwitz stability, while the two-dimensional blocks vary smoothly between real- and complex-eigenvalue regimes, avoiding discrete model selection. The H-SSBP represents every real Hurwitz matrix diagonalizable over the complex numbers and has dense, full-measure support within the Hurwitz cone. Its block structure reduces the OU transition matrix, stationary and innovation covariances, and their reverse-mode derivatives to constant-size block or block-pair computations plus change-of-basis multiplications. Numerical experiments demonstrate dramatic speed-ups over alternatives, particularly for matrix-exponential adjoints and Lyapunov-equation kernels. Real-data analyses of asynchronous financial data and multivariate phylogenetic traits illustrate the proposed framework in practice.

Estimating high dimensional Generalized Structural Equation Models presents severe computational challenges. Traditional simultaneous estimators frequently suffer from numerical instability and prohibitive computational costs. Moreover, there are no tractable algorithms for families such as Poisson, negative binomial, and gamma. To overcome these limitations, this article introduces a Two Stage Quasi-Likelihood Expectation-Maximization framework. The proposed method isolates the structural model from the measurement model. First, it approximates the conditional distribution of the latent variables given the observed indicators. Second, it employs marginal quasi-likelihood estimating equations to evaluate the structural parameters, deriving the necessary conditional moments either exactly or through Monte Carlo integration. This approach completely avoids the need to evaluate the full joint likelihood. Extensive simulations demonstrate that our method drastically reduces computational runtime, providing a numerically stable framework that minimizes the mean squared error and structural bias to yield a scalable and flexible solution for analyzing complex latent variable models.

We present a computational framework for robust topology optimization that integrates polygonal finite-element discretizations, spatially correlated non-Gaussian material modeling, and non-intrusive polynomial-chaos surrogates. Spatial uncertainty in Young's modulus is represented as a homogeneous non-Gaussian random field obtained via a memoryless transformation of a truncated Karhunen-Loève expansion, ensuring physical admissibility through positivity of stiffness while preserving the prescribed autocovariance. Polygonal finite elements provide a stable discretization for density-based optimization on unstructured meshes and mitigate checkerboard artefacts and mesh bias, while the sparse polynomial-chaos expansion enables efficient estimation of low-order statistical moments required by the robust objective at a fraction of the cost of intrusive or Monte Carlo approaches. Numerical studies on a cantilever and a curved beam show that introducing non-Gaussian material variability leads to systematic load-path redistribution and a reallocation of 6-12% of the structural volume, together with a reduction in compliance scatter. The non-intrusive surrogate reproduces intrusive reference results within 3% using an order of magnitude fewer full finite-element analyses. These results demonstrate that the proposed framework offers a physically consistent and computationally efficient route to topology-optimized designs that remain reliable under realistic material uncertainty.

Understanding whether a digital marketing campaign is effective is central to designing effective customer engagement strategies. We analyze a large-scale, longitudinal promotional email campaign conducted by a U.S.\ retailer to evaluate how value-added incentives, such as free shipping, compare with traditional price discounts in influencing customer purchasing behavior. The analysis is complicated by non-compliance, due to not opening emails, multiple longitudinal mediators, and zero-inflated mediators and purchase outcomes. To address these challenges, we develop a Bayesian causal mediation framework based on enriched Dirichlet process mixture models and estimate the causal estimands using a scalable G-computation algorithm. We show that analyses ignoring email-opening behavior substantially attenuate estimated effects. Value-added incentives consistently outperform price discounts, yielding higher estimated potential purchase amounts, with benefits accumulating over time. We design an individualized sequential emailing strategy that optimizes expected purchase count in the observed data.

We propose a methodology based on the standard ReLU Deep Neural Networks (DNN) to make predictions and quantify their uncertainty. Classically, people rely on linear, non-linear, or non-parametric kernel methods to fit and then predict the time series. As the universal approximation ability was revealed for DNN, its application has become more and more popular for prediction tasks in various scientific areas. However, the corresponding uncertainty quantification has not been studied thoroughly. Particularly, the uncertainty in prediction will consist of two parts: (1) the future variability; (2) the estimation variability within training data. To capture both variabilities, we build the so-called pertinent prediction interval (PPI) with the DNN model estimator. We first explore the consistency property of the DNN estimator with beta-mixing dependent data. Subsequently, we show that the implied forward bootstrap series is still beta-mixing and possesses the same stationary distribution as the original time series in probability, which is a key condition to enable the PPI. Lastly, the desired PPI is built after imposing minimal conditions on the limiting distribution of predictive roots. Simulations and real-data analysis are deployed to challenge our approach with standard non-parametric methods.

Decorrelating a multivariate system need not destroy source-variable identity. We introduce Constrained Oblique Rotation with Anchored Loadings (CORAL), which minimizes residual cross-correlation while guaranteeing a declared minimum correlation between each transformed variable and its designated source. For a p-variable correlation matrix $R$, we show that every exact decorrelator can be written as $R^{-1/2}Q$ for some orthogonal matrix $Q$, and define $\rho_\star(R)$ as the maximum common source fidelity compatible with exact decorrelation. Constructive lower bounds and rigorous analytical upper bounds tightly bracket $\rho_\star$ at [0.972,0.976], [0.959,0.962], and [0.949,0.950] in simulations with p={6,18,50}, respectively, compared with PCA's largest achievable minimum correlation between distinct principal components and matched source variables of 0.358, 0.329, and 0.172. Corresponding intervals are [0.751,0.758] for World Development Indicators and [0.826,0.835] for wine chemistry data sets. Thus, loss of source-variable identity is not inherent to exact decorrelation but depends on the decorrelator selected. CORAL uses constrained Riemannian optimization and extends to exact support restrictions.

In this paper we show how to extend quantile-stratified sampling to produce simulations from various multivariate distributions. These simulations have desirable space-filling and coverage properties relative to simulation using IID sampling. We examine the coverage performance of these simulations against IID sampling by looking at plots of ordered log-density values from the simulations.

This paper presents a framework for constructing structure-preserving representations of complex systems through discrete approximation, and demonstrates its use in studying the vertical temperature and salinity structures in the mesopelagic zone across the global ocean using the ARGO dataset. Clustering serves as a means of organizing complexity into a finite set of structures that approximate the overall oceanic conditions, and a color encoding design then integrates these structures into a coherent map, with the three color components derived from interpretable geometric features of a profile: its initial level, its magnitude of variation, and its shape. Instead of focusing on specific depth levels or computing zonal averages within selected regions, our approach preserves the full vertical structure of individual profiles and incorporates each profile in the global ocean, capturing both fine-scale profile detail and large-scale spatial variability. By clustering over one million profiles collected over a decade, we identify and characterize representative profile shapes, which form the basis for a visualization strategy that provides an integrated, comprehensive, and interpretable presentation of the large-scale spatial distributions of these oceanic vertical patterns.

Three-player Auction Bridge is a finite imperfect-information game in which dynamic partnerships create a distinctive interaction between scoring, strategic incentives, and payoff distribution. This paper develops a unified statistical and game-theoretic framework to evaluate a traditional scoring rule against a modified mechanism designed to improve strategic incentives. We first identify a structural defect in the original scheme: a bid-invariant slam bonus can make lower contracts strictly more attractive than higher ones. A bid-dependent correction removes this incentive distortion. We then compare the two schemes using higher-moment analysis, fairness measures, Nash equilibrium analysis, and General-Sum Counterfactual Regret Minimization (GS-CFR). The results show significant differences in the shape of the payoff distributions and reveal that, although both schemes remain highly balanced across physical seats, the modified scheme substantially increases inequality across strategic roles. The game is shown to be neither zero-sum nor constant-sum, motivating a general-sum rather than minimax analysis. Full-game GS-CFR further indicates a substantial increase in the bidder's payoff under the modified scheme, accompanied by a reduction in defender payoff and a shift from separating to partially pooling bidding behavior. These findings demonstrate that scoring-rule design can fundamentally reshape incentives, payoff distribution, and information transmission in imperfect-information games. The proposed framework provides a systematic approach to evaluating such mechanisms from both statistical and strategic perspectives.

Area-level small area estimation (SAE) models, such as the Fay–Herriot (FH) model, borrow strength across domains through covariates and random effects, but they can struggle when the relationship between the covariates and the outcome is spatially heterogeneous, that is, when it changes across the spatial domain of interest. We propose a spatially-clustered FH (SC-FH) framework that simultaneously (i) estimates cluster-specific regression coefficients and random effects variances and (ii) generates spatially coherent partitions of the geographical domain. Estimation maximizes a penalized likelihood that augments the FH likelihood with a Potts-type spatial cohesion term over the areal adjacency graph, through an efficient strategy that alternates between sequential label updates and closed-form FH updates within clusters. In simulation experiments run on the real geography of the application, the method recovers the latent regimes almost exactly whenever they are separated in the covariate–response space and improves prediction accuracy over the standard FH benchmark, with the spatial penalty acting as a stabilizer of both classification and estimation. An empirical application to the average standard output of farms in the Po Valley (Northern Italy) identifies two spatially compact production regimes with significantly different cluster-wise coefficients, and shows that the clusterwise predictor improves on the direct estimates while avoiding the over-shrinkage of the pooled model.

We develop a new uniform drift condition and local minorization that implies a stronger weighted form of uniform ergodicity for Markov chains we call hyper-V uniform ergodicity. The convergence guarantees geometric decay of the bias towards the invariant measure independently of the initialization for all functions controlled by a dominating function V. A key advantage of the approach is that it bypasses the need to establish a global minorization condition, which is often substantially more difficult to verify in practice, while yielding stronger convergence guarantees than global minorization. Optimal convergence bounds in a minimax sense of the framework are established. The utility of the framework is demonstrated through applications to the P'olya-Gamma and Kolmogorov-Gamma Gibbs samplers. We also show qualitative hyper-V uniform ergodicity convergence for two-variable Gibbs samplers can be inferred by the form of the invariant measure, bypassing convergence analysis entirely.

Calibrating computational models to experimental data is a core task in applied statistics, especially in scientific domains, where physical experiments are costly and simulations play a central role in design and inference. Motivated by uncertainty quantification challenges in particle accelerator experiments, we develop and evaluate a Hierarchical Bayesian Calibration framework. In contrast to standard Bayesian calibration, certain inputs - such as beam injection amplitude - must be estimated separately for each experiment. We adopt the Kennedy-O'Hagan formulation and extend it with a hierarchical prior structure to model the distribution of experiment-specific calibration parameters, thus borrowing strength and improving generalisation across repeated experiments. A key methodological challenge arises from the need to evaluate a large number of forward simulations, which renders conventional Markov chain Monte Carlo approaches computationally prohibitive. To address this, we leverage the Bayesian Committee Machine as a scalable modelling strategy for Gaussian Process emulators. The BCM provides a principled divide-and-conquer approach, enabling parallel inference and reducing computational cost without requiring problem-specific tuning of the emulator approximation. Posterior sampling is performed using the No-U-Turn Sampler, supported by automatic differentiation in Julia, which removes the need for analytic gradient derivation and facilitates flexible model specification. We assess the proposed framework using established benchmark problems and simulated data from the Argonne Wakefield Accelerator. The results demonstrate substantial computational savings and robust calibration performance, highlighting the applicability of the method to large-scale scientific modelling problems.

We present a hierarchical Bayesian framework for non-homogeneous pairwise interaction Gibbs point process models, where the global and local effect functions are modeled via basis function expansions. We further propose a testing procedure in order to assess complete spatial randomness. The proposed methodology is exemplified through two real benchmark data examples involving water striders and forest fires.

Cramér's V, the effect size reported beside almost every chi-square test, is almost never accompanied by a confidence interval, because the interval is a hard nuisance-parameter problem: infinitely many tables share one effect-size value. The two intervals an analyst can reach for are unsatisfactory. The projection of a joint confidence region onto the effect size is guaranteed for every table but wide, over-covering at essentially 1.000 on larger tables; the noncentral inversion the software prints, and the bootstrap, are narrow but do not cover (median coverage 0.31 on this study's grid). This paper supplies an interval that is both valid and informative. Conditioning on the observed margins makes the exact conditional distribution of the Pearson statistic under a non-null association computable with no asymptotics and no Monte Carlo; inverting it yields a near-nominal confidence interval for Cramér's V that is a third to a half the width of the projection, the advantage growing with table size. The one price, stated plainly, is a change of estimand to the effect size at the observed margins; the mid-p default undercovers small effects, where a conservative variant or the projection remains the fallback.

We introduce a simple and performant approach for rapidly and accurately performing polyharmonic spline (PHS) interpolation using a combination of the fast multipole method (FMM) from computational electrostatics and the Vecchia approximation from the Gaussian process and sparse approximate inverse literatures. Using basic properties about Hadamard products and low-rank matrices, we demonstrate that an FMM with two kernels, the logarithmic and distance kernels, results in fast PHS interpolation for all orders. Furthermore, we demonstrate the exceptional performance of sparse inverse approximation methods with the Matérn covariance model for preconditioning. Combined with careful management of disallowed subspaces, we describe a procedure for obtaining prediction weights using preconditioned conjugate gradient that converges in less than $15$ iterations, even for problem sizes with over one million points. As a result, thin-plate spline interpolation—a particularly popular method that does not require parameter tuning—that matches the fully dense $\mathcal{O}(n^3)$ computation in accuracy can be done at the cost of approximately $50-60$ FMMs. A high-performance software library for odd-order PHS interpolation in two dimensions is made available as a companion to this work.

Omics datasets, particularly single-cell RNA sequencing data, are high-dimensional, sparse, noisy, and dominated by zero values, making faithful low-dimensional representation challenging. Existing dimensionality-reduction methods may distort local neighbourhoods, global organization, or the cohesion of meaningful populations, with similar limitations arising in genealogical data. We introduce Contrastive Manifold Approximation and Projection (CosMAP), a graph-based unsupervised dimensionality-reduction method for producing faithful and interpretable embeddings. CosMAP extends the graph-based framework of UMAP by combining cosine-similarity neighbourhoods with temperature-normalized contrastive affinities, which are optimized in the embedding space using an attractive–repulsive objective. It further employs a two-phase refinement strategy: an intermediate higher-dimensional representation is first learned and then used to reconstruct the neighbourhood graph and initialize the final low-dimensional embedding. We evaluate CosMAP on MNIST and USPS handwritten-digit datasets, mouse retina and cortex single-cell RNA-sequencing datasets, and a large genealogical kinship dataset derived from BALSAC-CARTaGENE. Compared with state-of-the-art dimensionality-reduction methods, CosMAP produces more coherent visual representations, improves neighbourhood preservation, and provides clearer global organization of digit classes, biological cell populations, and regional genealogical patterns. These results indicate that CosMAP offers a robust framework for exploratory analysis of complex, sparse, high-dimensional data. The implementation is publicly available at https://github.com/FenosoaRandrianjatovo/CosMAP-dr.

Optimal stratification aggregates strata into a small number of final strata to minimise total sample size required to meet target precision constraints. This combinatorial objective, reformulated as a within-cluster dispersion surrogate, can be expressed as a quadratic unconstrained binary optimisation (QUBO) problem. This paper reports a comparative case study of four solvers for that surrogate, run under heterogeneous free-tier constraints on an identical twenty-stratum frame from the swissmunicipalities dataset: a D-Wave-formulated QUBO, solved here by simulated annealing (QPU access unavailable); a gate-based quantum processor (IBM Quantum, running QAOA); a photonic entropy-quantum-computing device (QCi Dirac-3); and a classical GPU-based Ising machine as control (Fixstars Amplify AE). Solvers differ in hardware, computational budget, iteration count, and, for IBM, the encoding itself, so the comparison is a case study rather than a controlled experiment. Against a genetic-algorithm benchmark of 129 sample units, the best-known QUBO solution (objective 501.0) maps, after Bethel-Chromy evaluation, to a sample size of 160, a gap of roughly 24%, tied to the best solution found rather than a certified optimum. IBM and Dirac-3 fall further short, at 262 and 286; the IBM result is consistent with SWAP-routing degradation on a sparse qubit lattice, though not isolated from other causes, while Dirac-3, despite genuine optimisation, settles short of what a free GPU reaches in seconds. The paper concludes that a classical genetic algorithm remains, at the scale tested, the method of choice for optimal stratification, reflecting (a) limits of the dispersion surrogate versus the true Bethel-Chromy objective, not encodable as a low-degree polynomial, and (b) solver- and platform-specific limits in the IBM and Dirac-3 runs. These findings do not establish global optimality for either result.

High-dimensional kernel density estimation suffers from the curse of dimensionality. This study proposes a hybrid density estimator defined as the product of a joint density over a pre-specified subset of dimensions and marginal densities for the remaining variables, instead of estimating the full high-dimensional density directly. Under this framework, given the observed data, we select a subset of variables for joint density construction to improve estimation accuracy, while modeling the remaining variables through their respective marginal densities. This construction involves a trade-off between the approximation error induced by simplifying the dependence structure and the variance reduction achieved by lowering the dimension of the joint density. We employ a genetic algorithm to efficiently identify such variable subsets. Theoretical and numerical results demonstrate that the proposed estimator can outperform conventional full-dimensional kernel density estimation when this trade-off is appropriately balanced.

Quantile regression is well suited to heterogeneous and heavy-tailed data, but computation becomes challenging for large, distributed data sets because the check loss is nonsmooth. We propose a distributed stochastic smoothing alternating direction method of multipliers (DSS-ADMM) for horizontally partitioned penalized quantile regression. Each worker computes a mini-batch gradient of a Huber-smoothed check loss, and a coordinator performs a single proximal aggregation step for the regularizer. Raw observations remain local, worker updates run in parallel, and the method requires no matrix inversion. For proper, closed, and convex penalties, stacking the local coefficient vectors yields a standard two-block stochastic ADMM formulation. With fixed smoothing, we establish an expected $O(\log K/\sqrt K)$ joint objective-feasibility bound and an explicit $\eps/4$ approximation term for the original check-loss objective; when the smooth block is strongly convex, the bound improves to $O(\log K/K)$. We also characterize the scope of an extension to the minimax concave penalty and the smoothly clipped absolute deviation penalty. Reproducible simulations consider both homogeneous worker partitions, in which observations are independently and identically distributed across workers, and heterogeneous partitions, in which worker-specific covariate distributions differ. Sensitivity studies and analyses of the diabetes and Engel data illustrate the trade-offs among per-observation gradient evaluations, communication, consensus, sparsity, and prediction.

Hamiltonian Monte Carlo (HMC) is a successful generic inference method in probabilistic programming, but in its ordinary formulation it needs gradients and finite-dimensional parameter spaces. In Haskell, lazy evaluation lets probabilistic programs express stochastic processes and other non-parametric Bayesian models over implicit infinite-dimensional spaces. This paper develops new formulations of gradient-based HMC for this infinite-dimensional setting, via lazy evaluation. For automatic differentiation, we provide an analysis based on a new notion of "piecewise analytic under cylindrical analytic partition" (PACAP), to show that even if a program is infinite-dimensional and defined lazily, the gradient of the likelihood function is finitely supported. For the Monte Carlo method itself, we develop several HMC variants and a No-U-Turn Sampler that operate over the infinite-dimensional parameter space but are still productive because of lazy evaluation. Experiments cover Gaussian mixture clustering, random walks, and piecewise-constant regression with Poisson-process changepoints.

Every chi-square-family statistic for a two-way contingency table (Pearson's X^2, the power-divergence members, the variance-stabilized T_root) is referred to an approximate null distribution (chi-square, moment-matched chi-square, saddlepoint, or bootstrap) that miscalibrates on sparse or heterogeneous tables, where the exact distribution is a lattice of atoms rather than a continuous curve. The exact reference conditions on both margins (the multivariate Fisher noncentral hypergeometric law) but is usually treated as uncomputable, because the number of tables sharing a margin is astronomical. We show that the exact conditional moments and the exact conditional distribution of any chi-square-family statistic are computable without enumerating tables, because the statistic is additive over cells and the margin-conditional law factorizes cell by cell: a dynamic program walks the table one cell at a time, carrying per row-capacity state either a few moment accumulators, a map from statistic value to probability, or one complex number (the characteristic function), at a cost set by the number of margin states rather than the number of tables. The moment engine computes the exact conditional moments of a 5x5 table at three per cell (about 79 billion tables on the margin fibre) in three seconds; the distribution engine returns the exact tail, verified against enumeration to machine precision, where a three-moment chi-square fit misses it by up to 0.28. The engine furnishes Monte-Carlo-free ground truth against which any approximate reference can be measured, extends to the non-null alternative, and is the computational substrate beneath the exact conditional interval and reporting tools of the companion papers.

Phenotyping conformation traits is important for dairy cattle breeding and management. Although large-scale phenotyping is possible with 3D imaging technologies, manual annotation of anatomical landmarks and long run times prevent full pipeline automation. In particular, the morphological heterogeneity of cow udders makes automated detection of teat landmarks challenging. To address this limitation, we propose and evaluate a method for teat length estimation from udder point clouds with a Gaussian process. We model the vertical coordinates as the sum of a Gaussian process representing the udder floor and an unknown function representing the teat. Since the udder floor process is smooth and has a significantly wider dependence lengthscale than the teat function, this model allows separating the two terms needed for teat landmark definition. To ensure computational feasibility, we implement a low-rank approximation of the covariance matrix, reducing the computational complexity of the method from $\mathcal{O}(n^3)$ to $\mathcal{O}(n)$. This approach is both faster and more accurate than existing methods, reducing RMSE by a factor of two. It is also more robust to uncommon udder morphologies, making it better suited for automated phenotyping of large numbers of individuals.