Search

Numerical modeling errors are unavoidable in finite element analysis. The presence of model errors inherently reflects both model accuracy and uncertainty. To date there have been few methods for explicitly quantifying errors at points of interest (e.g. at finite element nodes). The lack of explicit model error approximators has been addressed recently with the emergence of machine learning (ML), which closes the loop between numerical model features/solutions and explicit model error approximations. In this paper, we propose physics-informed neural networks (PINNs) for simultaneous numerical model error approximation and superresolution. To test our approach, numerical data was generated using finite element simulations on a two-dimensional elastic plate with a central opening. Four- and eight-node quadrilateral elements were used in the discretization to represent the reduced-order and higher-order models, respectively. It was found that the developed PINNs effectively predict model errors in both x and y displacement fields with small differences between predictions and ground truth. Our findings demonstrate that the integration of physics-informed loss functions enables neural networks (NNs) to surpass a purely data-driven approach for approximating model errors.

This paper describes the R package fdesigns that implements a methodology for identifying Bayesian optimal experimental designs for models whose factor settings are functions, known as profile factors. This type of experiments which involve factors that vary dynamically over time, presenting unique challenges in both estimation and design due to the infinite-dimensional nature of functions. The package fdesigns implements a dimension reduction method leveraging basis functions of the B-spline basis system. The package fdesigns contains functions that effectively reduce the design problem to the optimisation of basis coefficients for functional linear functional generalised linear models, and it accommodates various options. Applications of the fdesigns package are demonstrated through a series of examples that showcase its capabilities in identifying optimal designs for functional linear and generalised linear models. The examples highlight how the package's functions can be used to efficiently design experiments involving both profile and scalar factors, including interactions and polynomial effects.

Global variable importance measures are commonly used to interpret machine learning model results. Local variable importance techniques assess how variables contribute to individual observations rather than the entire dataset. Current methods typically fail to accurately reflect locally dependent relationships between variables and instead focus on marginal importance values. Additionally, they are not natively adapted for multi-class classification problems. We propose a new model-agnostic method for calculating local variable importance, CLIQUE, that captures locally dependent relationships, contains improvements over permutation-based methods, and can be directly applied to multi-class classification problems. Simulated and real-world examples show that CLIQUE emphasizes locally dependent information and properly reduces bias in regions where variables do not affect the response.

Computing expected information gain (EIG) from prior to posterior (equivalently, mutual information between candidate observations and model parameters or other quantities of interest) is a fundamental challenge in Bayesian optimal experimental design. We formulate flexible transport-based schemes for EIG estimation in general nonlinear/non-Gaussian settings, compatible with both standard and implicit Bayesian models. These schemes are representative of two-stage methods for estimating or bounding EIG using marginal and conditional density estimates. In this setting, we analyze the optimal allocation of samples between training (density estimation) and approximation of the outer prior expectation. We show that with this optimal sample allocation, the MSE of the resulting EIG estimator converges more quickly than that of a standard nested Monte Carlo scheme. We then address the estimation of EIG in high dimensions, by deriving gradient-based upper bounds on the mutual information lost by projecting the parameters and/or observations to lower-dimensional subspaces. Minimizing these upper bounds yields projectors and hence low-dimensional EIG approximations that outperform approximations obtained via other linear dimension reduction schemes. Numerical experiments on a PDE-constrained Bayesian inverse problem also illustrate a favorable trade-off between dimension truncation and the modeling of non-Gaussianity, when estimating EIG from finite samples in high dimensions.

Stability selection is a widely adopted resampling-based framework for high-dimensional structure estimation and variable selection. However, the concept of 'stability' is often narrowly addressed, primarily through examining selection frequencies, or 'stability paths'. This paper seeks to broaden the use of an established stability estimator to evaluate the overall stability of the stability selection framework, moving beyond single-variable analysis. We suggest that the stability estimator offers two advantages: it can serve as a reference to reflect the robustness of the outcomes obtained and help identify an optimal regularization value to improve stability. By determining this value, we aim to calibrate key stability selection parameters, namely, the decision threshold and the expected number of falsely selected variables, within established theoretical bounds. Furthermore, we explore a novel selection criterion based on this regularization value. With the asymptotic distribution of the stability estimator previously established, convergence to true stability is ensured, allowing us to observe stability trends over successive sub-samples. This approach sheds light on the required number of sub-samples addressing a notable gap in prior studies. The 'stabplot' package is developed to facilitate the use of the plots featured in this manuscript, supporting their integration into further statistical analysis and research workflows.

Existing guarantees for algorithms sampling from nonlogconcave measures on $\mathbb{R}^d$ are generally inexplicit or unscalable. Even for the class of measures with logdensities that have bounded Hessians and are strongly concave outside a Euclidean ball of radius $R$, no available theory is comprehensively satisfactory with respect to both $R$ and $d$. In this paper, it is shown that complete polynomial complexity can in fact be achieved if $R\leq c\sqrt{d}$, whilst an exponential number of point evaluations is generally necessary for any algorithm as soon as $R\geq C\sqrt{d}$ for constants $C>c>0$. A simple importance sampler with tail-matching proposal achieves the former, owing to a blessing of dimensionality. On the other hand, if strong concavity outside a ball is replaced by a distant dissipativity condition, then sampling guarantees must generally scale exponentially with $d$ in all parameter regimes.

Understanding complex interactions within microbiomes is essential for exploring their roles in health and disease. However, constructing reliable microbiome networks often poses a challenge due to variations in the output of different network inference algorithms. To address this issue, we present CMiNet, an R package designed to generate a consensus microbiome network by integrating results from multiple established network construction methods. CMiNet incorporates nine widely used algorithms, including Pearson, Spearman, Biweight Midcorrelation (Bicor), SparCC, SpiecEasi, SPRING, GCoDA, and CCLasso, along with a novel algorithm based on conditional mutual information (CMIMN). By combining the strengths of these algorithms, CMiNet generates a single, weighted consensus network that provides a more stable and comprehensive representation of microbial interactions. The package includes customizable functions for network construction, visualization, and analysis, allowing users to explore network structures at different threshold levels and assess connectivity and reliability. CMiNet is designed to handle both quantitative and compositional data, ensuring broad applicability for researchers aiming to understand the intricate relationships within microbiome communities. Availability: Source code is freely available at https://github.com/solislemuslab/CMiNet.

This paper aims at improving the convergence to equilibrium of finite ergodic Markov chains via permutations and projections. First, we prove that a specific mixture of permuted Markov chains arises naturally as a projection under the KL divergence or the squared-Frobenius norm. We then compare various mixing properties of the mixture with other competing Markov chain samplers and demonstrate that it enjoys improved convergence. This geometric perspective motivates us to propose samplers based on alternating projections to combine different permutations and to analyze their rate of convergence. We give necessary, and under some additional assumptions also sufficient, conditions for the projection to achieve stationarity in the limit in terms of the trace of the transition matrix. We proceed to discuss tuning strategies of the projection samplers when these permutations are viewed as parameters. Along the way, we reveal connections between the mixture and a Markov chain Sylvester's equation as well as assignment problems, and highlight how these can be used to understand and improve Markov chain mixing. We provide two examples as illustrations. In the first example, the projection sampler (with a suitable choice of the permutation) improves upon Metropolis-Hastings in a discrete bimodal distribution with a reduced relaxation time from exponential to polynomial in the system size, while in the second example, the mixture of permuted Markov chain yields a mixing time that is logarithmic in system size (with high probability under random permutation), compared to a linear mixing time in the Diaconis-Holmes-Neal sampler., Comment: 42 pages, 2 figures

We present the R package MSTest, which implements hypothesis testing procedures to identify the number of regimes in Markov switching models. These models have wide-ranging applications in economics, finance, and numerous other fields. The MSTest package includes the Monte Carlo likelihood ratio test procedures proposed by Rodriguez-Rondon and Dufour (2024), the moment-based tests of Dufour and Luger (2017), the parameter stability tests of Carrasco, Hu, and Ploberger (2014), and the likelihood ratio test of Hansen (1992). Additionally, the package enables users to simulate and estimate univariate and multivariate Markov switching and hidden Markov processes, using the expectation-maximization (EM) algorithm or maximum likelihood estimation (MLE). We demonstrate the functionality of the MSTest package through both simulation experiments and an application to U.S. GNP growth data.

Sequential sampling models (SSMs) are a widely used framework describing decision-making as a stochastic, dynamic process of evidence accumulation. SSMs popularity across cognitive science has driven the development of various software packages that lower the barrier for simulating, estimating, and comparing existing SSMs. Here, we present a software tool, SequentialSamplingModels.jl (SSM.jl), designed to make SSM simulations more accessible to Julia users, and to integrate with the Julia ecosystem. We demonstrate the basic use of SSM.jl for simulation, plotting, and Bayesian inference.

Parallel-in-time (PinT) techniques have been proposed to solve systems of time-dependent differential equations by parallelizing the temporal domain. Among them, Parareal computes the solution sequentially using an inaccurate (fast) solver, and then "corrects" it using an accurate (slow) integrator that runs in parallel across temporal subintervals. This work introduces RandNet-Parareal, a novel method to learn the discrepancy between the coarse and fine solutions using random neural networks (RandNets). RandNet-Parareal achieves speed gains up to x125 and x22 compared to the fine solver run serially and Parareal, respectively. Beyond theoretical guarantees of RandNets as universal approximators, these models are quick to train, allowing the PinT solution of partial differential equations on a spatial mesh of up to $10^5$ points with minimal overhead, dramatically increasing the scalability of existing PinT approaches. RandNet-Parareal's numerical performance is illustrated on systems of real-world significance, such as the viscous Burgers' equation, the Diffusion-Reaction equation, the two- and three-dimensional Brusselator, and the shallow water equation., Comment: Accepted at the 38th Conference on Neural Information Processing Systems (NeurIPS 2024)

A statistical framework we call CQUESST (Carbon Quantification and Uncertainty from Evolutionary Soil STochastics), which models carbon sequestration and cycling in soils, is applied to a long-running agricultural experiment that controls for crop type, tillage, and season. The experiment, known as the Millenium Tillage Trial (MTT), ran on 42 field-plots for ten years from 2000-2010; here CQUESST is used to model soil carbon dynamically in six pools, in each of the 42 agricultural plots, and on a monthly time step for a decade. We show how CQUESST can be used to estimate soil-carbon cycling rates under different treatments. Our methods provide much-needed statistical tools for quantitatively inferring the effectiveness of different experimental treatments on soil-carbon sequestration. The decade-long data are of multiple observation types, and these interacting time series are ingested into a fully Bayesian model that has a dynamic stochastic model of multiple pools of soil carbon at its core. CQUESST's stochastic model is motivated by the deterministic RothC soil-carbon model based on nonlinear difference equations. We demonstrate how CQUESST can estimate soil-carbon fluxes for different experimental treatments while acknowledging uncertainties in soil-carbon dynamics, in physical parameters, and in observations. CQUESST is implemented efficiently in the probabilistic programming language Stan using its MapReduce parallelization, and it scales well for large numbers of field-plots, using software libraries that allow for computation to be shared over multiple nodes of high-performance computing clusters.

Many workers in cities across the world, who have been teleworking because of the COVID-19 pandemic, are expected to be back to their commutes. As this process is believed to be gradual and telecommuting is likely to remain an option for many workers, hybrid model and flexible schedules might become the norm in the future. This variable work schedules allows employees to commute outside of traditional rush hours. Moreover, many studies showed that commuters might be skeptical of using trains, buses, and carpools and could turn to personal vehicles to get to work, which might increase congestion and crashes in the roads. This study attempts to provide information on the safest time to commute to Washington, DC area analyzing historical traffic crash data before the COVID-19 pandemic. It also aims to advance our understanding of traffic crashes and other relating factors such as weather in the Washington, DC area. We created a model to predict crashes by time of the day, using a negative binomial regression after rejecting a Poisson regression, and additionally explored the validity of a Random Forest regression. Our main consideration for an eventual application of this study is to reduce crashes in Washington DC, using this tool that provides people with better options on when to commute and when to telework, if available. The study also provides policymakers and researchers with real-world insights that decrease the number of traffic crashes to help achieve the goals of The Vision Zero Initiative adopted by the district.

Monte Carlo methods, such as Markov chain Monte Carlo (MCMC), remain the most regularly-used approach for implementing Bayesian inference. However, the computational cost of these approaches usually scales worse than linearly with the dimension of the parameter space, limiting their use for models with a large number of parameters. However, it is not uncommon for models to have identifiability problems. In this case, the likelihood is not informative about some subspaces of the parameter space, and hence the model effectively has a dimension that is lower than the number of parameters. Constantine et al. (2016) and Schuster et al. (2017) introduced the concept of directly using a Metropolis-Hastings (MH) MCMC on the subspaces of the parameter space that are informed by the likelihood as a means to reduce the dimension of the parameter space that needs to be explored. The same paper introduces an approach for identifying such a subspace in the case where it is a linear transformation of the parameter space: this subspace is known as the active subspace. This paper introduces sequential Monte Carlo (SMC) methods that make use of an active subspace. As well as an SMC counterpart to MH approach of Schuster et al. (2017), we introduce an approach to learn the active subspace adaptively, and an SMC$^{2}$ approach that is more robust to the linearity assumptions made when using active subspaces.

Gaussian process are a widely-used statistical tool for conducting non-parametric inference in applied sciences, with many computational packages available to fit to data and predict future observations. We study the use of the Greta software for Bayesian inference to apply Gaussian process regression to spatio-temporal data of infectious disease outbreaks and predict future disease spread. Greta builds on Tensorflow, making it comparatively easy to take advantage of the significant gain in speed offered by GPUs. In these complex spatio-temporal models, we show a reduction of up to 70\% in computational time relative to fitting the same models on CPUs. We show how the choice of covariance kernel impacts the ability to infer spread and extrapolate to unobserved spatial and temporal units. The inference pipeline is applied to weekly incidence data on tuberculosis in the East and West Midlands regions of England over a period of two years.

The Gaussian process (GP) is a widely used probabilistic machine learning method for stochastic function approximation, stochastic modeling, and analyzing real-world measurements of nonlinear processes. Unlike many other machine learning methods, GPs include an implicit characterization of uncertainty, making them extremely useful across many areas of science, technology, and engineering. Traditional implementations of GPs involve stationary kernels (also termed covariance functions) that limit their flexibility and exact methods for inference that prevent application to data sets with more than about ten thousand points. Modern approaches to address stationarity assumptions generally fail to accommodate large data sets, while all attempts to address scalability focus on approximating the Gaussian likelihood, which can involve subjectivity and lead to inaccuracies. In this work, we explicitly derive an alternative kernel that can discover and encode both sparsity and nonstationarity. We embed the kernel within a fully Bayesian GP model and leverage high-performance computing resources to enable the analysis of massive data sets. We demonstrate the favorable performance of our novel kernel relative to existing exact and approximate GP methods across a variety of synthetic data examples. Furthermore, we conduct space-time prediction based on more than one million measurements of daily maximum temperature and verify that our results outperform state-of-the-art methods in the Earth sciences. More broadly, having access to exact GPs that use ultra-scalable, sparsity-discovering, nonstationary kernels allows GP methods to truly compete with a wide variety of machine learning methods.

Stochastic simulation models are generative models that mimic complex systems to help with decision-making. The reliability of these models heavily depends on well-calibrated input model parameters. However, in many practical scenarios, only output-level data are available to learn the input model parameters, which is challenging due to the often intractable likelihood of the stochastic simulation model. Moreover, stochastic simulation models are frequently inexact, with discrepancies between the model and the target system. No existing methods can effectively learn and quantify the uncertainties of input parameters using only output-level data. In this paper, we propose to learn differentiable input parameters of stochastic simulation models using output-level data via kernel score minimization with stochastic gradient descent. We quantify the uncertainties of the learned input parameters using a frequentist confidence set procedure based on a new asymptotic normality result that accounts for model inexactness. The proposed method is evaluated on exact and inexact G/G/1 queueing models., Comment: 31 pages, 12 tables, 4 figures

In aims to uncover insights into medical decision-making embedded within observational data from clinical settings, we present a novel application of Inverse Reinforcement Learning (IRL) that identifies suboptimal clinician actions based on the actions of their peers. This approach centers two stages of IRL with an intermediate step to prune trajectories displaying behavior that deviates significantly from the consensus. This enables us to effectively identify clinical priorities and values from ICU data containing both optimal and suboptimal clinician decisions. We observe that the benefits of removing suboptimal actions vary by disease and differentially impact certain demographic groups., Comment: 13 pages, 4 figures

Generalized linear mixed models (GLMMs) are a widely used tool in statistical analysis. The main bottleneck of many computational approaches lies in the inversion of the high dimensional precision matrices associated with the random effects. Such matrices are typically sparse; however, the sparsity pattern resembles a multi partite random graph, which does not lend itself well to default sparse linear algebra techniques. Notably, we show that, for typical GLMMs, the Cholesky factor is dense even when the original precision is sparse. We thus turn to approximate iterative techniques, in particular to the conjugate gradient (CG) method. We combine a detailed analysis of the spectrum of said precision matrices with results from random graph theory to show that CG-based methods applied to high-dimensional GLMMs typically achieve a fixed approximation error with a total cost that scales linearly with the number of parameters and observations. Numerical illustrations with both real and simulated data confirm the theoretical findings, while at the same time illustrating situations, such as nested structures, where CG-based methods struggle.

Communities on the web rely on open conversation forums for a number of tasks, including governance, information sharing, and decision making. However these forms of collective deliberation can often result in biased outcomes. A prime example are Articles for Deletion (AfD) discussions on Wikipedia, which allow editors to gauge the notability of existing articles, and that, as prior work has suggested, may play a role in perpetuating the notorious gender gap of Wikipedia. Prior attempts to address this question have been hampered by access to narrow observation windows, reliance on limited subsets of both biographies and editorial outcomes, and by potential confounding factors. To address these limitations, here we adopt a competing risk survival framework to fully situate biographical AfD discussions within the full editorial cycle of Wikipedia content. We find that biographies of women are nominated for deletion faster than those of men, despite editors taking longer to reach a consensus for deletion of women, even after controlling for the size of the discussion. Furthermore, we find that AfDs about historical figures show a strong tendency to result into the redirecting or merging of the biography under discussion into other encyclopedic entries, and that there is a striking gender asymmetry: biographies of women are redirected or merged into biographies of men more often than the other way round. Our study provides a more complete picture of the role of AfD in the gender gap of Wikipedia, with implications for the governance of the open knowledge infrastructure of the web.

Today, cheap numerical hardware offers huge amounts of parallel computing power, much of which is used for the task of fitting neural networks to data. Adoption of this hardware to accelerate statistical Markov chain Monte Carlo (MCMC) applications has been much slower. In this chapter, we suggest some patterns for speeding up MCMC workloads using the hardware (e.g., GPUs, TPUs) and software (e.g., PyTorch, JAX) that have driven progress in deep learning over the last fifteen years or so. We offer some intuitions for why these new systems are so well suited to MCMC, and show some examples (with code) where we use them to achieve dramatic speedups over a CPU-based workflow. Finally, we discuss some potential pitfalls to watch out for.

This study presents a scalable Bayesian estimation algorithm for sparse estimation in exploratory item factor analysis based on a classical Bayesian estimation method, namely Bayesian joint modal estimation (BJME). BJME estimates the model parameters and factor scores that maximize the complete-data joint posterior density. Simulation studies show that the proposed algorithm has high computational efficiency and accuracy in variable selection over latent factors and the recovery of the model parameters. Moreover, we conducted a real data analysis using large-scale data from a psychological assessment that targeted the Big Five personality traits. This result indicates that the proposed algorithm achieves computationally efficient parameter estimation and extracts the interpretable factor loading structure.

The Minimum Volume Covering Ellipsoid (MVCE) problem, characterised by $n$ observations in $d$ dimensions where $n \gg d$, can be computationally very expensive in the big data regime. We apply methods from randomised numerical linear algebra to develop a data-driven leverage score sampling algorithm for solving MVCE, and establish theoretical error bounds and a convergence guarantee. Assuming the leverage scores follow a power law decay, we show that the computational complexity of computing the approximation for MVCE is reduced from $\mathcal{O}(nd^2)$ to $\mathcal{O}(nd + \text{poly}(d))$, which is a significant improvement in big data problems. Numerical experiments demonstrate the efficacy of our new algorithm, showing that it substantially reduces computation time and yields near-optimal solutions., Comment: 20 pages, 3 figures

Cluster-weighted factor analyzers (CWFA) are a versatile class of mixture models designed to estimate the joint distribution of a random vector that includes a response variable along with a set of explanatory variables. They are particularly valuable in situations involving high dimensionality. This paper enhances CWFA models in two notable ways. First, it enables the prediction of multiple response variables while considering their potential interactions. Second, it identifies factors associated with disjoint groups of explanatory variables, thereby improving interpretability. This development leads to the introduction of the multivariate cluster-weighted disjoint factor analyzers (MCWDFA) model. An alternating expectation-conditional maximization algorithm is employed for parameter estimation. The effectiveness of the proposed model is assessed through an extensive simulation study that examines various scenarios. The proposal is applied to crime data from the United States, sourced from the UCI Machine Learning Repository, with the aim of capturing potential latent heterogeneity within communities and identifying groups of socio-economic features that are similarly associated with factors predicting crime rates. Results provide valuable insights into the underlying structures influencing crime rates which may potentially be helpful for effective cluster-specific policymaking and social interventions., Comment: 6 Tables, 7 Figures, 49 pages

The use of mediation analysis has become increasingly popular in various research fields in recent years. The primary objective of mediation analysis is to examine the indirect effects along the pathways from exposure to outcome. Meanwhile, the advent of data collection technology has sparked a surge of interest in the burgeoning field of big data analysis, where mediation analysis of streaming data sets has recently garnered significant attention. The enormity of the data, however, results in an augmented computational burden. The present study proposes an online updating approach to address this issue, aiming to estimate and test mediation effects in the context of linear and logistic mediation models with massive data streams. The proposed algorithm significantly enhances the computational efficiency of Sobel test, adjusted Sobel test, joint significance test, and adjusted joint significance test. We conduct a substantial number of numerical simulations to evaluate the performance of the renewable method. Two real-world examples are employed to showcase the practical applicability of this approach.

In this paper, we present a theoretical and computational workflow for the non-parametric Bayesian inference of drift and diffusion functions of autonomous diffusion processes. We base the inference on the partial differential equations arising from the infinitesimal generator of the underlying process. Following a problem formulation in the infinite-dimensional setting, we discuss optimization- and sampling-based solution methods. As preliminary results, we showcase the inference of a single-scale, as well as a multiscale process from trajectory data.

This article introduces a new symbolic regression algorithm based on the SPINEX (Similarity-based Predictions with Explainable Neighbors Exploration) family. This new algorithm (SPINEX_SymbolicRegression) adopts a similarity-based approach to identifying high-merit expressions that satisfy accuracy- and structural similarity metrics. We conducted extensive benchmarking tests comparing SPINEX_SymbolicRegression to over 180 mathematical benchmarking functions from international problem sets that span randomly generated expressions and those based on real physical phenomena. Then, we evaluated the performance of the proposed algorithm in terms of accuracy, expression similarity in terms of presence operators and variables (as compared to the actual expressions), population size, and number of generations at convergence. The results indicate that SPINEX_SymbolicRegression consistently performs well and can, in some instances, outperform leading algorithms. In addition, the algorithm's explainability capabilities are highlighted through in-depth experiments.

General state-space models (SSMs) are widely used in statistical machine learning and are among the most classical generative models for sequential time-series data. SSMs, comprising latent Markovian states, can be subjected to variational inference (VI), but standard VI methods like the importance-weighted autoencoder (IWAE) lack functionality for streaming data. To enable online VI in SSMs when the observations are received in real time, we propose maximising an IWAE-type variational lower bound on the asymptotic contrast function, rather than the standard IWAE ELBO, using stochastic approximation. Unlike the recursive maximum likelihood method, which directly maximises the asymptotic contrast, our approach, called online sequential IWAE (OSIWAE), allows for online learning of both model parameters and a Markovian recognition model for inferring latent states. By approximating filter state posteriors and their derivatives using sequential Monte Carlo (SMC) methods, we create a particle-based framework for online VI in SSMs. This approach is more theoretically well-founded than recently proposed online variational SMC methods. We provide rigorous theoretical results on the learning objective and a numerical study demonstrating the method's efficiency in learning model parameters and particle proposal kernels.

Robust estimation provides essential tools for analyzing data that contain outliers, ensuring that statistical models remain reliable even in the presence of some anomalous data. While robust methods have long been available in R, users of Python have lacked a comprehensive package that offers these methods in a cohesive framework. RobPy addresses this gap by offering a wide range of robust methods in Python, built upon established libraries including NumPy, SciPy, and scikit-learn. This package includes tools for robust preprocessing, univariate estimation, covariance matrices, regression, and principal component analysis, which are able to detect outliers and to mitigate their effect. In addition, RobPy provides specialized diagnostic plots for visualizing casewise and cellwise outliers. This paper presents the structure of the RobPy package, demonstrates its functionality through examples, and compares its features to existing implementations in other statistical software. By bringing robust methods to Python, RobPy enables more users to perform robust data analysis in a modern and versatile programming language.

Rahimi and Recht [31] introduced the idea of decomposing shift-invariant kernels by randomly sampling from their spectral distribution. This famous technique, known as Random Fourier Features (RFF), is in principle applicable to any shift-invariant kernel whose spectral distribution can be identified and simulated. In practice, however, it is usually applied to the Gaussian kernel because of its simplicity, since its spectral distribution is also Gaussian. Clearly, simple spectral sampling formulas would be desirable for broader classes of kernel functions. In this paper, we propose to decompose spectral kernel distributions as a scale mixture of $\alpha$-stable random vectors. This provides a simple and ready-to-use spectral sampling formula for a very large class of multivariate shift-invariant kernels, including exponential power kernels, generalized Mat\'ern kernels, generalized Cauchy kernels, as well as newly introduced kernels such as the Beta, Kummer, and Tricomi kernels. In particular, we show that the spectral densities of all these kernels are scale mixtures of the multivariate Gaussian distribution. This provides a very simple way to modify existing Random Fourier Features software based on Gaussian kernels to cover a much richer class of multivariate kernels. This result has broad applications for support vector machines, kernel ridge regression, Gaussian processes, and other kernel-based machine learning techniques for which the random Fourier features technique is applicable., Comment: 16 pages, 16 figures

Most methods in explainable AI (XAI) focus on providing reasons for the prediction of a given set of features. However, we solve an inverse explanation problem, i.e., given the deviation of a label, find the reasons of this deviation. We use a Bayesian framework to recover the ``true'' features, conditioned on the observed label value. We efficiently explain the deviation of a label value from the mode, by identifying and ranking the influential features using the ``distances'' in the ANOVA functional decomposition. We show that the new method is more human-intuitive and robust than methods based on mean values, e.g., SHapley Additive exPlanations (SHAP values). The extra costs of solving a Bayesian inverse problem are dimension-independent.

Efficient sampling from un-normalized target distributions is pivotal in scientific computing and machine learning. While neural samplers have demonstrated potential with a special emphasis on sampling efficiency, existing neural implicit samplers still have issues such as poor mode covering behavior, unstable training dynamics, and sub-optimal performances. To tackle these issues, in this paper, we introduce Denoising Fisher Training (DFT), a novel training approach for neural implicit samplers with theoretical guarantees. We frame the training problem as an objective of minimizing the Fisher divergence by deriving a tractable yet equivalent loss function, which marks a unique theoretical contribution to assessing the intractable Fisher divergences. DFT is empirically validated across diverse sampling benchmarks, including two-dimensional synthetic distribution, Bayesian logistic regression, and high-dimensional energy-based models (EBMs). Notably, in experiments with high-dimensional EBMs, our best one-step DFT neural sampler achieves results on par with MCMC methods with up to 200 sampling steps, leading to a substantially greater efficiency over 100 times higher. This result not only demonstrates the superior performance of DFT in handling complex high-dimensional sampling but also sheds light on efficient sampling methodologies across broader applications.

Algorithms for generating random numbers that follow a gamma distribution with shape parameter less than unity are proposed. Acceptance-rejection algorithms are developed, based on the generalized exponential distribution. The squeeze technique is applied to our method, and then piecewise envelope functions are further considered. The proposed methods are excellent in acceptance efficiency and promising in speed., Comment: 13 pages, 2 figures; To appear in Economics Bulletin

To understand collective network behavior in the complex human brain, pairwise correlation networks alone are insufficient for capturing the high-order interactions that extend beyond pairwise interactions and play a crucial role in brain network dynamics. These interactions often reveal intricate relationships among multiple brain networks, significantly influencing cognitive processes. In this study, we explored the correlation of correlation networks and topological network analysis with resting-state fMRI to gain deeper insights into these higher-order interactions and their impact on the topology of brain networks, ultimately enhancing our understanding of brain function. We observed that the correlation of correlation networks highlighted network connections while preserving the topological structure of correlation networks. Our findings suggested that the correlation of correlation networks surpassed traditional correlation networks, showcasing considerable potential for applications in various areas of network science. Moreover, after applying topological network analysis to the correlation of correlation networks, we observed that some high-order interaction hubs predominantly occurred in primary and high-level cognitive areas, such as the visual and fronto-parietal regions. These high-order hubs played a crucial role in information integration within the human brain., Comment: 4 pages, 2 figures, 1 table; Submitted to IEEE International Symposium on Biomedical Imaging (ISBI 2025)

The human brain dynamically integrated and configured information to adapt to the environment. To capture these changes over time, dynamic second-order functional connectivity was typically used to capture transient brain patterns. However, dynamic second-order functional connectivity typically ignored interactions beyond pairwise relationships. To address this limitation, we utilized dynamic triple interactions to investigate multiscale network interactions in the brain. In this study, we evaluated a resting-state fMRI dataset that included individuals with psychotic disorders (PD). We first estimated dynamic triple interactions using resting-state fMRI. After clustering, we estimated cohort-specific and cohort-common states for controls (CN), schizophrenia (SZ), and schizoaffective disorder (SAD). From the cohort-specific states, we observed significant triple interactions, particularly among visual, subcortical, and somatomotor networks, as well as temporal and higher cognitive networks in SZ. In SAD, key interactions involved temporal networks in the initial state and somatomotor networks in subsequent states. From the cohort-common states, we observed that high-cognitive networks were primarily involved in SZ and SAD compared to CN. Furthermore, the most significant differences between SZ and SAD also existed in high-cognitive networks. In summary, we studied PD using dynamic triple interaction, the first time such an approach has been used to study PD. Our findings highlighted the significant potential of dynamic high-order functional connectivity, paving the way for new avenues in the study of the healthy and disordered human brain., Comment: 4 pages, 3 figures; Submitted to IEEE International Symposium on Biomedical Imaging (ISBI 2025)

Gibbs sampling is a common procedure used to fit finite mixture models. However, it is known to be slow to converge when exploring correlated regions of a parameter space and so blocking correlated parameters is sometimes implemented in practice. This is straightforward to visualize in contexts like low-dimensional multivariate Gaussian distributions, but more difficult for mixture models because of the way latent variable assignment and cluster-specific parameters influence one another. Here we analyze correlation in the space of latent variables and show that latent variables of outlier observations equidistant between component distributions can exhibit significant correlation that is not bounded away from one, suggesting they can converge very slowly to their stationary distribution. We provide bounds on convergence rates to a modification of the stationary distribution and propose a blocked sampling procedure that significantly reduces autocorrelation in the latent variable Markov chain, which we demonstrate in simulation.

Bayesian Kernel Machine Regression (BKMR) has emerged as a powerful tool to detect negative health effects from exposure to complex multi-pollutant mixtures. However, its performance is degraded when data deviate from normality. In this comprehensive simulation analysis, we show that BKMR's power and test size vary under different distributions and covariance matrix structures. Our results demonstrate specifically that BKMR's robustness is influenced by the response's coefficient of variation (CV), resulting in reduced accuracy to detect true effects when data are skewed. Test sizes become uncontrolled (> 0.05) as CV values increase, leading to inflated false detection rates. However, we find that BKMR effectively utilizes off-diagonal covariance information corresponding to predictor interdependencies, increasing statistical power and accuracy. To achieve reliable and accurate results, we advocate for scrutiny of data skewness and covariance before applying BKMR, particularly when used to predict cognitive decline from blood/urine heavy metal concentrations in environmental health contexts., Comment: This article is submitted to the Journal of Statistical Computation and Simulation and currently under review

Nudging is a popular algorithmic strategy in numerical filtering to deal with the problem of inference in high-dimensional dynamical systems. We demonstrate in this paper that general nudging techniques can also tackle another crucial statistical problem in filtering, namely the misspecification of the transition model. Specifically, we rely on the formulation of nudging as a general operation increasing the likelihood and prove analytically that, when applied carefully, nudging techniques implicitly define state-space models (SSMs) that have higher marginal likelihoods for a given (fixed) sequence of observations. This provides a theoretical justification of nudging techniques as data-informed algorithmic modifications of SSMs to obtain robust models under misspecified dynamics. To demonstrate the use of nudging, we provide numerical experiments on linear Gaussian SSMs and a stochastic Lorenz 63 model with misspecified dynamics and show that nudging offers a robust filtering strategy for these cases.

We develop EigenVI, an eigenvalue-based approach for black-box variational inference (BBVI). EigenVI constructs its variational approximations from orthogonal function expansions. For distributions over $\mathbb{R}^D$, the lowest order term in these expansions provides a Gaussian variational approximation, while higher-order terms provide a systematic way to model non-Gaussianity. These approximations are flexible enough to model complex distributions (multimodal, asymmetric), but they are simple enough that one can calculate their low-order moments and draw samples from them. EigenVI can also model other types of random variables (e.g., nonnegative, bounded) by constructing variational approximations from different families of orthogonal functions. Within these families, EigenVI computes the variational approximation that best matches the score function of the target distribution by minimizing a stochastic estimate of the Fisher divergence. Notably, this optimization reduces to solving a minimum eigenvalue problem, so that EigenVI effectively sidesteps the iterative gradient-based optimizations that are required for many other BBVI algorithms. (Gradient-based methods can be sensitive to learning rates, termination criteria, and other tunable hyperparameters.) We use EigenVI to approximate a variety of target distributions, including a benchmark suite of Bayesian models from posteriordb. On these distributions, we find that EigenVI is more accurate than existing methods for Gaussian BBVI., Comment: 25 pages, 9 figures. Advances in Neural Information Processing Systems (NeurIPS), 2024

We introduce a novel method for obtaining a wide variety of moments of a random variable with a well-defined moment-generating function (MGF). We derive new expressions for fractional moments and fractional absolute moments, both central and non-central moments. The new moment expressions are relatively simple integrals that involve the MGF, but do not require its derivatives. We label the new method CMGF because it uses a complex extension of the MGF and can be used to obtain complex moments. We illustrate the new method with three applications where the MGF is available in closed-form, while the corresponding densities and the derivatives of the MGF are either unavailable or very difficult to obtain.

We study multiproposal Markov chain Monte Carlo algorithms, such as Multiple-try or generalised Metropolis-Hastings schemes, which have recently received renewed attention due to their amenability to parallel computing. First, we prove that no multiproposal scheme can speed-up convergence relative to the corresponding single proposal scheme by more than a factor of $K$, where $K$ denotes the number of proposals at each iteration. This result applies to arbitrary target distributions and it implies that serial multiproposal implementations are always less efficient than single proposal ones. Secondly, we consider log-concave distributions over Euclidean spaces, proving that, in this case, the speed-up is at most logarithmic in $K$, which implies that even parallel multiproposal implementations are fundamentally limited in the computational gain they can offer. Crucially, our results apply to arbitrary multiproposal schemes and purely rely on the two-step structure of the associated kernels (i.e. first generate $K$ candidate points, then select one among those). Our theoretical findings are validated through numerical simulations., Comment: 19 pages, 1 figure

Computational capability often falls short when confronted with massive data, posing a common challenge in establishing a statistical model or statistical inference method dealing with big data. While subsampling techniques have been extensively developed to downsize the data volume, there is a notable gap in addressing the unique challenge of handling extensive reliability data, in which a common situation is that a large proportion of data is censored. In this article, we propose an efficient subsampling method for reliability analysis in the presence of censoring data, intending to estimate the parameters of lifetime distribution. Moreover, a novel subsampling method for subsampling from severely censored data is proposed, i.e., only a tiny proportion of data is complete. The subsampling-based estimators are given, and their asymptotic properties are derived. The optimal subsampling probabilities are derived through the L-optimality criterion, which minimizes the trace of the product of the asymptotic covariance matrix and a constant matrix. Efficient algorithms are proposed to implement the proposed subsampling methods to address the challenge that optimal subsampling strategy depends on unknown parameter estimation from full data. Real-world hard drive dataset case and simulative empirical studies are employed to demonstrate the superior performance of the proposed methods., Comment: 28 pages, 3 figures, to be published in Statistics and Its Interface

We generalize the non-negative matrix factorization algorithm of Lee and Seung to accept a weighted norm, and to support ridge and Lasso regularization. We recast the Lee and Seung multiplicative update as an additive update which does not get stuck on zero values. We apply the companion R package rnnmf to the problem of finding a reduced rank representation of a database of cocktails., Comment: 6 figures

Variance is a basic metric to evaluate the degree of data dispersion, and it is also frequently used in the realm of statistics. However, due to the computing variance and the large dataset being time-consuming, there is an urge to accelerate this computing process. The paper suggests a new method to reduce the time of this computation, it assumes a scenario in which we already know the variance of the original dataset, and the whole variance of this merge dataset could be expressed in the form of addition between the original variance and a remainder term. When we want to calculate the total variance after this adds up, the method only needs to calculate the remainder to get the result instead of recalculating the total variance again, which we named this type of method as PKA(Prior Knowledge Acceleration). The paper mathematically proves the effectiveness of PKA in variance calculation, and the conditions for this method to accelerate properly.

Stability selection is a versatile framework for structure estimation and variable selection in high-dimensional setting, primarily grounded in frequentist principles. In this paper, we propose an enhanced methodology that integrates Bayesian analysis to refine the inference of inclusion probabilities within the stability selection framework. Traditional approaches rely on selection frequencies for decision-making, often disregarding domain-specific knowledge and failing to account for the inherent uncertainty in the variable selection process. Our methodology uses prior information to derive posterior distributions of inclusion probabilities, thereby improving both inference and decision-making. We present a two-step process for engaging with domain experts, enabling statisticians to elucidate prior distributions informed by expert knowledge while allowing experts to control the weight of their input on the final results. Using posterior distributions, we offer Bayesian credible intervals to quantify uncertainty in the variable selection process. In addition, we highlight how selection frequencies can be uninformative or even misleading when covariates are correlated with each other, and demonstrate how domain expertise can alleviate such issues. Our approach preserves the versatility of stability selection and is suitable for a broad range of structure estimation challenges.

We consider an infinite server queue where the arrival and the service rates are both modulated by a stochastic environment governed by an $S$-valued stochastic process $X$ that is ergodic with a limiting measure $\pi\in \mathcal{P}(S)$. Under certain conditions when $X$ is semi-Markovian and satisfies the renewal regenerative property, long-term behavior of the total counts of people in the queue (denoted by $Y:=(Y_{t}:t\ge 0)$) becomes explicit and the limiting measure of $Y$ can be described through a well-studied affine stochastic recurrence equation (SRE) $X\stackrel{d}{=}CX+D,\,\, X\perp\!\!\!\perp (C, D)$. We propose a sampling scheme from that limiting measure with explicit convergence diagnostics. Additionally, one example is presented where the stochastic environment makes the system transient, in absence of a `no-feedback' assumption., Comment: Any comments are welcome at palmabhishek@gmail.com

This paper presents a variant of the Multinomial mixture model tailored for the unsupervised classification of short text data. Traditionally, the Multinomial probability vector in this hierarchical model is assigned a Dirichlet prior distribution. Here, however, we explore an alternative prior--the Beta-Liouville distribution--which offers a more flexible correlation structure than the Dirichlet. We examine the theoretical properties of the Beta-Liouville distribution, focusing on its conjugacy with the Multinomial likelihood. This property enables the derivation of update equations for a CAVI (Coordinate Ascent Variational Inference) variational algorithm, facilitating the approximate posterior estimation of model parameters. Additionally, we propose a stochastic variant of the CAVI algorithm that enhances scalability. The paper concludes with data examples that demonstrate effective strategies for setting the Beta-Liouville hyperparameters., Comment: 32 pages, 4 figures. Submitted

Joint models (JMs) for longitudinal and time-to-event data are an important class of biostatistical models in health and medical research. When the study population consists of heterogeneous subgroups, the standard JM may be inadequate and lead to misleading results. Joint latent class models (JLCMs) and their variants have been proposed to incorporate latent class structures into JMs. JLCMs are useful for identifying latent subgroup structures, obtaining a more nuanced understanding of the relationships between longitudinal outcomes, and improving prediction performance. We consider the generic form of JLCM, which poses significant computational challenges for both frequentist and Bayesian approaches due to the numerical intractability and multimodality of the associated model's likelihood or posterior. Focusing on the less explored Bayesian paradigm, we propose a new Bayesian inference framework to tackle key limitations in the existing method. Our algorithm leverages state-of-the-art Markov chain Monte Carlo techniques and parallel computing for parameter estimation and model selection. Through a simulation study, we demonstrate the feasibility and superiority of our proposed method over the existing approach. Our simulations also generate important computational insights and practical guidance for implementing such complex models. We illustrate our method using data from the PAQUID prospective cohort study, where we jointly investigate the association between a repeatedly measured cognitive score and the risk of dementia and the latent class structure defined from the longitudinal outcomes.

Traditional implicit generative models are capable of learning highly complex data distributions. However, their training involves distinguishing real data from synthetically generated data using adversarial discriminators, which can lead to unstable training dynamics and mode dropping issues. In this work, we build on the \textit{invariant statistical loss} (ISL) method introduced in \cite{de2024training}, and extend it to handle heavy-tailed and multivariate data distributions. The data generated by many real-world phenomena can only be properly characterised using heavy-tailed probability distributions, and traditional implicit methods struggle to effectively capture their asymptotic behavior. To address this problem, we introduce a generator trained with ISL, that uses input noise from a generalised Pareto distribution (GPD). We refer to this generative scheme as Pareto-ISL for conciseness. Our experiments demonstrate that Pareto-ISL accurately models the tails of the distributions while still effectively capturing their central characteristics. The original ISL function was conceived for 1D data sets. When the actual data is $n$-dimensional, a straightforward extension of the method was obtained by targeting the $n$ marginal distributions of the data. This approach is computationally infeasible and ineffective in high-dimensional spaces. To overcome this, we extend the 1D approach using random projections and define a new loss function suited for multivariate data, keeping problems tractable by adjusting the number of projections. We assess its performance in multidimensional generative modeling and explore its potential as a pretraining technique for generative adversarial networks (GANs) to prevent mode collapse, reporting promising results and highlighting its robustness across various hyperparameter settings.

Black-box variational inference (BBVI) scales poorly to high dimensional problems when it is used to estimate a multivariate Gaussian approximation with a full covariance matrix. In this paper, we extend the batch-and-match (BaM) framework for score-based BBVI to problems where it is prohibitively expensive to store such covariance matrices, let alone to estimate them. Unlike classical algorithms for BBVI, which use gradient descent to minimize the reverse Kullback-Leibler divergence, BaM uses more specialized updates to match the scores of the target density and its Gaussian approximation. We extend the updates for BaM by integrating them with a more compact parameterization of full covariance matrices. In particular, borrowing ideas from factor analysis, we add an extra step to each iteration of BaM -- a patch -- that projects each newly updated covariance matrix into a more efficiently parameterized family of diagonal plus low rank matrices. We evaluate this approach on a variety of synthetic target distributions and real-world problems in high-dimensional inference.