Papers

We propose sliced Orlicz-Wasserstein (SOW) distance which is a generalization of sliced Wasserstein (SW) distance. SOW replaces the $L^p$ norm in SW with a Luxemburg norm cost induced by an Orlicz function $\phi$. First, we prove that SOW distance is a metric on the space of measures with finite Orlicz norm, and show that it recovers the SW distance when the Orlicz function is $\phi(x)=x^p$. Next, we derive the topological properties of the SOW distance. In particular, we show that convergence under SOW implies weak convergence, and the converse is true under the compact support condition. We then present the theoretical results for estimating the SOW distance. We derive sample complexity for both the distance itself and the powered functional of the distance, and prove their minimax optimality. In addition, we discuss the computational algorithm for approximating the SOW distance by Monte-Carlo estimation and bisection search, as well as the associated approximation error and computational complexity analysis. Our experimental results reveal the superior computational efficiency of SOW compared with Orlicz-Wasserstein (OW) distance. Also, in the experiments, we demonstrate the favorable flexibility of SOW distance over SW in detecting differences between distributions by comparing their performance in two-sample tests and evaluating generative models on image datasets.

Amortized neural likelihoods enable computationally expensive inference for models with analytically intractable or unspecified likelihoods, but their black-box nature limits interpretability. We introduce a symbolic distillation pipeline that converts trained neural likelihoods into explicit, interpretable expressions optimized for efficient parameter estimation. Our approach uses a recovery-directed objective to guide symbolic regression toward expressions that preserve parameter-recovery accuracy rather than merely approximating the likelihood function. Candidate expressions are evaluated on held-out datasets and selected using a criterion that jointly accounts for expression complexity, parameter-recovery performance, and distributional distance from the learned likelihood. We evaluate the pipeline on the diffusion decision model, a classical cognitive model, whose analytically tractable likelihood provides ground truth for controlled evaluation. The proposed recovery-directed objective improves parameter recovery over standard symbolic-regression objectives. The resulting symbolic likelihoods enable over 100 times faster parameter evaluation than both neural likelihoods and, when available, the exact likelihood, while maintaining a manageable loss in precision. We further demonstrate these computational benefits in Bayesian hierarchical inference on empirical data. Our pipeline provides a lightweight interface for integrating symbolic distillation with existing neural-likelihood estimation methods and can be adapted to a range of simulation-based inference settings.

Nested sampling estimates the evidence by assigning prior volumes to an ordered sequence of likelihood contours. Markov chain constrained-prior samplers, which are used in high-dimensional problems, generate many intermediate states before producing the next classic sample. These intermediate states are discarded because their correlation prevents them from being inserted into the ordered NS sequence without changing its order-statistic law. This paper introduces a novel method of using them to improve evidence estimation, by formulating NS in a Bayesian way and conditioning on phantom samples as Monte Carlo observations. We then present the open-source software package, JAXNS v3, and its implementation choices. We validate the approach on a set of problems, and identify its limitations via ablation. In our experiments, when problem structure is well resolved, conditioning on all retained phantom samples reduces log-evidence RMSE by at least $30\%$, with larger improvements at higher dimensionality. For the tested problems with unresolved structure, full phantom conditioning produces no detectable improvement or deterioration in evidence accuracy. Phantom conditioning produces overconfident evidence uncertainties. We also introduce two dynamic nested sampling allocation schemes. Evidence-improving allocation approximately halves the number of likelihood evaluations required to achieve comparable evidence accuracy relative to uniform allocation. Posterior-improving allocation doubles the classic posterior's Kish effective sample size for $17.5\%$ additional likelihood evaluations.

Deep learning has substantially accelerated the calibration of complex stochastic-volatility models, but neural point calibration alone does not capture the uncertainty remaining after an implied-volatility (IV) surface has been observed. We develop a simulation-based inference framework for rough Heston (rHeston) calibration that learns the posterior distribution of the model parameters conditional on an IV surface. Using neural ratio estimation, we obtain calibrated posterior samples that can be propagated through heteroscedastic neural surrogate pricers for path-dependent exotic options. The resulting posterior-predictive distributions combine residual parameter uncertainty with conditional surrogate uncertainty and yield uncertainty-aware price intervals. We further introduce Hellinger-SHAP, an information-theoretic explainability method for posterior inference. Rather than attributing a single parameter point estimate, it applies local-background Kernel SHAP to a posterior-information functional measuring contraction from the prior to the posterior. This identifies maturity–moneyness regions associated with posterior information gain for individual rHeston parameters. In a simulation study, posterior-predictive intervals provide calibrated or conservative coverage across forward-start, barrier, and realized-variance claims, while point plug-in prices can be materially unreliable for selected contract regimes. Together, the UQ and XAI analyses provide a transparent framework for uncertainty-aware neural calibration and downstream exotic pricing under the specified prior-predictive model.

Tree-Adjoining Grammars (TAGs) have recently been introduced to Nonlinear System Identification (NLSI) as a means of encoding an entire model class as a finite set of grammatical rules, from which candidate models are assembled as trees. Existing TAG-based identifiers rely on evolutionary optimisation and return point estimates of the model structure. This paper instead proposes the TAG framework within a Bayesian setting. A generative prior is defined over tree structures and their parameters, and a Reversible-Jump MCMC sampler with structure-preserving tree moves is used to infer the joint posterior over model structure, parameters and predictions. Two training objectives are considered; that is, a one-step-ahead objective with conjugate parameter proposals, and a simulation-based objective handled by likelihood-free inference. The approach is validated on a simulated polynomial NARX system, the Silverbox benchmark, and wave-loading data from the Christchurch Bay Tower, where embedding Morison's equation as a fixed initial tree yields a grey-box model that outperforms the physics-driven baseline. The results demonstrate that Bayesian TAGs are well suited to quantifying uncertainty in equation discovery for dynamical systems and to fitting physics-informed models.

Many studies expose units to several binary treatments at once and record several correlated outcomes, yet analysts usually estimate one treatment's average effect on one outcome at a time – discarding how treatments interact and how their effects vary across units. We introduce the factorial multivariate Bayesian causal forest, a Bayesian nonparametric model that decomposes the factorial response surface over the treatment lattice into a prognostic sum-of-trees plus one sum-of-trees per main or interaction effect, each with correlated multivariate leaf parameters, estimated jointly with coherent uncertainty. A single indicator-weighted kernel samples every component – the prognostic term is the special case whose indicator is unity – so the two-treatment model, the single-treatment multivariate causal forest, and a general order-r truncation are one and the same sampler. We give the ANOVA/Möbius identification of each estimand, an efficient Rcpp engine, and an interpretability layer using value-suppressing uncertainty maps. In simulations the average-effect estimators are unbiased with nominal coverage, consistent, and robust under misspecification, failing only under unmeasured confounding, which we flag. We illustrate the method on clinical, agricultural and economic data and on a deeper application to the NHANES survey. Software is provided as an R package.

Trimming methods for robust clusterwise regression discard a fixed fraction of the data. Too low a level breaks the fit; too generous a level can trim away a small group. We first propose a group-recovery step that can follow any trimming or flagging method: it searches the discarded units for a line, tests whether those near it form a peak rather than a band, and restores the line as a group when the likelihood of Gaussian groups plus uniform noise improves by a margin like that of the Bayesian information criterion. In simulations it repaired the failures of a generous TCLUST-REG level with unequal groups, but not with three or four groups. The second proposal, ESF (exact-subsample flagging), is a flagging procedure without a trimming level: it solves the clusterwise least-squares problem exactly on small subsamples, flags units far from the best fit, and draws later subsamples from the rest. Two constants stand in for the level: a subsample size, set from a lower bound on the smallest group proportion, and a cap on the flagged set. It is meant for data about whose contamination nothing is known: TCLUST-REG at the fixed level 0.30, followed by reweighting and the recovery step, was as accurate as ESF on average up to a fifth of outliers, and higher levels were more accurate beyond. The flagged fraction estimates the contamination only when the errors are close to Gaussian. On taxi fares with known tariffs, ESF found both in every sample.

We consider stochastic nonlinear optimization problems with deterministic equality constraints. While unconstrained stochastic optimization is well understood, the interplay between optimality and feasibility in the constrained setting poses significant challenges. Moreover, existing theoretical guarantees for constrained stochastic methods predominantly rely on bounded-variance assumptions, leaving the heavy-tailed noise regime largely unexplored. To address this gap, we propose a novel trust-region method within the stochastic sequential quadratic programming framework, termed TR-SSQP. Our method employs a normal-tangential decomposition in the step computation to balance optimality and feasibility. In addition, we incorporate a normalization mechanism in the design of the trust-region radius, together with Polyak momentum for gradient estimation, ensuring stable updates without gradient clipping. When the trust-region radius and the momentum parameter decay at appropriate rates, we establish global almost-sure convergence of the method. To the best of our knowledge, this is the first asymptotic convergence result for constrained stochastic optimization under heavy-tailed noise. We demonstrate the promising performance of the proposed method through extensive numerical experiments, including comparisons among its variants and with existing constrained stochastic optimization methods.

We consider the problem of learning the underlying causal directed acyclic graph (DAG) structure corresponding to a structural equation model (SEM) with non-Gaussian errors. Motivated by an intentionally misspecified non-Gaussian SEM with all Laplace errors, we first introduce the mean absolute residual risk, defined over the space of all real matrices, and show that, asymptotically, the risk of the true weighted causal DAG matrix is strictly smaller than that of any other matrix. Nevertheless, to enhance generality and account for high-dimensional and finite-sample settings, we further incorporate row-specific sparsity penalties along with a soft DAG constraint to derive a continuous score function over the space of real matrices. Accordingly, we propose a score-based DAG learning method, named MARCEDES, formulated as an unconstrained score minimization problem, which can be efficiently solved using gradient-based optimization techniques, thereby circumventing the challenges associated with constrained optimization. Furthermore, we develop a computational algorithm to handle the non-smoothness of the score objective and to enable optimal tuning of row-specific sparsity penalties under a generalized Bayes framework. Finally, we demonstrate the efficiency and improved performance of the proposed method over existing approaches through an extensive simulation study.

Readable annotations are essential for interpreting statistical graphics, yet conventional single-color lines and labels can lose visibility when they cross backgrounds with heterogeneous luminance. We introduce a contrast-aware annotation framework implemented in the ggtwotone R package for ggplot2. The framework combines dual-stroke rendering, adaptive text-color selection, and perceptually guided highlight palettes to improve annotation visibility across light and dark regions. A shared contrast-adjustment engine supports WCAG- and APCA-based contrast criteria, reducing the need for manual color adjustment. The package provides contrast-aware geoms for segments, curves, paths, mathematical functions, regression overlays, and text while remaining compatible with standard ggplot2 workflows. A simulation-based evaluation across heterogeneous background colors demonstrates improved worst-case contrast for dual-stroke annotations and adaptive text selection, while also identifying trade-offs in perceptual separation as the number of requested highlight colors increases. Applications to statistical graphics and scientific images illustrate the framework in practical visualization settings. Together, these tools provide a reproducible approach for incorporating contrast considerations directly into graphical annotation.

In this A.I. era, Computational Taxonomy is proposed to study complex systems by analyzing their databases under taxonomic hierarchies abiding the Principle of Science by providing "good explanations". Comparisons among branches or classes are carried out by Scientific Data Analysis (SDA) paradigm that explores all potential associative patterns, including interacting effects of all high orders, and evaluate finite sample precisions for all information pieces individually by effectively making use of all variables' categorical nature. Under each comparison, all confirmed information pieces are collected and displayed along row-axis of a heatmap with all involved study-subjects on the column-axis. Each comparison's heatmap individually characterizes participating classes and study-subjects and simultaneously provides a scientific basis for outlier detection upon all non-participants. All these heatmaps then collectively constitutes so-called Class-informatics that offers good explanations based on characteristic of all classes and study-subjects. Computational Taxonomy's Class-informatics indeed resolves multiple fundamental issues: Tukey's more than 60 years outlier detection problem, issue of self-correction annotation, and a crucial check on assumption of information-content equality between testing and training data sets in Machine Learning. A showcase of Computational Taxonomy is exclusively illustrated on Iris data.

The evidence behind a scientific or operational decision is often held by hospitals, banks, or registries that cannot pool individual observations. Cross-silo federated learning moves computation to the data and exchanges agreed summaries. Targeted maximum likelihood estimation (TMLE) refines a flexible initial fit, yielding plug-in estimators that respect the model and support efficient inference. TMLE itself, however, has remained a fully centralized procedure. To fill this gap, our paper introduces the first federated TMLE algorithm. We federate targeting itself, for an arbitrary target, loss, and fluctuation family, through two complementary frameworks. FedTMLE-G aggregates local gradients and reproduces centralized targeting step for step. FedTMLE-L lets each institution complete its own fluctuation fit before a single exchange of fitted updates, trading synchronized fidelity for local autonomy. For gradient aggregation, we develop a finite-precision protocol that transmits changes rather than values and certifies targeting accuracy within explicit bounds on exchanges and bits. A description-length analysis of the accepted updates then shows that this finite communication leaves numerical targeting error negligible against sampling uncertainty. The cost of computing an estimator is thus distinct from the complexity of selecting it. Our analysis also indicates that keeping data local is not itself a privacy guarantee of TMLE, since instability of full-record reconstruction need not prevent recovery of a specified sensitive attribute. For a personalized version of local averaging, institutions retain their own estimates and leave once local targeting is complete. A nonconvex convergence bound charges the improvement forfeited through averaging to disagreement among local fits and exposes a tradeoff between equal institutional influence and the sampling variability of small silos.

Thin plate splines are highly attractive smoothers. However, they have cubic computational cost, which severely limits their use in practice. As a remedy, Wood (2003) suggested thin plate regression splines (TPRS), which provide a low rank approximation. The key step of the TPRS approximation is a truncated eigendecomposition of the radial basis function (RBF) design matrix. However, as Wood (2003) writes, the optimality of the TPRS approximation is a slightly weak one. This is because the RBF coefficients are subject to orthogonality constraints and the TPRS approximation is only optimal if these constraints are ignored. To address this shortcoming, we suggest a slightly different low rank approximation. The suggested approximation is based on a truncated Demmler-Reinsch basis (TDRB), which provides a best low rank approximation of the smoother matrix in terms of Frobenius and spectral norm. We prove that the TDRB smoother achieves the optimal rate of convergence and suggest an efficient algorithm for its construction. This algorithm is based on a truncated Karhunen-Loève (KL) expansion of the equivalent Bayesian smoothness prior and it has the same computational cost as required for TPRS. We demonstrate the applicabilty of our approach through simulations and a real data example. We find that the performance is very similar to that of TPRS but the suggested approach has some advantages.

Sparse direct Cholesky solvers fix one data structure for an entire matrix, but symmetric positive definite systems range from nearly dense to irregular, sometimes mixing both within one matrix. We let the data structure follow the sparsity structure, across matrices and across tiles within a matrix. Before numerical work starts, a lightweight selector captures the sparsity pattern of the Cholesky factor and routes the matrix, on one static shared-memory schedule, to one of three regimes: dense, sparse, or an intermediate semisparse regime. Dense tiles are stored in full; the semisparse regime uses a new active-column tile that keeps only the columns the factorization will touch, so banded or arrowhead-shaped structures common in spatial models still reach BLAS-3 efficiency. The approach suits the integrated nested Laplace approximation (INLA): one sparsity pattern factorized thousands of times with different values, so the one-time analysis cost is amortized and every factorization saving compounds. We evaluate the technique on 60 SPD matrices, comparing against MUMPS, PaStiX, CHOLMOD, symPACK, and Intel oneMKL PARDISO on Intel Xeon and AMD EPYC nodes; the selector alone achieves the lowest total factorization time in every regime. Summed over the suite, it beats the best fixed single-structure mode by 1.6 to 2.6x, and every alternative by 1.8 to 12.5x on Intel and 2.6 to 10.1x on AMD, with the largest gains on the most expensive factorizations. It trades more one-time analysis for less time per factorization, pulling ahead by the third factorization of a given pattern. As a first GPU extension, the dense route on one NVIDIA A100, with the factor resident on the device, runs 1.2 to 6.3x faster than on the faster CPU node, the margin widening with factor size. Solver, Python/R/Julia interfaces, benchmark suite, and results are open at https://github.com/esmail-abdulfattah/sTiles.

We introduce GRACE (the General Regime-Aware Changepoint Estimator), a general, globally exact multiple-changepoint detection and estimation framework in which arbitrary numbers of model parameters may change at or persist across estimated changepoint locations, and all permissible changepoint types and patterns of parameter sharing compete within a single $\ell_0$-penalized optimization. The analyst specifies which types of changepoints are permitted and which parameters each type must preserve, without specifying which changes actually occur, where they occur, or in what order. GRACE jointly estimates the number and locations of changepoints, the type of change at each location, the model fitted between changepoints, and all continuous parameters. As a direct consequence, GRACE can estimate models in which the slope persists across a mean shift. The generality of the framework is achieved through an exact dynamic-programming recursion that retains the accumulated fitting cost conditional on whichever parameters a subsequent change may preserve, tracking multiple parameters jointly when necessary. We establish the global exactness of the recursion and develop functional, penalty-based, and objective-bound pruning rules that eliminate candidates without sacrificing optimality. Ordinary optimal partitioning, FPOP, and CPOP are recovered as special cases, while piecewise-polynomial and seasonal models illustrate the broader range of permissible segment structures and transitions. The resulting framework supports heterogeneous changepoint types, parameter sharing across segments, and substantial parallel computation within a unified approach to exact changepoint estimation.

Local geometric mixing localizes geometric mixing by requiring geometric convergence to equilibrium in total variation only over finitely many transitions. It accommodates local convergence rates and captures rapid local equilibration, even when global mixing is much slower. We establish and discuss local geometric mixing bounds through Dobrushin contraction. We then apply this approach to Diffusion Path Monte Carlo, a recently proposed Markov chain Monte Carlo method, aimed at leveraging advances in score-based modeling, whose ideal transitions coincide with those of the Proximal Sampler. Our analysis covers both the ideal method and its implementable Metropolis-adjusted counterpart, providing mixing guarantees under minimal assumptions. For the ideal method, these guarantees complement recent spectral gap estimates, which we develop into mixing time bounds.

We employ Computational Taxonomy (CT) to reanalyze the penguin data set penguins_lter by validating and addressing two biological issues: Sexual Size Dimorphism (SSD) and mate-selection criteria. Via Scientific Data Analysis (SDA) computing, CT constructs a Taxonomic Hierarchy by splitting Species first and then Sex, without involving Island, to achieve less complexity. This Taxonomic Hierarchy validates SSD as a branch comparison: (Species, Sex = Male)-vs-(Species, Sex = Female), upon which SDA explores all potential pieces of associative information from all covariate feature-sets, including interacting effects from order-2 to order-4, and then confirms them via their idiosyncratic reliability checks. The collective of confirmed information pieces are displayed on a heatmap platform to manifest underlying dynamics of SSD with explicit block-structured heterogeneity found within males and females. SSD dynamics is explained through mechanistic dependence pertaining to one chief factor consisting of up to 8 feature-sets: Body-Mass coupled by combinations of {Culmen-length,Culmen-depth, Flipper-length}, and two minor factors consisting of low-order combinations of {Culmen-length,Culmen-depth, Flipper-length}. Such Intra-Sex heterogeneity invalidates all Logistic regression modeling on SSD in the original paper. Further, we explore potential mate-selection criteria through the data-frame of Nest-ID within-species homogeneity.

Bayesian Federated Inference (BFI) estimates statistical models from multicenter data when individual-level observations cannot be combined across centers. We present \pkg{BFI}, an \proglang{R} package that implements this methodology for Gaussian, binomial logistic, and survival regression models. Each center performs a Bayesian maximum a posteriori analysis and sends only parameter estimates and curvature information to a central server, where these summaries are combined to approximate the analysis of the combined data. The package supports prior specification, structured forms of between-center heterogeneity, several parametric and flexible baseline-hazard models for survival analysis, and treatment-effect estimation for observational and randomized studies. We describe the software design and the information exchanged between centers, give reproducible workflows for these analyses, and compare the federated results with pooled-data analyses for Gaussian, logistic, and survival models.

We consider $K$-dimensional Bayesian vector autoregressions (BVARs) with Cholesky stochastic volatility (SV), in which the innovation covariance matrix is a lower-triangular linear transform of $K$ independent univariate SV processes. Such models are widely used in empirical macroeconomics to capture time-varying uncertainty and improve forecast accuracy, but the cost of posterior simulation is the binding constraint on the size of the system, and a major bottleneck for empirical work. We introduce a Markov chain Monte Carlo (MCMC) kernel that mixes better than existing samplers at the same computational complexity. It rests on a reparametrisation that makes the $K$ volatility trajectories conditionally independent, and on Particle Gibbs to update each trajectory. The kernel targets the exact posterior, rather than an approximation of it, and in an application with $K = 15$ it raises the mean effective sample size per second, relative to the benchmark corrected triangular algorithm, by a factor of approximately $14$ for the VAR coefficients and $3.4$ for the volatilities. We also introduce a a Sequential Monte Carlo squared (\smcsq{}) sampler, which used our MCMC kernel as a building block, and which delivers at every $t$ the one-step-ahead predictive density and the marginal likelihood of the data up to $t$, and hence on-line forecasting and model choice. To our knowledge, this is the first algorithm that delivers sequential marginal likelihoods for Cholesky-SV BVARs with static contemporaneous coefficients and a non-conjugate prior. We illustrate both on US monthly macroeconomic data, with $K=15$ for posterior inference and $K=6$ for the sequential sampler and model choice.

This paper contains new 80-digit positive-weight Gauss-Hermite and positive-weight interior-node Gauss-Legendre quadrature rules for up to five dimensions and varying polynomial degree accuracy (depending on quadrature type and dimension). Some of these rules improve on the best available rules in the literature and some offer rules where none (other than the tensor product) existed. The results were produced by combining methodology developed by previous researchers with new approaches. The full precision rules themselves are at https://doi.org/10.5281/zenodo.22881864. Software in Julia, Python, and R providing the rules is available via a package registry and/or GitHub: Quadriceps.jl (both 64-bit and 128-bit; https://github.com/NittanyLion/Quadriceps.jl), quadriceps-py (64-bit only; https://github.com/NittanyLion/quadriceps-py), and quadriceps-r (64-bit only; https://github.com/NittanyLion/quadriceps-r). The software used to create these rules is available via PositiveWeightQuadratureSolvers.jl (https://github.com/NittanyLion/PositiveWeightQuadratureSolvers.jl). The replication package is at QuadricepsReplicationPackage.jl (https://github.com/NittanyLion/QuadricepsReplicationPackage.jl). A snapshot of all five packages is archived at https://doi.org/10.5281/zenodo.22883240.

We develop an exact simulation algorithm for scalar diffusion paths and diffusion bridges when the Poisson potential is unbounded in both tails. The method reconstructs the realized range of a Brownian bridge proposal by sampling its maximum and location, together with the maxima and locations of the two adjacent restricted Brownian meanders. Conditional on this finite information, the remaining path decomposes into four conditionally independent interval-constrained Brownian bridges, which can be sampled exactly at the Poisson times required by the rejection test. In contrast to constructions based on an enclosing range layer, the proposed representation retains the exact extrema and their locations. Our algorithm returns an exact finite-dimensional skeleton without time-discretization error and permits exact post-acceptance refinement at arbitrary finite collections of times. Numerical experiments validate the resulting finite-dimensional laws and identify the restricted-meander extremum simulation as the principal computational cost in the nonlinear example.

Principal component analysis (PCA) is a fundamental tool to reduce the dimensionality of the data in many applications. PCA finds a few signal directions that contain most of the variability of the data by computing the eigenvectors of the sample covariance matrix. In this work, we focus on the spiked covariance model, in which the data vectors are defined by a few orthogonal signals plus an isotropic Gaussian noise, and our goal is to estimate one or more of the leading signals. Our main theoretical finding is that the subspace spanned by several leading eigenvectors of the sample covariance matrix contains significant information about the desired signals long before the individual eigenvectors converge to the population principal components. To prove this, we derive a posteriori bounds for the angle between the subspace spanned by the desired population signals and the subspace obtained from the sample using perturbation theory for singular vectors. This leads to a new algorithm, SuperPCA (SUbsPace subsamplER PCA), which capitalizes on an approximate eigenspace of the sample covariance matrix to find the leading signals far more efficiently and accurately than classical PCA in the high-dimensional, multi-signal setting. SuperPCA exploits only a small number of subsampled coordinates of the data, which can lead to tremendous savings in data acquisition cost, especially when the signals are approximately sparse. For the same number of measurements, SuperPCA can offer a factor $10$ improvement in accuracy compared to the classical PCA method.

Physical experiments frequently result in inverse problems, where the recorded response is used to deduce information about the sample, typically yielding only partial recovery. This work presents a reverse Monte Carlo modeling method for solving such problems involving statistical averaging, e.g. at the nanoscale. The proposed algorithm operates in an input subspace of the forward problem, with dimensions reduced from those of the original input space. The dimensionality reduction focuses on features that influence the response, as irrelevant information is lost in the process. We first study the performance of the method in tailored test cases that also reveal information loss when several input characteristics produce similar responses. Finally, we apply the method to a previously published simulated dataset of X-ray absorption, X-ray emission, and X-ray photoelectron spectra of the water molecule.

Essential Data Points (EDPs) - the vertices of the convex hull of a bilinear data matrix $\mathbf{D}$ in its row space, column space, or both - are widely used in chemometrics to reduce the size of large data sets while nominally preserving their underlying geometric structure. Using simulated three- and two-component chromatographic/spectral data sets and a real source-apportionment data set ($\mathbf{D = C\, A^\mathsf{T}}$), together with Borgen-Rajkó plots, Procrustes analysis, and variance-covariance comparisons, we show that this preservation is only partial: row-wise EDP reduction preserves the row-space geometry (inner and outer polygons) exactly while distorting the column-space geometry, and column-wise reduction shows the opposite behavior; joint row-and-column reduction distorts both. We then give a rigorous, general proof - based on the four fundamental subspaces of a matrix and its singular value decomposition $\mathbf{D = U\,S\,V^\mathsf{T}}$ - that the subspace which is not being reduced is always preserved exactly, up to an orthogonal rotation, whereas the subspace whose ambient dimension shrinks is related to the original only through a general, non-orthogonal isomorphism. This distinction is confirmed numerically to machine precision ($\sim 10^{-14}$-$10^{-16}$) on the real data set, and a deliberate negative control confirms that the "ambient-shrinking" map is genuinely non-orthogonal (residual $\approx 1$). These results demonstrate that the apparent rotation of an EDP-reduced polygon relative to the original is not, in general, a rigid rotation, and that visual or numerical comparisons between an EDP-reduced data set and the original data require an explicit, mode-dependent change-of-basis correction before any geometric or statistical conclusion can be drawn. A MATLAB implementation of this correction is provided.

Large language models are increasingly deployed as course-specific tutors, but their usefulness depends on grounding in vetted instructional materials that are often revised mid-semester. Our prior work built a multimodal retrieval-augmented generation (RAG) system over an authentic machine learning course corpus (Foundations of Machine Learning) and found that retrieval improved contextual grounding, but that fixed retrieval strategies were suboptimal. That motivates a different question: whether how a corpus is structured at ingest time matters more than how much is retrieved at query time. We present a controlled head-to-head comparison of two knowledge representations over an identical classroom corpus: (A) vector RAG, replicating the best-performing configuration from our prior study, and (B) an LLM-compiled wiki (Karpathy framework), in which the corpus is synthesized at ingest into linked concept pages with explicit cross-references and citations back to source materials. We evaluate 59 questions spanning single-fact recall, cross-unit concept linking, synthesis and explanation, and currency after a syllabus revision, scored by an LLM judge against a human-authored rubric. Both representations answered single-fact questions about equally well (9.33 vs. 9.96 of 10), but diverged sharply on questions requiring links across course units. The compiled wiki remained accurate and grounded (9.93; 100% grounded in cited sources), while retrieval scored lower and was markedly less grounded (8.14; 64%). The wiki's citations let students and instructors trace any claim back to the lecture that introduced it, adding a layer of dynamic retrieval that machine learning courses require. While further testing is needed, instructors using AI to support learning in ML courses should consider wiki-based structure for its potential to support foundational elements of best practice.

Bias due to uncontrolled confounding can be analyzed using several external adjustment formulas. Like other bias expressions, the risk difference bias formula involving an uncontrolled polytomous confounder requires only three externally obtainable parameters: confounder prevalences among the exposed and unexposed and the confounder-outcome risk difference. This paper presents two-parameter bounds on the risk difference bias expression when only one of the three required bias parameters is unavailable. Each bound is sharp over the full range of the unknown parameter. Some numerical examples with graphs are provided.

Calibrating a simulation model involves estimating its parameters by comparing model outputs with experimental data, so that simulation results faithfully reproduce the experimental observations. When the outputs are functions of time, there are multiple ways to quantify the discrepancy between experimental and simulated curves. A recent approach based on elastic functional data analysis decomposes a functional output into two components: a function temporally aligned to a template, and the corresponding warping function. This decomposition splits the problem into two independent calibration tasks, thereby addressing functional misalignment. However, it assumes that experimental and simulated curves share the same temporal support, an assumption often violated in practice when initial or end times are themselves uncertain or depend on the calibration parameters. In this work, we reinterpret the decomposition step as an approximation to a more general Bayesian calibration problem that incorporates an error term on the time axis. This perspective allows us to naturally extend the framework to a broader family of time warpings with varying initial or end times, using partial elastic alignment. We illustrate the method on a synthetic test case, comparing it with existing Bayesian calibration methods and demonstrating improved surrogate performance and error modeling. We then apply the proposed approach to the calibration of an equation of state (a thermodynamic equation relating the state variables of a material).

Comparisons of confounding-control methods are informative only when simulated data contain the treatment-covariate relationship those methods address. We assessed whether Plasmode 0.1.0 preserves that relationship, uses the returned exposure to generate outcomes, and reports the effect implied by its generating model. We inspected the source and ran the generator under identical seeds with source rows as supplied, sorted by observed exposure, or using an alignment correction. We ran 200 null-effect replicates per strategy in a synthetic cohort and the SUPPORT/Right Heart Catheterisation cohort. A sampled subject usually received another subject's exposure probability. In the synthetic cohort, source-model AUC was 0.503 before and 0.841 after alignment. Crude risk-difference bias moved from +0.0051 to +0.3037. Other orders preserved or reversed the association. After alignment, deliberately misspecified estimators retained risk-difference biases of 0.10-0.11. Correctly specified estimators were nearly unbiased. With the unmodified generator, all five appeared similarly unbiased. At an injected odds ratio of 2, the conditional log odds ratio appeared on the observed exposure (+0.696; target 0.693), not the returned exposure (+0.017). Misalignment can weaken, preserve, or reverse designed confounding and thereby distort method comparisons. Alignment correction does not repair outcome-exposure decoupling, disagreement between the reported and generating effects, or marginal miscalibration.

We develop a stochastic nested fixed point (SNFP) estimator for random coefficients logit demand models that updates model parameters using stochastic gradients and performs demand inversion one market at a time. Relative to the conventional nested fixed point (NFP) estimator, SNFP substantially reduces memory requirements and computational cost, making estimation feasible in very large datasets. We establish the large-$T$ (number of markets) asymptotic properties of the estimator under regularity conditions. We also characterize the effect of sharing one block of simulation draws across markets and show how to correct for it. Monte Carlo simulations show that the SNFP estimator achieves statistical accuracy comparable to the NFP estimator, and in our benchmark a single online pass estimates a model with 100 million markets in about 5.5 hours. An empirical application using scanner data further demonstrates the practical advantages of SNFP for large-scale demand estimation.

The temporal dependence inherent in time series models poses a fundamental challenge for distributed Bayesian inference, as it precludes embarrassingly parallel algorithms based on naive independence assumptions in the time domain. We propose a frequency-domain framework for scalable Bayesian inference in stationary time series that exploits the asymptotic independence underlying the Whittle likelihood. To exploit parallel computing resources, we develop a distributed fast Fourier transform and integrate it with embarrassingly parallel Markov chain Monte Carlo (MCMC) algorithms within a modern cluster-computing framework. This enables the analysis of time series that exceed the memory capacity of a single computational node or for which computation time is a bottleneck. The proposed methodology is compatible with a broad class of existing divide-and-conquer algorithms for independent data by applying them to frequency-domain rather than time-domain partitions. We establish that the error of the spectral divide-and-conquer MCMC posterior approximation relative to the exact time-domain posterior converges to zero in probability in a shrinking neighbourhood of the full-data Whittle posterior mode. The corresponding convergence rate is also derived. Across several experiments, we demonstrate that our approach provides accurate approximations to the full-data Whittle posterior. The proposed method is shown to outperform the current state-of-the-art time-domain divide-and-conquer methodology, particularly for highly persistent processes. The methodology is further illustrated by fitting a semi-long range model to a long meteorological time series.