Papers

We develop a finite-sample, design-based theory for random forests in which each tree is a randomized conditional predictor acting on fixed covariates and the forest is their Monte Carlo average. An exact variance identity separates Monte Carlo error from a covariance floor that persists under infinite aggregation. The floor arises through two mechanisms: observation reuse, where the same training outcomes receive weight across multiple trees, and partition alignment, where independently generated trees discover similar conditional prediction rules. We prove the floor is strictly positive under minimal conditions and show that alignment persists even when sample splitting eliminates observation overlap entirely. We introduce procedure-aligned synthetic resampling (PASR) to estimate the covariance floor, decomposing the total prediction uncertainty of a deployed forest into interpretable components. For continuous outcomes, resulting prediction intervals achieve nominal coverage with a theoretically guaranteed conservative bias direction. For classification forests, the PASR estimator is asymptotically unbiased, providing the first pointwise confidence intervals for predicted conditional probabilities from a deployed forest. Nominal coverage is maintained across a range of design configurations for both outcome types, including high-dimensional settings. The underlying theory extends to any tree-based ensemble with an exchangeable tree-generating mechanism.

We propose a probabilistic extension of Wiener-Laguerre models for causal operator learning. Classical Wiener-Laguerre models parameterize stable linear dynamics using orthonormal Laguerre bases and apply a static nonlinear map to the resulting features. While structurally efficient and interpretable, they provide only deterministic point estimates. We reinterpret the nonlinear component through the lens of Barron function approximation, viewing two-layer networks, random Fourier features, and extreme learning machines as discretizations of integral representations over parameter measures. This perspective naturally admits Bayesian inference on the nonlinear map and yields posterior predictive uncertainty. By combining Laguerre-parameterized causal dynamics with probabilistic Barron-type nonlinear approximators, we obtain a structured yet expressive class of causal operators equipped with uncertainty quantification. The resulting framework bridges classical system identification and modern measure-based function approximation, providing a principled approach to time-series modeling and nonlinear systems identification.

The ability of Flow Matching (FM) to model complex conditional distributions has established it as the state-of-the-art for prediction tasks (e.g., robotics, weather forecasting). However, deployment in safety-critical settings is hindered by a critical extrapolation hazard: driven by smoothness biases, flow models yield plausible outputs even for off-manifold conditions, resulting in silent failures indistinguishable from valid predictions. In this work, we introduce Diverging Flows, a novel approach that enables a single model to simultaneously perform conditional generation and native extrapolation detection by structurally enforcing inefficient transport for off-manifold inputs. We evaluate our method on synthetic manifolds, cross-domain style transfer, and weather temperature forecasting, demonstrating that it achieves effective detection of extrapolations without compromising predictive fidelity or inference latency. These results establish Diverging Flows as a robust solution for trustworthy flow models, paving the way for reliable deployment in domains such as medicine, robotics, and climate science.

Granger causality identifies predictive dependencies in multivariate time series. In distributed settings where parties cannot share data, federated causal learning enables joint analysis. Most federated causal methods assume that clients observe the same features and infer causal relationships as point estimates, with little formal uncertainty quantification. These assumptions do not hold in many industrial systems, where clients observe different features, and the objective is to estimate cross-client dependencies (edges). These dependencies must be estimated indirectly through repeated client-server iterations. Uncertainty from client data and model parameters propagates through this process, making point estimates alone insufficient for assessing cross-client edges. This paper characterizes this uncertainty propagation and uses edge-specific variances to distinguish genuine cross-client dependencies from spurious estimated edges. We consider aleatoric uncertainty from client data variability and epistemic uncertainty from model parameters. We derive closed-form variance recursions and steady-state variances for the client-server iterations. We prove that the propagated contribution of the initial model-parameter uncertainty vanishes asymptotically. These variances enable statistically principled selection of cross-client edges. Synthetic experiments show that our approach improves cross-client edge recovery over competing baselines. On real-world industrial datasets, it achieves high root-cause identification accuracy while yielding interpretable dependency structures.

In many randomized trials, outcomes such as essays or open-ended responses must be manually scored as a preliminary step to impact analysis, a process that is costly and limiting. Model-assisted estimation offers a way to combine surrogate outcomes generated by machine learning or large language models with a human-coded subset, yet typical implementations use simple random sampling and therefore overlook systematic variation in surrogate prediction error. We extend this framework by incorporating stratified sampling to more efficiently allocate human coding effort. We derive the exact variance of the stratified model-assisted estimator, characterize conditions under which stratification improves precision, and identify a Neyman-type optimal allocation rule that oversamples strata with larger residual variance. We evaluate our methods through a comprehensive simulation study to assess finite-sample performance. Overall, we find stratification consistently improves efficiency when surrogate prediction errors exhibit structured bias or heteroskedasticity. We also present two empirical applications, one using data from an education RCT and one using a large observational corpus, to illustrate how these methods can be implemented in practice using ChatGPT-generated surrogate outcomes. Overall, this framework provides a practical design-based approach for leveraging surrogate outcomes and strategically allocating human coding effort to obtain unbiased estimates with greater efficiency. While motivated by text-as-data applications, the methodology applies broadly to any setting where outcome measurement is costly or prohibitive, and can be applied to comparisons across groups or estimating the mean of a single group.

The one-dimensional Dickman distribution arises in various stochastic models across number theory, combinatorics, physics, and biology. Recently, a definition of the multidimensional Dickman distribution has appeared in the literature, together with its application to approximating the small jumps of multidimensional Lévy processes. In this paper, we extend this definition to a class of vector-valued random elements, which we characterise as fixed points of a specific affine transformation involving a random matrix obtained from the matrix exponential of a uniformly distributed random variable. We prove that these new distributions possess the key properties of infinite divisibility and operator selfdecomposability. Furthermore, we identify several cases where this new distribution arises as a limiting distribution.

Neuroimaging has profoundly enhanced our understanding of the human brain by characterizing its structure, function, and connectivity through modalities like MRI, fMRI, EEG, and PET. These technologies have enabled major breakthroughs across the lifespan, from early brain development to neurodegenerative and neuropsychiatric disorders. Despite these advances, the brain is a complex, multiscale system, and neuroimaging measurements are correspondingly high-dimensional. This creates major statistical challenges, including measurement noise, motion-related artifacts, substantial inter-subject and site/scanner variability, and the sheer scale of modern studies. This paper explores statistical opportunities and challenges in neuroimaging across four key areas: (i) brain development from birth to age 20, (ii) the adult and aging brain, (iii) neurodegeneration and neuropsychiatric disorders, and (iv) brain encoding and decoding. After a quick tutorial on major imaging technologies, we review cutting-edge studies, underscore data and modeling challenges, and highlight research opportunities for statisticians. We conclude by emphasizing that close collaboration among statisticians, neuroscientists, and clinicians is essential for translating neuroimaging advances into improved diagnostics, deeper mechanistic insight, and more personalized treatments.

We propose a training-free conditional sampling method for flow matching models based on importance sampling. Because a naïve application of importance sampling suffers from weight degeneracy in high-dimensional settings, we modify and incorporate a resampling technique in sequential Monte Carlo (SMC) during intermediate stages of the generation process. To encourage generated samples to diverge along distinct trajectories, we derive a stochastic flow with adjustable noise strength to replace the deterministic flow at the intermediate stage. Our framework requires no additional training, while providing theoretical guarantees of asymptotic accuracy. Experimentally, our method significantly outperforms existing approaches on conditional sampling tasks for MNIST and CIFAR-10. We further demonstrate the applicability of our approach in higher-dimensional, multimodal settings through text-to-image generation experiments on CelebA-HQ.

Mode collapse, the failure to capture one or more modes when targetting a multimodal distribution, is a central challenge in modern variational inference. In this work, we provide a mathematical analysis of annealing based strategies for mitigating mode collapse in a tractable setting: learning a Gaussian mixture, where mode collapse is known to arise. Leveraging a low dimensional summary statistics description, we precisely characterize the interplay between the initial temperature and the annealing rate, and derive a sharp formula for the probability of mode collapse. Our analysis shows that an appropriately chosen annealing scheme can robustly prevent mode collapse. Finally, we present numerical evidence that these theoretical tradeoffs qualitatively extend to neural network based models, RealNVP normalizing flows, providing guidance for designing annealing strategies mitigating mode collapse in practical variational inference pipelines.

The multi objective bandit setting has traditionally been regarded as more complex than the single objective case, as multiple objectives must be optimized simultaneously. In contrast to this prevailing view, we demonstrate that when multiple good arms exist for multiple objectives, they can induce a surprising benefit, implicit exploration. Under this condition, we show that simple algorithms that greedily select actions in most rounds can nonetheless achieve strong performance, both theoretically and empirically. To our knowledge, this is the first study to introduce implicit exploration in both multi objective and parametric bandit settings without any distributional assumptions on the contexts. We further introduce a framework for effective Pareto fairness, which provides a principled approach to rigorously analyzing fairness of multi objective bandit algorithms.

Continuous and strictly positive data that exhibit skewness and outliers frequently arise in many applied disciplines. Log-symmetric distributions provide a flexible framework for modeling such data. In this article, we develop new goodness-of-fit tests for log-symmetric distributions based on a recent characterization. These tests utilize the characteristic function as a novel tool and are constructed using an $L^2$-type weighted distance measure. The asymptotic properties of the resulting test statistic are studied. The finite-sample performance of the proposed method is assessed via Monte Carlo simulations and compared with existing procedures. The results under a range of alternative distributions indicate superior empirical power, while the proposed test also exhibits substantial computational efficiency compared to existing methods. The methodology is further illustrated using real data sets to demonstrate practical applicability.

Unimodal univariate distributions can be characterized as piecewise convex-concave cumulative distribution functions. In this note we transfer this shape constraint characterization to the quantile function. We show that this characterization comes with the upside that the quantile function of a unimodal distribution is always absolutely continuous and consequently unimodality is equivalent to the quasi-convexity of its Radon-Nikodym derivative, i.e., the quantile density. Our analysis is based on the theory of generalized inverses of non-decreasing functions and relies on a version of the inverse function rule for non-decreasing functions.

Auxiliary information is increasingly available from administrative and other data sources, but it is often incomplete and of non-probability origin. We propose a two-step small area estimation approach in which the first step relies on design-based model calibration and exploits a large non-probability source providing a noisy proxy of the study variable for only part of the population. A unit-level measurement-error working model is fitted on the linked overlap between the probability survey and the external source, and its predictions are incorporated through domain-specific model-calibration constraints to obtain approximately design-unbiased domain totals. These totals and their variance estimates are then used in a Fay-Herriot area-level model with exactly known covariates to produce empirical best linear unbiased predictors. The approach is demonstrated in three enterprise survey settings from official statistics by integrating probability sample data with (i) administrative records, (ii) a cut-off data source, and (iii) web-scraped online information. Empirical comparisons show consistent improvements in domain-level precision over direct estimation and over a Fay-Herriot benchmark that directly incorporates the proxy information as an error-prone covariate. These gains are achieved without modeling the selection mechanism of the non-probability sample.

Directional measurements such as wind directions are often recorded in a finite number of angular categories rather than as exact angles. When two such measurements are observed jointly, the resulting bivariate observations lie on a discrete torus. Commonly used bivariate circular models are formulated for continuous angular variables. Applying these models to categorical observations requires integrating their densities over regions corresponding to observed category pairs. We propose two parametric models defined directly on the discrete torus, with interpretable parameters for marginal locations and concentrations, and for dependence between the two circular variables. The models provide closed-form probability mass functions and trigonometric moments, which are used to show that, under certain conditions, the dependence parameter characterizes circular–circular correlation. Parameters are estimated by maximum likelihood, and the finite-sample performance is investigated through simulation. The proposed models are applied to three datasets of paired wind direction measurements recorded in 16 equally spaced compass directions at stations in India and compared with discretized versions of established continuous bivariate circular models. They provide competitive fits while allowing likelihood evaluation directly on the observed discrete support. The fitted models are also used to assess the dependence between the paired wind directions in each dataset.

Matrix valued time series (MaTS) and global vector autoregressive (GVAR) models both impose restrictions on the general VAR for multidimensional data sets, in order to bring down the number of parameters. Both models are motivated from a different viewpoint such that on first sight they do not have much in common. When investigating the models more closely, however, one notices many connections between the two model sets. This paper investigates the relations between the restrictions imposed by the two models. We show that under appropriate restrictions in both models we obtain a joint framework allowing to gain insight into the nature of GVARs from the viewpoint of MaTS.

We assume that we have multiple ordinal time series and we would like to specify their joint distribution. In general it is difficult to create multivariate distribution that can be easily used to jointly model ordinal variables and the problem becomes even more complex in the case of time series, since we have to take into consideration not only the autocorrelation of each time series and the dependence between time series, but also cross-correlation. Starting from the simplest case of two ordinal time series, we propose using copulas to specify their joint distribution. We extend our approach in higher dimensions, by approximating full likelihood with composite likelihood and especially conditional pairwise likelihood, where each bivariate model is specified by copulas. We suggest maximizing each bivariate model independently to avoid computational issues and synthesize individual estimates using weighted mean. Weights are related to the Hessian matrix of each bivariate model. Simulation studies showed that model fits well under different sample sizes. Forecasting approach is also discussed. A small real data application about unemployment state of different countries of European Union is presented to illustrate our approach.

We reformulate Optimal Transport Conditional Flow Matching (OT-CFM), a class of dynamical generative models, showing that it admits an exact proximal formulation via an extended Brenier potential, without assuming that the target distribution has a density. In particular, the mapping to recover the target point is exactly given by a proximal operator, which yields an explicit proximal expression of the vector field. We also discuss the convergence of minibatch OT-CFM to the population formulation as the batch size increases. Finally, using second epi-derivatives of convex potentials, we prove that, for manifold-supported targets, OT-CFM is terminally normally hyperbolic: after time rescaling, the dynamics contracts exponentially in directions normal to the data manifold while remaining neutral along tangential directions.

Estimating prognosis conditional on surviving an initial high-risk period is crucial in clinical research. Yet, standard metrics such as hazard ratios are often difficult to interpret, while mean-based summaries are sensitive to outliers and censoring. We propose a formal causal framework for estimating quantiles of residual lifetime among individuals surviving to a landmark time $t_0$. Our primary estimand, the "Observed Survivor Quantile Contrast" (OSQC), targets pragmatic prognostic differences within the observed survivor population. To estimate the OSQC, we develop a doubly robust estimator that combines propensity scores, outcome regression, and inverse probability of censoring weights, ensuring consistency under confounding and informative censoring provided that the censoring model is correctly specified and at least one additional nuisance model is correctly specified. Recognizing that the OSQC conflates causal efficacy and compositional selection, we also introduce a reweighting-based supplementary estimator for the "Principal Survivor Quantile Contrast" (PSQC) to disentangle these mechanisms under stronger assumptions. Extensive simulations demonstrate the robustness of the proposed estimators and clarify the role of post-treatment selection. We illustrate the framework using data from the SUPPORT study to assess the impact of right heart catheterization on residual lifetime among intensive care unit survivors, and from the NSABP B-14 trial to examine post-surgical prognosis under adjuvant tamoxifen therapy across multiple landmark times.

The rule of thumb regarding the relationship between the bias-variance tradeoff and model size plays a key role in classical machine learning, but is now well-known to break down in the overparameterized setting as per the double descent curve. In particular, minimum-norm interpolating estimators can perform well, suggesting the need for new tradeoff in these settings. Accordingly, we propose a regularization-sharpness tradeoff for overparameterized linear regression with an $\ell^p$ penalty. Inspired by the interpolating information criterion, our framework decomposes the selection penalty into a regularization term (quantifying the alignment of the regularizer and the interpolator) and a geometric sharpness term on the interpolating manifold (quantifying the effect of local perturbations), yielding a tradeoff analogous to bias-variance. Building on prior analyses that established this information criterion for ridge regularizers, this work first provides a general expression of the interpolating information criterion for $\ell^p$ regularizers where $p \ge 2$. Subsequently, we extend this to the LASSO interpolator with $\ell^1$ regularizer, which induces stronger sparsity. Empirical results on real-world datasets with random Fourier features and polynomials validate our theory, demonstrating how the tradeoff terms can distinguish performant linear interpolators from weaker ones.

In this paper, we consider procedures for testing hypotheses on the dimension of the linear span generated by a growing number of $p\times p$ covariance matrices from independent $q$ populations. Under a proper limiting scheme where all the parameters, $q$, $p$, and the sample sizes from the $q$ populations, are allowed to increase to infinity, we derive the asymptotic normality of the proposed test statistics. The proposed test procedures show satisfactory performance in finite samples under both the null and the alternative. We also apply the proposed many-sample dimensionality test to investigate a matrix-valued gene dataset from the Mouse Aging Project and gain some new knowledge about its covariance structures.

We present Unified Latent Dynamics (ULD), a novel reinforcement learning algorithm that unifies the efficiency of model-free methods with the representational strengths of model-based approaches, without incurring planning overhead. By embedding state-action pairs into a latent space in which the true value function is approximately linear, our method supports a single set of hyperparameters across diverse domains – from continuous control with low-dimensional and pixel inputs to high-dimensional Atari games. We prove that, under mild conditions, the fixed point of our embedding-based temporal-difference updates coincides with that of a corresponding linear model-based value expansion, and we derive explicit error bounds relating embedding fidelity to value approximation quality. In practice, ULD employs synchronized updates of encoder, value, and policy networks, auxiliary losses for short-horizon predictive dynamics, and reward-scale normalization to ensure stable learning under sparse rewards. Evaluated on 80 environments spanning Gym locomotion, DeepMind Control (proprioceptive and visual), and Atari, our approach matches or exceeds the performance of specialized model-free and general model-based baselines – achieving cross-domain competence with minimal tuning and a fraction of the parameter footprint. These results indicate that value-aligned latent representations alone can deliver the adaptability and sample efficiency traditionally attributed to full model-based planning.

Differential privacy provides a formal framework for releasing statistical estimators that limit how much any single observation can influence the output, by injecting calibrated random noise. We study differentially private estimation in two-stage procedures common in causal inference and individualized treatment rule (ITR) learning, in which data-dependent weights are first estimated to enforce covariate balance and a parameter of interest is then obtained by weighted empirical risk minimization. We propose Differentially Private Two-Stage Empirical Risk Minimization (DP-2ERM), which privatizes the final estimator directly through objective perturbation calibrated to the data-dependent sensitivity of the full pipeline. The analysis combines deterministic weight-perturbation bounds for several covariate-balancing methods (inverse propensity weighting, entropy balancing weighting, and maximum mean discrepancy weighting) with probabilistic sensitivity bounds for the second-stage solution. The resulting calibration is sharper than the natural stage-wise composition baseline, which the same sensitivity analysis supplies as a byproduct. Simulation studies and a benchmark application to ITR learning demonstrate the improved privacy–utility trade-off.

Self-attention is often viewed as probabilistic query-key lookup, motivating designs that preserve normalized attention scores and fixed positional semantics. We advocate a simpler and more unified perspective: an autoregressive attention head can be viewed as a dynamic two-layer MLP whose weights are instantiated from the context history. From this view, attention scores form an ever-growing hidden representation, and standard MLP activations such as ReLU or GLU naturally implement input-conditioned selection over a context-dependent memory pool rather than a probability distribution. Based on this formulation, we introduce HyperMLP and HyperGLU, which learn dynamic mixing in both feature space and sequence space, using a reverse-offset (lag) layout to align temporal mixing with autoregressive semantics. We provide theoretical characterizations of the expressivity and implications of this structure, and empirically show that HyperMLP/HyperGLU consistently outperform strong softmax-attention baselines under matched parameter budgets.

A powerful robust mean estimator introduced by Catoni (2012) allows for mean estimation of heavy-tailed data while achieving the performance characteristics of classical mean estimator for sub-Gaussian data. While Catoni's framework has been widely extended across statistics, stochastic algorithms, and machine learning, fundamental asymptotic questions regarding the Central Limit Theorem and rare event deviations remain largely unaddressed. In this paper, we investigate Catoni-type robust estimators in two contexts: (i) mean estimation for heavy-tailed data, and (ii) linear regression with heavy-tailed innovations. For the first model, we establish the Berry–Esseen bound and moderate deviation principles, addressing both known and unknown variance settings. For the second model, we demonstrate that the associated estimator is consistent and satisfies a multi-dimensional Berry-Esseen bound.

Heterogeneity in multinomial choice data is often accounted for using logit models with random coefficients. Such models are called "mixed", but they can be difficult to estimate for large datasets. We review current Bayesian variational inference (VI) methods that can do so, and propose a new VI method that scales more effectively. The key innovation is a step that updates efficiently a Gaussian approximation to the conditional posterior of the random coefficients, addressing a bottleneck within the variational optimization. The approach is used to estimate three types of mixed logit models: standard, nested and bundle variants. We first demonstrate the improvement of our new approach over existing VI methods using simulations. Our method is then applied to a large scanner panel dataset of pasta choice. We find consumer response to price and promotion variables exhibits substantial heterogeneity at the grocery store and product levels. Store size, premium and geography are found to be drivers of store level estimates of price elasticities. Extension to bundle choice with pasta sauce improves model accuracy further. Predictions from the mixed models are more accurate than those from fixed coefficients equivalents, and our VI method provides insights in circumstances which other methods find challenging.

In network epidemic models, controlling the spread of a disease often requires targeted interventions such as vaccinating high-risk individuals based on network structure. However, typical approaches assume complete knowledge of the underlying contact network, which is often unavailable. While network structure can be learned from observed epidemic dynamics, existing methods require long observation windows that may delay critical interventions. In this work, we show that full network reconstruction may not be necessary: control-relevant features, such as high-degree vertices (super-spreaders), can be learned far more efficiently than the complete structure. Specifically, we develop an algorithm to identify such vertices from the dynamics of a Susceptible-Infected-Susceptible (SIS) process. We prove that in an $n$-vertex graph, vertices of degree at least $n^\alpha$ can be identified over an observation window of size $\Omega (1/\alpha)$, for any $\alpha \in (0,1)$. In contrast, existing methods for exact network reconstruction requires an observation window that grows linearly with $n$. Simulations demonstrate that our approach accurately identifies super-spreaders and enables effective epidemic control.

In truncated linear regression, samples $(x,y)$ are shown only when the outcome $y$ falls inside a certain survival set $S^\star$ and the goal is to estimate the unknown $d$-dimensional regressor $w^\star$. This problem has a long history of study in Statistics and Machine Learning going back to the works of (Galton, 1897; Tobin, 1958) and more recently in, e.g., (Daskalakis et al., 2019; 2021; Lee et al., 2023; 2024). Despite this long history, however, most prior works are limited to the special case where $S^\star$ is precisely known. The more practically relevant case, where $S^\star$ is unknown and must be learned from data, remains open: indeed, here the only available algorithms require strong assumptions on the distribution of the feature vectors (e.g., Gaussianity) and, even then, have a $d^{\mathrm{poly} (1/\varepsilon)}$ run time for achieving $\varepsilon$ accuracy. In this work, we give the first algorithm for truncated linear regression with unknown survival set that runs in $\mathrm{poly} (d/\varepsilon)$ time, by only requiring that the feature vectors are sub-Gaussian. Our algorithm relies on a novel subroutine for efficiently learning unions of a bounded number of intervals using access to positive examples (without any negative examples) under a certain smoothness condition. This learning guarantee adds to the line of works on positive-only PAC learning and may be of independent interest.

Conditional Value-at-Risk (CoVaR) quantifies systemic financial risk by measuring the loss quantile of one asset, conditional on another asset experiencing distress. We develop a Transformer-based methodology that integrates financial news articles directly with market data to improve CoVaR estimates. Unlike approaches that use predefined sentiment scores, our method incorporates raw text embeddings generated by a large language model (LLM). We prove explicit error bounds for our Transformer CoVaR estimator, showing that accurate CoVaR learning is possible even with small datasets. Using U.S. market returns and Reuters news items from 2006–2013, our out-of-sample results show that textual information impacts the CoVaR forecasts. With better predictive performance, we identify a pronounced negative dip during market stress periods across several equity assets when comparing the Transformer-based CoVaR to both the CoVaR without text and the CoVaR using traditional sentiment measures. Our results show that textual data can be used to effectively model systemic risk without requiring prohibitively large data sets.

Randomized Kaczmarz (RK) is a simple and fast solver for consistent overdetermined systems, but it is known to be fragile under noise. We study overdetermined $m\times n$ linear systems with a sparse set of corrupted equations, $ {\bf A}{\bf x}^\star = {\bf b}, $where only $\tilde{\bf b} = {\bf b} + \boldsymbol{\varepsilon}$ is observed with $\|\boldsymbol{\varepsilon}\|_0 \le \beta m$. The recently introduced QuantileRK (QRK) algorithm addresses this issue by testing residuals against a quantile threshold, but computing a per-iteration quantile across many rows is costly. In this work we (i) reanalyze QRK and show that its convergence rate improves monotonically as the corruption fraction $\beta$ decreases; (ii) propose a simple online detector that flags and removes unreliable rows, which reduces the effective $\beta$ and speeds up convergence; and (iii) make the method practical by estimating quantiles from a small random subsample of rows, preserving robustness while lowering the per-iteration cost. Simulations on imaging and synthetic data demonstrate the efficiency of the proposed method.

Learning Gibbs distributions using only sufficient statistics has long been recognized as a computationally hard problem. On the other hand, computationally efficient algorithms for learning Gibbs distributions rely on access to full sample configurations generated from the model. For many systems of interest that arise in physical contexts, expecting a full sample to be observed is not practical, and hence it is important to look for computationally efficient methods that solve the learning problem with access to only a limited set of statistics. We examine the trade-offs between the power of computation and observation within this scenario, employing the Ising model as a paradigmatic example. We demonstrate that it is feasible to reconstruct the model parameters for a model with $\ell_1$ width $\gamma$ by observing statistics up to an order of $O(\gamma)$. This approach allows us to infer the model's structure and also learn its couplings and magnetic fields. We also discuss a setting where prior information about structure of the model is available and show that the learning problem can be solved efficiently with even more limited observational power.