Papers

In this paper we focus on jointly estimating the edge probabilities for multi-layer networks. We define a novel multi-layer graphon, a ternary function in contrast to the bivariate graphon function in the literature by introducing an additional latent layer position parameter, which is model-free and covers a wide range of multi-layer networks. We develop a computationally efficient two-step neighborhood smoothing algorithm to estimate the edge probabilities of multi-layer networks, which requires little tuning and fully utilize the similarity across both network layers and nodes. Numerical experiments demonstrate the advantages of our method over the existing state-of-the-art ones. A real Worldwide Food Import/Export Network dataset example is analyzed to illustrate the better performance of the proposed method over benchmark methods in terms of link prediction.

Finite mixture models are widely used in econometric analyses to capture unobserved heterogeneity. This paper shows that maximum likelihood estimation of finite mixtures of parametric densities can suffer from substantial finite-sample bias in all parameters under mild regularity conditions. The bias arises from the influence of outliers in component densities with unbounded or large support and increases with the degree of overlap among mixture components. I show that maximizing the classification-mixture likelihood function, equipped with a consistent classifier, yields parameter estimates that are less biased than those obtained by standard maximum likelihood estimation (MLE). I then derive the asymptotic distribution of the resulting estimator and provide conditions under which oracle efficiency is achieved. Monte Carlo simulations show that conventional mixture MLE exhibits pronounced finite-sample bias, which diminishes as the sample size or the statistical distance between component densities tends to infinity. The simulations further show that the proposed estimation strategy generally outperforms standard MLE in finite samples in terms of both bias and mean squared errors under relatively weak assumptions. An empirical application to latent group panel structures using health administrative data shows that the proposed approach reduces out-of-sample prediction error by approximately 17.6% relative to the best results obtained from standard MLE procedures.

We study online change point detection for multivariate inhomogeneous Poisson point process time series. This setting arises commonly in applications such as earthquake seismology, climate monitoring, and epidemic surveillance, yet remains underexplored in the machine learning and statistics literature. We propose a method that uses low-rank matrices to represent the multivariate Poisson intensity functions, resulting in an adaptive nonparametric detection procedure. Our algorithm is single-pass and requires only constant computational cost per new observation, independent of the elapsed length of the time series. We provide theoretical guarantees to control the overall false alarm probability and characterize the detection delay under temporal dependence. We also develop a new Matrix Bernstein inequality for temporally dependent Poisson point process time series, which may be of independent interest. Numerical experiments demonstrate that our method is both statistically robust and computationally efficient.

We study concentration inequalities for structured weighted sums of random data, including (i) tensor inner products and (ii) sequential matrix sums. We are interested in tail bounds and concentration inequalities for those structured weighted sums under exchangeability, extending beyond the classical framework of independent terms. We develop Hoeffding and Bernstein bounds provided with structure-dependent exchangeability. Along the way, we recover known results in weighted sum of exchangeable random variables and i.i.d. sums of random matrices to the optimal constants. Notably, we develop a sharper concentration bound for combinatorial sum of matrix arrays than the results previously derived from Chatterjee's method of exchangeable pairs. For applications, the richer structures provide us with novel analytical tools for estimating the average effect of multi-factor response models and studying fixed-design sketching methods in federated averaging. We apply our results to these problems, and find that our theoretical predictions are corroborated by numerical evidence.

Breiman and Cutler's original Random Forest was designed as a unified ML engine – not merely an ensemble predictor. Their implementation included classification, regression, unsupervised learning, proximity-based similarity, outlier detection, missing value imputation, and visualization – capabilities that modern libraries like scikit-learn never implemented. RFX-Fuse (Random Forests X [X=compression] – Forest Unified Learning and Similarity Engine) delivers Breiman and Cutler's complete vision with native GPU/CPU support. Modern ML pipelines require 5+ separate tools – XGBoost for prediction, FAISS for similarity, SHAP for explanations, Isolation Forest for outliers, custom code for importance. RFX-Fuse provides a 1 to 2 model object alternative – a single set of trees grown once. Novel Contributions: (1) Proximity Importance – native explainable similarity: proximity measures that samples are similar; proximity importance explains why. (2) Dataset-specific imputation validation for general tabular data – ranking imputation methods by how real the imputed data looks, without ground truth labels.

In the analysis of multivariate spatial and univariate spatio-temporal data, it is commonly recognized that asymmetric dependence may exist, which can be addressed using an asymmetric (matrix or space-time, respectively) covariance function within a Gaussian process framework. This paper introduces a new paradigm for constructing asymmetric space-time covariances, which we refer to as "reflective asymmetric," by leveraging recently-introduced models for multivariate spatial data. We first provide new results for reflective asymmetric multivariate spatial models that extends their applicability. We then propose their asymmetric space-time extension, which come from a substantially different perspective than Lagrangian asymmetric space-time covariances. There are fewer parameters in the new models, one controls both the spatial and temporal marginal covariances, and the standard separable model is a special case. In simulation studies and analysis of the frequently-studied Irish wind data, these new models also improve model fit and prediction performance, and they can be easier to estimate. These features indicate broad applicability for improved analysis in environmental and other space-time data.

Since its introduction, Facebook Prophet has attracted positive attention from both classical statisticians and the Bayesian statistics community. The model provides two built-in inference methods: maximum a posteriori estimation using the L-BFGS-B algorithm, and Markov Chain Monte Carlo (MCMC) sampling via the No-U-Turn Sampler (NUTS). While exploring various time-series forecasting problems using Bayesian inference with Prophet, we encountered limitations stemming from the inability to apply alternative inference techniques beyond those provided by default. Additionally, the fluent API design of Facebook Prophet proved insufficiently flexible for implementing our custom modeling ideas. To address these shortcomings, we developed a complete reimplementation of the Prophet model in PyMC, which enables us to extend the base model and evaluate and compare multiple Bayesian inference methods. In this paper, we present our PyMC-based implementation and analyze in detail the implementation of different Bayesian inference techniques. We consider full MCMC techniques, MAP estimation and Variational inference techniques on a time-series forecasting problem. We discuss in details the sampling approach, convergence diagnostics, forecasting metrics as well as their computational efficiency and detect possible issues which will be addressed in our future work.

We prove an exponential separation in sample complexity between Euclidean and hyperbolic representations for learning on hierarchical data under standard Lipschitz regularization. For depth-$R$ hierarchies with branching factor $m$, we first establish a geometric obstruction for Euclidean space: any bounded-radius embedding forces volumetric collapse, mapping exponentially many tree-distant points to nearby locations. This necessitates Lipschitz constants scaling as $\exp(\Omega(R))$ to realize even simple hierarchical targets, yielding exponential sample complexity under capacity control. We then show this obstruction vanishes in hyperbolic space: constant-distortion hyperbolic embeddings admit $O(1)$-Lipschitz realizability, enabling learning with $n = O(mR \log m)$ samples. A matching $\Omega(mR \log m)$ lower bound via Fano's inequality establishes that hyperbolic representations achieve the information-theoretic optimum. We also show a geometry-independent bottleneck: any rank-$k$ prediction space captures only $O(k)$ canonical hierarchical contrasts.

Metastatic castration-resistant prostate cancer (mCRPC) is a highly aggressive disease with poor prognosis and heterogeneous treatment response. In this work, we developed and externally validated a visit-level 180-day mortality risk model using longitudinal data from two Phase III cohorts (n=526 and n=640). Only visits with observable 180-day outcomes were labeled; right-censored cases were excluded from analysis. We compared five candidate architectures: Long Short-Term Memory, Gated Recurrent Unit (GRU), Cox Proportional Hazards, Random Survival Forest (RSF), and Logistic Regression. For each dataset, we selected the smallest risk-threshold that achieved an 85% sensitivity floor. The GRU and RSF models showed high discrimination capabilities initially (C-index: 87% for both). In external validation, the GRU obtained a higher calibration (slope: 0.93; intercept: 0.07) and achieved an PR-AUC of 0.87. Clinical impact analysis showed a median time-in-warning of 151.0 days for true positives (59.0 days for false positives) and 18.3 alerts per 100 patient-visits. Given late-stage frailty or cachexia and hemodynamic instability, permutation importance ranked BMI and systolic blood pressure as the strongest associations. These results suggest that longitudinal routine clinical markers can estimate short-horizon mortality risk in mCRPC and support proactive care planning over a multi-month window.

Standard Bayesian Optimization (BO) assumes uniform smoothness across the search space an assumption violated in multi-regime problems such as molecular conformation search through distinct energy basins or drug discovery across heterogeneous molecular scaffolds. A single GP either oversmooths sharp transitions or hallucinates noise in smooth regions, yielding miscalibrated uncertainty. We propose RAMBO, a Dirichlet Process Mixture of Gaussian Processes that automatically discovers latent regimes during optimization, each modeled by an independent GP with locally-optimized hyperparameters. We derive collapsed Gibbs sampling that analytically marginalizes latent functions for efficient inference, and introduce adaptive concentration parameter scheduling for coarse-to-fine regime discovery. Our acquisition functions decompose uncertainty into intra-regime and inter-regime components. Experiments on synthetic benchmarks and real-world applications, including molecular conformer optimization, virtual screening for drug discovery, and fusion reactor design, demonstrate consistent improvements over state-of-the-art baselines on multi-regime objectives.

Randomized controlled experiments assess new policy impacts on performance metrics to inform launch decisions. Traditional approaches evaluate metrics independently despite correlations, and mixed results (e.g., positive revenue impact, negative customer experience) require manual judgment, hindering scalability. We propose a Bayesian decision-theoretic framework that systematically incorporates multiple objectives and trade-offs by comparing expected risks across decisions. Our approach combines experimenter-defined loss functions with observed evidence, using hierarchical models to leverage historical experiment learnings for prior information on treatment effects. Through real and simulated Amazon supply chain experiments, we demonstrate that compared to null hypothesis statistical testing, our method increases estimation efficiency via informative hierarchical priors and simplifies decision-making by systematically incorporating business preferences and costs for comprehensive, scalable decisions.

We consider the problem of graph matching for a sequence of graphs generated under a time-dependent Markov chain noise model. Our edgelighter error model, a variant of the classical lamplighter random walk, iteratively corrupts the graph $G_0$ with edge-dependent noise, creating a sequence of noisy graph copies $(G_t)$. Much of the graph matching literature is focused on anonymization thresholds in edge-independent noise settings, and we establish novel anonymization thresholds in this edge-dependent noise setting when matching $G_0$ and $G_t$. Moreover, we also compare this anonymization threshold with the mixing properties of the Markov chain noise model. We show that when $G_0$ is drawn from an Erdős-Rényi model, the graph matching anonymization threshold and the mixing time of the edgelighter walk are both of order $\Theta(n^2\log n)$. We further demonstrate that for more structured model for $G_0$ (e.g., the Stochastic Block Model), graph matching anonymization can occur in $O(n^\alpha\log n)$ time for some $\alpha<2$, indicating that anonymization can occur before the Markov chain noise model globally mixes. Through extensive simulations, we verify our theoretical bounds in the settings of Erdős-Rényi random graphs and stochastic block model random graphs, and explore our findings on real-world datasets derived from a Facebook friendship network and a European research institution email communication network.

This paper introduces a version of decoupling and randomization to establish concentration inequalities for double-indexed permutation statistics. The results yield, among other applications, a new combinatorial Hanson-Wright inequality and a new combinatorial Bennett inequality. Several illustrative examples from rank-based statistics, graph-based statistics, and causal inference are also provided.

Industrial image anomaly detection is a challenging problem owing to extreme class imbalance and the scarcity of labeled defective samples, particularly in few-shot settings. We propose BayPrAnoMeta, a Bayesian generalization of Proto-MAML for few-shot industrial image anomaly detection. Unlike existing Proto-MAML approaches that rely on deterministic class prototypes and distance-based adaptation, BayPrAnoMeta replaces prototypes with task-specific probabilistic normality models and performs inner-loop adaptation via a Bayesian posterior predictive likelihood. We model normal support embeddings with a Normal-Inverse-Wishart (NIW) prior, producing a Student-$t$ predictive distribution that enables uncertainty-aware, heavy-tailed anomaly scoring and is essential for robustness in extreme few-shot settings. We further extend BayPrAnoMeta to a federated meta-learning framework with supervised contrastive regularization for heterogeneous industrial clients and prove convergence to stationary points of the resulting nonconvex objective. Experiments on the MVTec AD benchmark demonstrate consistent and significant AUROC improvements over MAML, Proto-MAML, and PatchCore-based methods in few-shot anomaly detection settings.

Deep neural networks (DNNs) achieve remarkable performance on a wide range of tasks, yet their mathematical analysis remains fragmented: stability and generalization are typically studied in disparate frameworks and on a case-by-case basis. Architecturally, DNNs rely on the recursive application of parametrized functions, a mechanism that can be unstable and difficult to train, making stability a primary concern. Even when training succeeds, there are few rigorous results on how well such models generalize beyond the observed data, especially in the generative setting. In this work, we leverage the theory of stochastic Iterated Function Systems (IFS) and show that two important deep architectures can be viewed as, or canonically associated with, place-dependent IFS. This connection allows us to import results from random dynamical systems to (i) establish the existence and uniqueness of invariant measures under suitable contractivity assumptions, and (ii) derive a Wasserstein generalization bound for generative modeling. The bound naturally leads to a new training objective that directly controls the collage-type approximation error between the data distribution and its image under the learned transfer operator. We illustrate the theory on a controlled 2D example and empirically evaluate the proposed objective on standard image datasets (MNIST, CelebA, CIFAR-10).

Bayesian evidence evaluation becomes computationally prohibitive in high dimensions due to the curse of dimensionality and the sequential nature of sampling-based methods. We introduce SunBURST, a deterministic GPU-native algorithm for Bayesian evidence calculation that replaces global volume exploration with mode-centric geometric integration. The pipeline combines radial mode discovery, batched L-BFGS refinement, and Laplace-based analytic integration, treating modes independently and converting large batches of likelihood evaluations into massively parallel GPU workloads. For Gaussian and near-Gaussian posteriors, where the Laplace approximation is exact or highly accurate, SunBURST achieves numerical agreement at double-precision tolerance in dimensions up to 1024 in our benchmarks, with sub-linear wall-clock scaling across the tested range. In multimodal Gaussian mixtures, conservative configurations yield sub-percent accuracy while maintaining favorable scaling. SunBURST is not intended as a universal replacement for sampling-based inference. Its design targets regimes common in physical parameter estimation and inverse problems, where posterior mass is locally well approximated by Gaussian structure around a finite number of modes. In strongly non-Gaussian settings, the method can serve as a fast geometry-aware evidence estimator or as a preprocessing stage for hybrid workflows. These results show that high-precision Bayesian evidence evaluation can be made computationally tractable in very high dimensions through deterministic integration combined with massive parallelism.

The first law of geography is a cornerstone of spatial analysis, emphasizing that nearby and related locations tend to be more similar, however, defining what constitutes "near" and "related" remains challenging, as different phenomena exhibit distinct spatial patterns. Traditional local regression models, such as Geographically Weighted Regression (GWR) and Multiscale GWR (MGWR), quantify spatial relationships solely through geographic proximity. In an era of globalization and digital connectivity, however, geographic proximity alone may be insufficient to capture how locations are interconnected. To address this limitation, we propose a new multiscale local regression framework, termed M-SGWR, which characterizes spatial interaction across two dimensions: geographic proximity and attribute (variable) similarity. For each predictor, geographic and attribute-based weight matrices are constructed separately and then combined using an optimized parameter, alpha, which governs their relative contribution to local model fitting. Analogous to variable-specific bandwidths in MGWR, the optimal alpha varies by predictor, allowing the model to flexibly account for geographic, mixed, or non-spatial (remote similarity) effects. Results from two simulation experiments and one empirical application demonstrate that M-SGWR consistently outperforms GWR, SGWR, and MGWR across all goodness-of-fit metrics.

This study investigates component wise estimation of ordered variances of scale mixture of two normal distributions. For this study two special loss functions are considered namely squared error loss function and entropy loss function. We have derived the general improvement results and based on these results the estimators that outperform BAEE are obtained. Moreover under certain sufficient conditions a class of improved estimators is proposed for both loss functions. As a special case of scale mixture of normal distribution the results are applied to the multivariate t-distribution and obtained the improvement results. For this case a detailed numerical comparison is carried out which validates our theoretical findings.

Beta Rank Function (BRF) is a two-sided distribution characterized by a smooth peak and double powerlaw decay, widely used to model empirical data exhibiting deviations from pure power laws. In this paper, we introduce a novel two-step generative process that produces data exactly following the BRF distribution. The first step involves any mechanism generating a power-law distribution, while the second step applies a regressive redistribution process that reallocates resources from poorer to richer entities, thereby amplifying inequality. This approach represents the first analytic derivation of an exact BRF distribution from a generative mechanism. We validate the model through applications to income and urban population distributions. Beyond exact generation, this framework offers new insights into the systemic origins of deviations from power laws frequently observed in complex systems, linking rank distributions to underlying feedback and redistribution dynamics.

Network Meta-Analysis (NMA) is an increasingly popular evidence synthesis tool that can provide a ranking of competing treatments, also known as a treatment hierarchy. Treatment-Covariate Interactions (TCIs) can be included in NMA models to allow relative treatment effects to vary with covariate values. We show that in an NMA model that includes TCIs, treatment hierarchies should be created with a particular covariate profile in mind. We outline the typical approach for creating a treatment hierarchy in standard Bayesian NMA and show how a treatment hierarchy for a particular covariate profile can be created from an NMA model that estimates TCIs. We demonstrate our methods using a real network of studies for treatments of major depressive disorder.

Intricate interactions among firms, institutions, and spatial structures shape urban economic systems. In this study, we propose a framework based on three structural dimensions – abundance, diversity, and longevity (ADL) of economic units – as proxies of urban economic complexity and resilience. Using a decade of georeferenced firm-level data from Mexico City, we analyze the relationships among ADL variables using regression, spatial correlation, and time-series clustering. Our results reveal nonlinear dynamics across urban space, with powerlaw behavior in central zones and logarithmic saturation in peripheral areas, suggesting differentiated growth regimes. Notably, firm longevity modulates the relationship between abundance and diversity, particularly in periurban transition zones. These spatial patterns point to an emerging polycentric restructuring within a traditionally monocentric metropolis. By integrating economic complexity theory with spatial analysis, our approach provides a scalable method to assess the adaptive capacity of urban economies. This has implications for understanding informality, designing inclusive urban policies, and navigating structural transitions in rapidly urbanizing regions.

Processing high-volume, streaming data is increasingly common in modern statistics and machine learning, where batch-mode algorithms are often impractical because they require repeated passes over the full dataset. This has motivated incremental stochastic estimation methods, including the incremental stochastic Expectation-Maximization (EM) algorithm formulated via stochastic approximation. In this work, we revisit and analyze an incremental stochastic variant of the Majorization-Minimization (MM) algorithm, which generalizes incremental stochastic EM as a special case. Our approach relaxes key EM requirements, such as explicit latent-variable representations, enabling broader applicability and greater algorithmic flexibility. We establish theoretical guarantees for the incremental stochastic MM algorithm, proving consistency in the sense that the iterates converge to a stationary point characterized by a vanishing gradient of the objective. We demonstrate these advantages on a softmax-gated mixture of experts (MoE) regression problem, for which no stochastic EM algorithm is available. Empirically, our method consistently outperforms widely used stochastic optimizers, including stochastic gradient descent, root mean square propagation, adaptive moment estimation, and second-order clipped stochastic optimization. These results support the development of new incremental stochastic algorithms, given the central role of softmax-gated MoE architectures in contemporary deep neural networks for heterogeneous data modeling. Beyond synthetic experiments, we also validate practical effectiveness on two real-world datasets, including a bioinformatics study of dent maize genotypes under drought stress that integrates high-dimensional proteomics with ecophysiological traits, where incremental stochastic MM yields stable gains in predictive performance.

Indicator variograms and madograms are structural tools used in many disciplines of the natural sciences and engineering to describe random sets and random fields. To date, several necessary conditions are known for a function to be a valid indicator variogram but, except for intractable corner-positive inequalities, a complete characterization of indicator variograms is missing. Likewise, only partial characterizations of madograms are known. This paper provides novel necessary and sufficient conditions for a given function to be the variogram of an indicator random field with constant mean value or to be the madogram of a random field, and establishes under which conditions these two families of functions coincide. Our results apply to any set of points where the random field is defined and rely on distance geometry and Gaussian random field theory.

We study grokking, the onset of generalization long after overfitting, in a classical ridge regression setting. We prove end-to-end grokking results for learning over-parameterized linear regression models using gradient descent with weight decay. Specifically, we prove that the following stages occur: (i) the model overfits the training data early during training; (ii) poor generalization persists long after overfitting has manifested; and (iii) the generalization error eventually becomes arbitrarily small. Moreover, we show, both theoretically and empirically, that grokking can be amplified or eliminated in a principled manner through proper hyperparameter tuning. To the best of our knowledge, these are the first rigorous quantitative bounds on the generalization delay (which we refer to as the "grokking time") in terms of training hyperparameters. Lastly, going beyond the linear setting, we empirically demonstrate that our quantitative bounds also capture the behavior of grokking on non-linear neural networks. Our results suggest that grokking is not an inherent failure mode of deep learning, but rather a consequence of specific training conditions, and thus does not require fundamental changes to the model architecture or learning algorithm to avoid.

Cure rate models address survival data in which a proportion of individuals will never experience the event of interest. Existing parametric approaches are predominantly based on finite mixtures, which impose restrictive assumptions on both the cure mechanism and the distribution of susceptible event times. A cure model based on phase-type distributions is introduced, leveraging their latent Markov jump process representation to allow immunity to occur either at baseline or dynamically during follow-up. This structure yields a flexible and interpretable formulation of long-term survival while encompassing classical mixture cure models as special cases. A unified regression framework is developed for covariate effects on both the cure rate and the susceptible survival distribution, and the proposed model class is dense, reducing the impact of parametric misspecification. Estimation is performed via expectation-maximization algorithms, accompanied by an automatic model selection strategy. Simulation studies and a real-data example demonstrate the practical advantages of the approach.

The empirical success of deep learning is often attributed to deep networks' ability to exploit hierarchical structure in data, constructing increasingly complex features across layers. Yet despite substantial progress in deep learning theory, most optimization results sill focus on networks with only two or three layers, leaving the theoretical understanding of hierarchical learning in genuinely deep models limited. This leads to a natural question: can we prove that deep networks, trained with gradient-based methods and standard input-label pairs, can efficiently exploit hierarchical structure? In this work, we consider Random Hierarchy Models – a hierarchical context-free grammar introduced by arXiv:2307.02129 and conjectured to separate deep and shallow networks. We prove that, under mild conditions, a deep convolutional network can be efficiently trained to learn this function class. Our proof builds on a general observation: if intermediate layers can receive clean signal from the labels and the relevant features are weakly identifiable, then layerwise training each individual layer suffices to hierarchically learn the target function.

We propose a framework to construct practical kernel-based two-sample tests from the family of $f$-divergences. The test statistic is computed from the witness function of a regularized variational representation of the divergence, which we estimate using kernel methods. The proposed test is adaptive over hyperparameters such as the kernel bandwidth and the regularization parameter. We provide theoretical guarantees for statistical test power across our family of $f$-divergence estimates. While our test covers a variety of $f$-divergences, we bring particular focus to the Hockey-Stick divergence, motivated by its applications to differential privacy auditing and machine unlearning evaluation. For two-sample testing, experiments demonstrate that different $f$-divergences are sensitive to different localized differences, illustrating the importance of leveraging diverse statistics. For machine unlearning, we propose a relative test that distinguishes true unlearning failures from safe distributional variations.

Estimating health indicators for restricted sub-populations is a recurring challenge in epidemiology and public health. When survey data are used, Small Area Estimation (SAE) methods can improve precision by borrowing strength across domains. In many applications, however, outcomes are self-reported and affected by coarsening mechanisms, such as rounding and digit preference, that reduce data resolution and may bias inference. This paper addresses both issues by developing a Bayesian unit-level SAE framework for semi-continuous, coarsened responses. Motivated by the 2019 Italian European Health Interview Survey, we estimate smoking indicators for domains defined by the cross-classification of Italian regions and age groups, capturing both smoking prevalence and intensity. The model adopts a two-part structure: a logistic component for smoking prevalence and a flexible mixture of Lognormal distributions for average cigarette consumption, coupled with an explicit model for coarsening and topcoding. Simulation studies show that ignoring coarsening can yield biased and unstable domain estimates with poor interval coverage, whereas the proposed model improves accuracy and achieves near-nominal coverage. The empirical application provides a detailed picture of smoking patterns across region-age domains, helping to characterize the dynamics of the phenomenon and inform targeted public health policies.

Finding effective ways to exploit parallel computing to accelerate Markov chain Monte Carlo methods is an important problem in Bayesian computation and related disciplines. In this paper, we consider the zeroth-order setting where the unnormalized target distribution can be evaluated but its gradient is unavailable for theoretical, practical, or computational reasons. We also assume access to $m$ parallel processors to accelerate convergence. The proposed approach is inspired by modern zeroth-order optimization methods, which mimic gradient-based schemes by replacing the gradient with a zeroth-order stochastic gradient estimator. Our contribution is twofold. First, we show that a naive application of popular zeroth-order stochastic gradient estimators within Markov chain Monte Carlo methods leads to algorithms with poor dependence on $m$, both for unadjusted and Metropolis-adjusted schemes. We then propose a simple remedy to this problem, based on a random-slice perspective, as opposed to a stochastic gradient one, obtaining a new class of zeroth-order samplers that provably achieve a polynomial speed-up in $m$. Theoretical findings are supported by numerical studies.

Constructing efficient portfolios requires balancing expected returns with risk through optimal stock selection, while accounting for investor preferences. In a recent work by Paul and Kundu (2026), the fractional-order entropy due to Ubriaco was introduced as an uncertainty measure to capture varying investor attitudes toward risk. Building on this foundation, we introduce a novel normalized fractional order entropy aligned with investors' risk preferences that combines normalized fractional entropy with expected utility and variance. Risk sensitivity is modeled through the fractional parameter, interpolating between conservative or risk aversion and adventurous or high risk tolerance attitudes. Furthermore, the robustness and statistical significance of the fractional order entropy-based risk measure, termed normalized expected utility-fractional entropy (NEU-FE) and normalized expected utility-fractional entropy-variance (NEU-FEV) risk measures are explained with the help of machine learning tools, including Random forest, Ridge regression, Lasso Regression and artificial neural networks by using Indian stock market (NIFTY50). The results confirm that the proposed decision models support investors in making high-quality portfolio investments.