Papers

Nonlinear Mixed Effects (NLME) models are widely used in pharmacometrics and related fields to analyze hierarchical and longitudinal data. However, as the number of parameters and random effects increases, traditional methods for maximizing the marginal likelihood become computationally expensive. This paper explores the Variational Expectation Maximization (VEM) algorithm, a scalable alternative for fitting NLME models. Originally introduced in the context of probabilistic graphical models and later popularized through variational autoencoders, VEM has not been extensively applied to NLME modeling. By leveraging flexible variational families and reverse-mode automatic differentiation, VEM can efficiently maximize the marginal likelihood, scaling to NLME models with over 15,000 population parameters. This work provides a detailed description of VEM, compares it to other NLME fitting algorithms, and highlights its scalability through computational experiments. Using the Pumas statistical software, we fit two test models: 1) a standard warfarin model, and 2) an unnecessarily over-parameterized DeepNLME Friberg model with 15,410 population parameters and 16 random effects. The warfarin model was fitted to completion to demonstrate the correctness of VEM, while the DeepNLME Friberg model instead demonstrates VEM's scalability on a toy but large model. VEM improves the log likelihood steadily over hundreds of iterations at a practical per-iteration cost, while FOCE fails to complete even one iteration within a day. The model is deliberately over-parameterized for its small dataset and over-fits it, so what this experiment establishes is that VEM optimizes the objective of a model of this size at a practical cost. Applying VEM to large models that are genuinely useful is left to future work.

We consider learning from labeled data collected across multiple environments, where the data distribution may vary across these environments. This problem is commonly approached from a causal perspective, seeking invariant representations that retain causal factors while discarding spurious ones. However, this framework assumes that the environment has no direct effect on the target. In contrast, we consider settings in which this assumption fails, but still aim to learn representations that support robust prediction on average across previously unseen environments. To this end, we study representations learned by explicitly modeling variation across environments and then marginalizing that variation out. We analyze the resulting representations and characterize when they are preferable to those learned by causal invariant-representation methods. We propose a concrete method based on generalized random-intercept models, a class of predictors in which such marginalization is possible, and study their generalization properties. Empirically, we show that these models outperform invariant-learning methods across a range of challenging settings.

This paper proposes a deep reinforcement learning (DRL)-based event-triggered controller design for networked artificial pancreas (AP) systems. Although existing DRL-based AP controllers typically assume periodic control updates, networked control systems (NCSs) require a reduction in communication frequency to achieve energy-efficient operation, which is directly tied to control updates. However, jointly learning both insulin dosing and update timing significantly increases the complexity of the learning problem. To alleviate this complexity, we develop a practical DRL-based controller design that avoids explicitly learning update timing by introducing a rule-based criterion defined by changes in blood glucose. As a result, decision-making occurs at irregular intervals, and the problem is naturally formulated as a semi-Markov decision process (SMDP), for which we extend a standard DRL algorithm. Numerical experiments demonstrate that the proposed method improves communication efficiency while maintaining control performance.

Causal inference in continuous-time sequential decision problems is challenged by hidden confounding and partially observed states. We show that, under explicit structural assumptions, observability of the latent state enables identification of dynamic treatment effects through a continuous-time conditional front-door adjustment, even in the presence of hidden confounding. We derive a general adjustment formula and show that it reduces to a tractable state-space formula when unobserved contemporaneous disturbances are temporally uncorrelated. This formula expresses potential-outcome distributions under alternative treatment trajectories through the measurement model, latent dynamics, and the filtering distribution over latent states. We propose Observable Neural ODEs (ObsNODEs), Neural ODE models in observable normal form that implement this tractable adjustment for causal forecasting. ObsNODEs learn continuous-time dynamics with states reconstructible from observations, enabling outcome prediction under alternative treatment paths. Experiments on synthetic, semi-synthetic, and real-world clinical data demonstrate strong performance over recent sequence models, including external validation.

In this paper, we develop a novel inferential approach based on geometric records for estimating the tail index of heavy-tailed distributions. We construct a maximum likelihood estimator for the Pareto model and establish strong consistency and asymptotic normality, providing also an explicit expression for the asymptotic variance. These results are then extended to a broad class of Pareto-type distributions. The performance of the estimator is assessed via Monte Carlo simulation and compared with classical estimators from the literature. The proposed method is particularly well suited for settings where data arrive sequentially, as it yields smooth estimation trajectories. It is also especially advantageous in applications such as destructive testing, where measuring each item is costly. In this context, the estimator achieves a comparable level of estimation accuracy to Hill's estimator, but with a considerably lower number of fully measured items. An application to the analysis of the distribution of fluctuations of the Dow Jones Industrial Average (DJI) is also presented.

Evidence accumulation models (EAMs) provide a powerful framework for inferring latent cognitive processes from choice and reaction time data. While EAMs are traditionally limited to binary choices, recent developments have extended them to rotationally symmetric continuous responses via the circular diffusion model \citep{smith2016diffusion} and the spatially continuous diffusion model \citep{ratcliff2018decision}. Yet, such extensions are limited in scope, as many psychological constructs are measured on bounded non-rotational scales. In this paper, we bridge this gap by presenting and comparing two adaptations designed for bounded continuous data: the Half-Circular Diffusion Model (HCDM) and the Beta Drift Diffusion Model (BDDM). Because both models have intractable likelihoods, we fit them using Amortized Bayesian Inference (ABI) and compare them using Amortized Bayesian Model Comparison (ABMC). We demonstrate the complete workflow on an empirical affect dataset (N = 215), including parameter recovery, simulation-based calibration, posterior predictive checks, and model comparison. Both models accurately capture the joint distribution of responses and reaction times and yield interpretable parameters that can be reliably recovered. The model comparison further reveals a simple diagnostic for choosing between them: the dispersion of the rating distribution, with HCDM preferred for moderate spread and BDDM for highly concentrated or highly dispersed ratings. This work extends the EAM framework to a new application context, bounded continuous self-report data, and offers researchers a user-friendly toolkit for modeling the cognitive dynamics of continuous responses. We release fully documented Python code with both GPU and CPU implementations, along with example datasets.

In this paper we propose a semiparametric spatial autoregressive model that combines a linear covariate component with a nonparametrically estimated spatial term, allowing flexible dependence modeling without restrictive covariance structure while preserving interpretability. We establish asymptotic properties, including consistency and asymptotic normality, and evaluate performance through simulations and real data. Results show competitive predictive accuracy relative to geostatistical methods and improved interpretability compared to spatial econometric models.

We consider the problem of scalable sampling algorithms to fit Bayesian generalized linear mixed models on large datasets. Stochastic gradient Langevin dynamics, coupled with smooth re-parameterizations of variance parameters, produces divergent Markov chains and cannot be reliably used for sampling covariance parameters of random effects. We advocate the use of a mirror Langevin dynamics algorithm, propose the novel stochastic mirror Langevin dynamics based on data subsampling, and provide concrete guidelines for its use in a Bayesian inference framework. Based on an explicit Wasserstein distance error bound between the posterior and its algorithmic approximation, we propose a post-processing step that yields an asymptotic, order-wise correct estimation of the posterior variance, eliminating the irreducible posterior variance estimation bias due to subsampling. Empirical performance of the method is evaluated through simulated experiments and a longitudinal study of pain trajectories in a study of breast cancer survivors.

The marginal likelihood, also known as the evidence, is regarded as a mathematical embodiment of Occam's razor, enabling hyperparameter selection that avoids overfitting. The evidence lower bound (ELBO) objective from variational inference has also been used for similar purposes. In this paper, we define a signal-to-noise ratio (SNR) for Bayesian linear regression in terms of key hyperparameters controlling signal and noise variance. We then study how ELBO-based hyperparameter learning can under or overestimate SNR depending on the assumed covariance matrix structure in a Gaussian approximate posterior. We prove that assuming a diagonal or low-rank plus diagonal covariance matrix underestimates the SNR, while a low-rank covariance matrix overestimates the SNR. Underestimating the SNR leads to higher than optimal mean squared error on the training set (underfitting), while overestimating SNR leads to lower error (overfitting). Bayesian practitioners hoping to scale to large models should be cautious about how covariance assumptions needed for tractability may impact model selection.

Auxiliary information is frequently utilized in survey sampling to improve the efficiency of estimators of the finite population mean. However, the simultaneous use of multiple auxiliary variables often induces multicollinearity, which adversely affects the stability and performance of conventional estimators. To address this issue, the present study proposes a principal component analysis (PCA) based estimation approach for the finite population mean in the presence of multicollinearity between two auxiliary variables. The proposed methodology transforms the correlated auxiliary variables into a set of orthogonal principal components, thereby removing the effect of multicollinearity while preserving the essential information contained in the auxiliary variables. An efficient estimator is then constructed using these components under simple random sampling without replacement. The bias and mean square error (MSE) of the proposed estimator are derived up to the first order of approximation. The performance of the proposed estimator is evaluated through both empirical and simulation studies under varying correlation structures. Moreover, the presence of multicollinearity is evaluated using variance inflation factors, condition indices, and eigenvalues. The results from empirical and simulation studies demonstrate that the proposed PCA-based estimator outperforms several conventional estimators in terms of MSE and percentage relative efficiency (PRE) when multicollinearity exists, ensuring robust and efficient estimation of the population mean.

Deep learning models are widely deployed in safety-critical domains, but remain vulnerable to adversarial attacks. In this paper, we study the adversarial robustness of NTK neural networks in the context of nonparametric regression. We establish minimax optimal rates for adversarial regression in Sobolev spaces and then show that NTK neural networks, trained via gradient flow with early stopping, can achieve this optimal rate. However, in the overfitting regime, we prove that the minimum norm interpolant is vulnerable to adversarial perturbations.

Identity teacher forcing (ITF) enables stable training of deterministic recurrent surrogates for chaotic dynamical systems and has been highly effective for dynamical systems reconstruction (DSR) with recurrent neural networks (RNNs), including interpretable almost-linear RNNs (AL-RNNs). However, as an intervention-based prediction loss (and thus a generalized Bayes update), teacher forcing need not match the free-running model's marginal likelihood geometry. We compare the objective-induced curvatures of ITF and marginal likelihood in a probabilistic switching augmentation of AL-RNNs, estimating ambiguity-aware observed information via Louis' identity. In the switching setting studied here, conditioning on a single forced regime path (as ITF does) inflates curvature, while marginal likelihood curvature is reduced by a missing-information correction when multiple switching explanations remain plausible. In Lorenz-63 experiments, windowed evidence fine-tuning improves held-out evidence but can degrade dynamical quantities of interest (QoIs) relative to ITF-pretrained models.

In this paper we develop a consistent variable selection procedure for GARCH-X models that identifies the truly relevant exogenous covariates influencing volatility dynamics. The proposed method is based on a multiple hypothesis testing framework with Wald-type test statistics and the Benjamini-Yekutieli False Discovery Rate (FDR) procedure to control the proportion of false discoveries. We establish the consistency of the selection rule, showing that it asymptotically recovers the correct set of covariates as the sample size increases. Monte Carlo simulations across different distributions and dependence structures validate the method's accuracy and robustness. The procedure is applied to modeling the volatility of the SP 500 using macroeconomic and commodity indicators.

Training language models via reinforcement learning often relies on imperfect proxy rewards, since ground truth rewards that precisely define the intended behavior are rarely available. Standard metrics for assessing the quality of proxy rewards, such as ranking accuracy, treat incorrect rewards as strictly harmful. In this work, however, we highlight that not all deviations from the ground truth are equal. By theoretically analyzing which outputs attract probability during policy gradient optimization, we categorize reward errors according to their effect on the increase in ground truth reward. The analysis establishes that reward errors, though conventionally viewed as harmful, can also be benign or even beneficial by preventing the policy from stalling around outputs with mediocre ground truth reward. We then present two practical implications of our theory. First, for reinforcement learning from human feedback (RLHF), we develop reward model evaluation metrics that account for the harmfulness of reward errors. Compared to standard ranking accuracy, these metrics typically correlate better with the performance of a language model after RLHF, yet gaps remain in robustly evaluating reward models. Second, we provide insights for reward design in settings with verifiable rewards. A key theme underlying our results is that the effectiveness of a proxy reward function depends heavily on its interaction with the initial policy and learning algorithm.

This article extends weak convergence bounds of Markov transition kernels to convergence bounds on the variance of the Markov kernel applied to Lipschitz functions. In the reversible case, weak convergence rates of the transition kernels imply chi-squared divergence convergence bounds if the density of the initialization measure is Lipschitz. These results provide new tools to establish central limit theorems for Lipschitz functions used in Markov chain Monte Carlo simulations. Applications are explored to the stability of Metropolis-Hastings algorithms in high dimensions, stochastic gradient descent, and solutions to stochastic delay equations.

Label noise presents a fundamental challenge in modern machine learning, especially when large-scale datasets are generated via automated processes. An increasingly common and important data paradigm, particularly in domains like medical imaging, involves learning from a large dataset with coarse, noisy labels supplemented by a small, expert-verified, clean dataset. This setting constitutes a typical information transfer and fusion problem. However, the significant distribution shift between the noisy and clean data violates the core overall parametric similarity assumptions of existing statistical transfer learning methods, while their reliance on parametric models is ill-suited for complex data like images. To address these limitations, this paper develops a generic model-agnostic nonparametric framework for classification with label noise, which applies to a broad class of classifiers. Our approach leverages the small clean dataset to “purify” the large noisy one and carefully manages the remaining ambiguous samples. This framework is underpinned by a rigorous statistical theory. Its empirical performance is demonstrated through simulations and a real-world application to medical image analysis for pneumonia diagnosis.

We show that widely used bubble tests, most prominently the PSY framework, suffer severe size distortion when fundamentals incorporate general-purpose technology adoption. Embedding a hump-shaped technology shock in the Campbell-Shiller present-value model, we prove that the fundamental price becomes locally explosive during adoption, thereby altering the asymptotic null distribution of the test statistic and causing the standard bubble test to overreject. We propose a technology-adjusted diagnostic that removes an estimated technology component from measures of productivity, IT-investment, and patents before testing the residual. The adjustment is conservative: because a boom can itself raise these technology measures, a rejection remains robust to such feedback, whereas a non-rejection only bounds residual explosiveness. Dot-com residual explosiveness concentrates in December 1999-March 2000; the 2020-2025 AI rally shows no residual explosiveness in our sample across baseline and sensitivity checks.

Magnification shift is a major obstacle to robust histopathology classification, because models trained on one imaging scale often generalize poorly to another. Here, we evaluated this problem on the BreaKHis dataset using a strict patient-disjoint leave-one-magnification-out protocol, comparing supervised baseline, baseline augmented with DCGAN-generated patches, and a gradient-reversal domain-general model designed to preserve discriminative information while suppressing magnification-specific variation. Across held-out magnifications, the domain-general model achieved the strongest overall discrimination and its clearest gain was observed when 200X was held out. By contrast, GAN augmentation produced inconsistent effects, improving some folds but degrading others, particularly at 400X. The domain-general model also yielded the lowest Brier score at 0.063 vs 0.089 at baseline. Sparse embedding analysis further revealed that domain-general training reduced average signature size more than three-fold (306 versus 1,074 dimensions) while preserving equivalent predictive performance (AUC: 0.967 vs 0.965; F1: 0.930 vs 0.931). It also increased cross-fold signature reproducibility from near-zero Jaccard overlap in the baseline to 0.99 between the 100X and 200X folds. These findings show that calibrated, compact, and transferable representations can be learned without added architectural complexity, with clear implications for the reliable deployment of computational pathology models across heterogeneous acquisition settings.

The availability of data from multiple heterogeneous environments has motivated methods that remain reliable under distributional shifts. When the joint distribution of response and predictors varies across environments, the response may still depend on a subset of predictors through an invariant mechanism. Existing methods typically assess candidate invariant sets through pooled stability criteria, treating environmental variation as nuisance. In this paper, we propose a Bayesian framework that explicitly separates a shared response mechanism from environment-specific or response-dependent associations, exploiting heterogeneity as evidence for structure learning. A competitive spike-and-slab prior is designed to force each predictor to compete between invariant and non-invariant spurious effects. Under a tractable working model, we establish invariant model selection consistency and posterior contraction for invariant coefficients. We further study the presence of irrelevant predictors, characterize posterior concentration on an equivalent invariant class, and introduce a post-selection refinement that consistently recovers the minimal invariant model. Simulations and a real application illustrate the robustness and finite-sample efficiency of the proposed method.

We test the hypothesis that simulataneous linear contrasts of multiple variance components equal zero in a Gaussian variance components model via a parametric bootstrap. Applications include but are not limited to nested and crossed designs. The main technical contributions are a computationally efficient decomposition of the normalized residual log-likelihood that does not require the variance components to be non-negative or variance design matrices to be positive semi-definite, a modified Newton method for its minimization, and a method for efficient optimization and sampling under the null hypothesis that certain linear combinations of variance components equal zero. A special case of the proposed procedure is a test for multiple variance components simulataneously equalling zero, for which a likelihood ratio test was not previously available. However, the proposed procedure is significantly more general.

In the last few decades, Markov chain Monte Carlo (MCMC) methods have been widely applied to Bayesian updating of structural dynamic models in the field of structural health monitoring. Recently, several MCMC algorithms have been developed that incorporate neural networks to enhance their performance for specific Bayesian model updating problems. However, a common challenge with these approaches lies in the fact that the embedded neural networks often necessitate retraining when faced with new tasks, a process that is time-consuming and significantly undermines the competitiveness of these methods. This paper introduces a newly developed adaptive meta-learning stochastic gradient Hamiltonian Monte Carlo (AM-SGHMC) algorithm. The idea behind AM-SGHMC is to optimize the sampling strategy by training adaptive neural networks, and due to the adaptive design of the network inputs and outputs, the trained sampler can be directly applied to various Bayesian updating problems of the same type of structure without further training, thereby achieving meta-learning. Additionally, practical issues for the feasibility of the AM-SGHMC algorithm for structural dynamic model updating are addressed, and two examples involving Bayesian updating of multi-story building models with different model fidelity are used to demonstrate the effectiveness and generalization ability of the proposed method.

Our research is closely related to ontological studies in mathematics. It provides crucial insights into the nature of decisions and strategies characterized by Markov moments. In a stopping game, a holistic decision-maker would evaluate comprehensive information by assessing the probabilities of various outcomes and their associated payoffs. This involves understanding the current state, historical data, and potential future scenarios. Such a decision-maker must also consider strategic interactions by anticipating and accounting for the strategies of other players. They must be flexible in adapting their strategy as the game evolves and able to integrate uncertainty by incorporating risk preferences and tolerances. They would perform scenario analysis to evaluate the impact of different stopping times under varying conditions. The goal of this modeling and its implementation in psychological practice is to introduce a novel method for assessing the state of players, leveraging deviations from rational strategies as diagnostic indicators of their psychological and decision-making profiles. The details of other models will be subject to contributed papers. The article presents the theoretical basis for combining various factors when modeling decision-making processes. The original title is "Rationality, Deviation, and Diagnosis: A Holistic Approach to Stopping Games" and will be used when it is possible to describe and interpret the results of the experiments we write about in the last section of the paper.

Based on discrete observations, we develop a test to infer if the volatility function $\sigma(\cdot)$ within the nonparametric Gaussian white noise model $dY_t = \sigma(t)dW_t$ is constant. The testing procedure is shown to be minimax-optimal and adaptive for infill asymptotics and these results entail that a deviation from the null hypothesis of constancy is best measured in terms of the ratio of $\sigma(t)$ and its $L^2$-average. The derivation of optimal constants requires the construction of hypotheses with height $h(b)$, where the parameter $b$ solves $F_n(b)=0$ for given functions $F_n$. Proving this equation to be solvable for each $n\in\mathbb{N}$ and establishing quantitative bounds of the solutions is built upon the implicit function theorem.

Sparse Optimal Scoring (SOS) reformulates linear discriminant analysis to enable feature selection through elastic net regularization, making it well-suited for high-dimensional settings where the number of features exceeds observations. Most existing SOS methods use deflation-based strategies that compute discriminant vectors sequentially, which can propagate errors and produce suboptimal solutions. We propose a novel approach that estimates all discriminant vectors simultaneously under an explicit global orthogonality constraint, which we call Deflation-Free Sparse Optimal Scoring (DFSOS). DFSOS combines Bregman iteration with orthogonality-constrained optimization, decomposing the problem into tractable subproblems for scoring vectors, discriminant vectors, and orthogonality enforcement. We establish convergence to stationary points of the augmented Lagrangian under mild conditions. Extensive experiments using synthetic data and real-world time series data demonstrate that DFSOS achieves classification accuracy comparable to or better than existing deflation-based methods. These results indicate that deflation-free approaches offer a robust and effective framework for sparse discriminant analysis in high-dimensional problems.

Nonlinear dynamical systems with regime transitions are typically described by ordinary differential equations with jumping parameters parameters. Traditional methods often treat change-point detection and parameter estimation as separate tasks, ignoring the inherent coupling between them. To address this, we propose residual-loss anomaly analysis of physics-informed neural networks, a unified framework that leverages dynamical consistency within the physics-informed learning paradigm. This approach jointly infers piecewise parameters and transition points under a single set of constraints. The method follows a two-stage strategy: First, local physical residuals are analyzed through overlapping subinterval decomposition. When a subinterval spans a true transition point, the residual exhibits a distinct structural elevation in noise-free conditions, which has a non-zero lower bound, enabling effective localization of potential transition intervals. Second, within our framework, change-point locations and piecewise parameters are integrated into a unified physical loss function for joint optimization, enabling simultaneous identification. Experiments on benchmark nonlinear dynamical systems, including Malthusian and logistic growth models, Van der Pol oscillator, Lotka-Volterra model and Lorenz system, demonstrate that the proposed method outperforms traditional decoupled approaches in both change-point localization and parameter estimation accuracy. This study provides an efficient, unified solution for structurally coupled inverse problems in nonlinear dynamical systems with regime switching.

We study the problem of estimating locations in time at which the level of technology in an economy changes when given a sequence of time ordered inputs and outputs. We approach the problem through the lens of nonparametric frontier analysis with frontiers that expand sharply and globally over time, and develop an offline change point detection procedure which achieves the minimax localization rates for the problem at hand up to logarithmic factors. We additionally give a simple method for constructing confidence intervals for the unobserved change point locations. Finally, we explain how the procedure can be modified to accommodate local changes in technology, meaning that efficiency gains are only realized for certain combinations of inputs. Simulation studies and real data examples are also presented to illustrate the practical value of our methods.

Generating realistic synthetic citation, patent, or component dependency networks is essential for benchmarking community detection, graph visualisation, and network data mining algorithms. We present a systematic comparison of generators of directed graphs that are nearly acyclic and have a ground-truth community structure. We evaluate 12 methods across 7 real citation networks and 26 metrics. We propose the practice of reversing the directions of edges in static generators to break cycles and induce a citation-like flow, which markedly improves the performance of a degree-corrected Stochastic Block Model. Our methodological approach to evaluating community detection benchmarks distinguishes between endogenous and exogenous mesoscopic similarities, with the latter proving more important. This distinction reveals that high-parameter models gain much of their advantage by memorising planted community statistics, which limits their fidelity when the evaluation is blind to those labels. Finally, we introduce the Citation Seeder (CS) algorithm, an iterative generator grounded in the Price-Pareto model of citation networks, with interpretable parameters and O(N+E) generation. CS achieves competitive results against the best-performing baselines while using up to four orders of magnitude fewer parameters, providing an interpretable account of a network's structure and, on a temporal hold-out, extrapolating its future in-degree distribution more faithfully than descriptive baselines.

We propose the covariate-balanced-and-adjusted response-adaptive randomization (CBARA) procedure for adaptive design in clinical trials, which integrates the complementary strengths of covariate-adjusted response-adaptive randomization (CARA) and covariate-adaptive randomization (CAR). The CBARA procedure updates the target allocation ratio according to observed responses and patient covariate profiles without requiring a correctly specified model, thereby retaining CARA's ethical and efficiency considerations while improving robustness. In addition, the CBARA procedure extends the CAR principle from fixed target allocation ratios to covariate-adjusted adaptive target allocation ratios, yet still pursues balance in treatment allocation with respect to covariate features. This integration is enabled by a newly defined imbalance vector and three interrelated components: the allocation function, parameter estimation and update mechanism. We establish the asymptotic properties of covariate imbalance and the estimators under the CBARA procedure. The results demonstrate that the CBARA procedure can improve balance for both observed and unobserved covariates while preserving the consistency of the allocation ratio. The theoretical analysis is developed through a pseudo-Markov chain framework, where a new discrepancy measure for transition kernels is introduced to handle the continuity of Poisson equation solutions with respect to parameters.

This paper studies the distributed adaptiveestimation problems for stochastic large regression modelswith an infinite number of parameters. By constructing a re-cursive local cost function, we propose a novel distributedrecursive least squares algorithm to estimate the unknownsystem parameters, where the growth rate of regressors'dimension is characterized by a non-decreasing positivefunction. The almost sure convergence of the proposedalgorithm is established under a cooperative excitationcondition, which incorporates the temporal information andthe spatial information to reflect the cooperative effectamong multiple agents. Moreover, we analyze the predic-tion error by establishing the asymptotic upper boundof the accumulated regret without any excitation condi-tions. The main difficulty of theoretical analysis lies in howto analyze properties of the product of non-independentand non-stationary random matrices, whose dimensionschange over time simultaneously. Some techniques, suchas stochastic Lyapunov function, double-array martingaletheory and algebraic graph theory, are employed to dealwith the above issue. Our theoretical results are derivedwithout imposing independence or stationarity assump-tions on the regression vectors, thereby not excluding thecorrelated feedback signals.

Posterior inference for Dirichlet process mixture models is analytically intractable and typically relies on Markov chain Monte Carlo methods, which can become computationally prohibitive at moderate to large sample sizes. In this work, we investigate the performance of Laplace and skew-Laplace posterior approximations for density estimation in this setting. Through an extensive numerical study covering four simulation scenarios with sample sizes ranging from n = 20 to n = 2,000 and four standard real datasets, we compare the standard Laplace approximation, its skew-corrected extension, and a slice sampling benchmark, assessing accuracy through total variation distance and computational efficiency through runtime. Our results show that the Gaussian Laplace approximation is more effective in this setting than might be anticipated, and that the skew-Laplace approximation consistently improves posterior recovery while remaining substantially faster than state-of-the-art Markov chain Monte Carlo samplers across all settings considered. In particular, the use of skew-Laplace in place of the standard Laplace approximation is especially beneficial in more complex density structures, where we observe error reductions typically on the order of 30%.