arXane › Statistics › Statistics
Wednesday, 4 November 2026: 18 new papers, 7 cross-lists and 20 replacements in Methodology (stat.ME), in arXiv's order, with full abstracts. The same papers as arXiv's own listing, set as a newspaper; every title links to arXiv.
In this paper we study a regularized version of nearly-isotonic regression with respect to a general partial order. In this estimator the sum of squares is weighted, with some of the weights allowed to be zero, and an additional quadratic penalty is imposed. The only assumption we make is that the quadratic part of the objective function is strictly convex, and we show that this assumption is mild. First, we show that for all sufficiently large values of the penalty parameter the estimator coincides with the regularized isotonic regression, in which the penalty is replaced by the isotonic constraint. Further, we prove that the resulting problem can be written as a generalized lasso problem with a design matrix of full column rank and with the penalty matrix unchanged. Next, we derive the dual problem, which is a quadratic program with box constraints, and write it in a sparse form that can be solved by standard quadratic programming solvers for sparse problems.
Canopy cover is normally reported as a fraction of a horizontal map plane. When the domain of inference is the land surface, the appropriate target is the fraction of terrain surface area under canopy. We define this estimand on a C1 graph surface and show that it differs from vertical-projection cover by the planimetric covariance of cover with the area density, divided by the mean area density; the rumple index is a different quantity. A LiDAR pulse is recorded as a canopy return only when detection succeeds, so the naive return fraction estimates cover multiplied by the canopy detection probability. We study an inverse-probability-weighted estimator that uses that probability as the inclusion probability, and a weighted corner-averaged Cram\'er-von Mises lack-of-fit statistic for constant mean surface-area cover under independent canopy returns, calibrated by a parametric bootstrap. Synthetic validation on the flat plane and the Gaussian hill gives a naive bias near -0.40 and a weighted bias near 0.00; at 5000 pulses the test holds size at 0.02-0.06 across the simulated null cells and reaches power 0.79 under surface-uniform and 0.85 under plane-uniform sampling against broad slope-dependent cover, but only 0.12-0.13 against a compact hotspot. The definitions and their exact reductions are machine-checked in Lean 4. No field data yet.
Interval-censored multistate data arise when subjects are observed at irregular visit times, leaving transition times and intermediate states unobserved. Likelihood-based inference for continuous-time Markov models often requires computationally intensive matrix exponential evaluations. We propose a deterministic uniformization expectation-maximization (EM) algorithm for parameter estimation and latent trajectory inference. By representing the continuous-time process through an embedded Markov chain indexed by latent Poisson event counts, the method enables deterministic E-step calculations without repeated matrix exponential evaluations. The M-step reduces to transition-specific Poisson regressions with exposure offsets, accommodating covariate-dependent transition intensities. We establish uniform truncation error bounds and show that the procedure is an asymptotically monotone generalized EM algorithm. Under standard regularity conditions, the estimator is consistent and asymptotically normal. Simulations demonstrate accurate estimation and computational advantages over direct likelihood optimization, particularly for complex transition structures. We also develop an exact dynamic programming method to calculate posterior probabilities of candidate latent pathways and identify the most probable pathway between observations. Posterior expected transition counts and state occupancy times provide additional subject-specific summaries, while expected-count ratios offer a computational shortcut for classification in simple mutually exclusive pathway structures. The framework supports both population-level inference and individualized characterization of disease progression.
Sequential forecast intervals are often evaluated by aggregate coverage even when calibration memory is bounded. This can be misleading when economically important states recur after long absences: a recent-only buffer may satisfy its overall coverage target while no longer containing errors from earlier episodes of the returning state. We study this problem in next-day financial range forecasting, holding the forecasting model, conformity score, adaptive update, and memory budget constant while varying the historical residual information available at each query. Recent-only first-in, first-out (FIFO) calibration is near nominal in aggregate but undercovers substantially during high-volatility states in both an exchange-traded fund development panel and a fixed stock panel. An unbounded state-history reference yields a lower observed stress-period interval score; Rare replay yields lower high-VIX interval scores than the matched Recent-n and Hybrid-n controls in both panels, and the state-structured versus recent-only contrast is largest near stress re-entry. Bounded calibration memory is therefore part of the forecasting information design, and aggregate calibration alone can conceal state-specific reliability failures.
Model-based Small Area Estimation (SAE) has seen rapid growth over the past 25 years or so due to the rapidly increasing demand for reliable small area estimators of totals, means and more complex area parameters. Direct area-specific estimators are either unreliable (due to small area-specific sample sizes) or not feasible (due to zero area-specific sample sizes). It is therefore necessary to borrow strength from related areas through implicit or explicit models linking related areas. We focus on explicit area-level and unit-level models. After a brief historical sketch of the evolution of SAE, we focus on selective recent developments (mostly after 2015). Developments before 2015 are mostly covered in the comprehensive Wiley book by Rao and Molina (2015). Under a basic area-level model, we consider semi-parametric and robust estimation of small area means under model misspecifications. Under a basic unit level model, methods for estimating complex parameters, such as poverty indicators, will be presented, including robust estimation in the presence of outliers. A comparison of area-level and unit-level models when unit level data are available will also be presented. Estimation methods studied include empirical best linear unbiased prediction (EBLUP), empirical best (EB) and hierarchical Bayes (HB).
Hierarchical composite endpoints, analyzed through win statistics such as the net benefit (NB), win ratio (WR), and win odds (WO), are increasingly used in randomized trials to prioritize endpoints according to clinical importance. Although covariate-adjustment methods for win statistics have been developed, a unified model-assisted approach for mixed-type, time-restricted HCEs with simultaneous asymptotic no-loss efficiency guarantees remains of interest. We propose a unified, model-assisted covariate-adjustment framework for time-restricted marginal win statistics evaluated at a prespecified post-randomization restriction horizon, a clinically interpretable estimand particularly suited to acute and high-acuity trials with well-defined short-term follow-up windows. The Hoeffding--H'ajek projection of the two-sample win--loss U-statistics yields subject-level influence functions for the marginal win and loss probabilities. We apply treatment-specific linear covariate adjustment to these influence functions and propagate the adjustment to NB, log WR, and log WO. The resulting estimators are assumption-lean, remaining consistent and asymptotically valid without correct specification of the linear working model, and are asymptotically guaranteed efficiency gain over their unadjusted counterparts. Simulations with mixed-type HCEs demonstrate substantial efficiency and power gains with prognostic covariate adjustment. However, including additional non-prognostic covariates can impair finite-sample variance estimation and inflate Type I error, particularly in smaller trials. Bootstrap variance estimation improved inferential calibration in these settings. We illustrate the proposed method using a randomized acute kidney injury alert trial.
This paper develops upper L-statistic methods for testing and localizing a high-dimensional mean change. At each candidate split, the proposed statistics rank the squared studentized CUSUM values across coordinates and aggregate either a fixed number or a fixed proportion of the largest entries. This construction connects max-type and sum-type procedures and provides sensitivity to a wide range of sparsity regimes. We establish asymptotic null distributions under an independent-component representation allowing weak cross-sectional dependence. The limiting laws are Poisson tangent processes for fixed ranks and Gaussian-process or Gumbel laws for proportional ranks. To improve finite-sample performance, we develop analytic calibrations for marginal variance estimation, cross-sectional dependence, extreme clustering, and skewness, and establish their first-order validity. We combine calibrated p-values from a matched pair of fixed- and proportional-order statistics using the Cauchy combination rule to adapt to unknown sparsity. For localization, the estimated change point is selected from the component scan with the smallest p-value. Simulation studies demonstrate generally satisfactory size calibration, competitive detection power, and effective localization across different sparsity levels. An application to semiconductor manufacturing data further illustrates the practical utility of the proposed methods.
Listwise deletion (LWD) appears more common than imputation: our review of leading journals finds that about 50% of articles state or imply LWD, versus 7% describing an imputation model. Assuming oracle regression with all required values observed consistently estimates the population average treatment effect, we ask when complete-case regression remains consistent. With constant effects, we offer two checks: whether covariates and the missing indicator satisfy the adjustment criterion in an augmented causal diagram (DAG), and whether the estimator pools valid within-covariate comparisons appropriately, potentially using inverse treatment-propensity weighting (IPW). For heterogeneous effects, assuming conditional effects carry over to the population, we reweight comparisons to the population covariate distribution. Simulations and empirical illustrations from American politics, comparative politics, and international relations demonstrate the framework's use and show why a missing-at-random (MAR) classification, although useful for assessing imputation, does not by itself establish whether complete-case regression consistently estimates the target population ATE. Thus, we recommend adding the missing indicator and its causes to causal models, then checking adjustment and weighting for LWD.
Dirichlet kernels provide a natural tool for density estimation on the simplex, but data-driven bandwidth selection for them has been mostly limited to likelihood cross-validation (LCV). We develop bandwidth selectors for Dirichlet kernel estimators of densities and projected first derivatives for compositional data. Using a Dirichlet mixture as the reference model, we derive an exact representation of the mean integrated squared error (MISE) and define the exact MISE (EMI) selector. We also derive a closed-form least-squares cross-validation (LSCV) criterion and introduce the reference-smoothed least-squares (RSLS) selector, which replaces the unstable leave-one-out term of LSCV with the fitted Dirichlet mixture. For derivative estimation, we obtain projected derivatives of the Dirichlet density and their integrability conditions, extend the new bandwidth selectors, and carry out kernel mean shift clustering on the simplex. Simulations show that EMI is the most competitive selector for larger samples and Dirichlet-like targets, RSLS generally performs best for small samples and is more stable than LSCV for derivative estimation, and LCV remains competitive for non-Dirichlet targets in density estimation. Applications to wholesale customer expenditure and sectoral employment data illustrate density estimation and clustering, respectively.
This paper investigates the determinants of head coach turnover in the Italian Serie A over the period 2010-2019 by explicitly accounting for multiple exit mechanisms within a competing risks framework. Four mutually exclusive exit types are considered: end of contract, mutual consent, resignation, and sacking. To identify robust predictors of turnover, we combine classical cause-specific Cox models with a variable selection procedure based on forward selection, LASSO regression, and Random Survival Forests. The results show that performance-related variables are the primary drivers of coach turnover, while personal characteristics play a limited role once performance is controlled for. Short-term performance indicators, such as points per match until turnover, strongly influence the risk of dismissal. Conversely, cumulative poor performance substantially increases dismissal risk: a higher proportion of losses over a coach's tenure raises the hazard of sacking by more than four times. Season timing also plays a crucial role, as sackings and mutual-consent terminations are significantly more likely to occur during the season, whereas contract expirations concentrate at the end of the season. Overall, the findings highlight that clubs differentiate between short-term performance shocks and longer-term performance trajectories when deciding how and when to replace a head coach. By distinguishing among exit mechanisms and adopting a conservative variable selection strategy, the study provides a robust and nuanced assessment of managerial turnover dynamics in the most important Italian professional football league.
This paper aims to develop a statistical method for integrating multiple incomplete data sources to improve the estimation of the failure-rate probability distribution of fusion-specific devices. The data aggregation approach proposed in this work is to support the Reliability, Availability, Maintainability, and Inspectability (RAMI) studies for fusion-oriented research facilities. This work focuses on developing a procedure for handling missing or incomplete data from multiple sources, as well as formulating a statistical framework for combining data from different sources that refer to the same equipment type and failure mode. Specifically, a Bayesian result is adopted under the assumption of a constant failure rate over time. This framework motivates the use of the Gamma distributions to model the failure-rate parameter of each device. Subsequently, two alternative procedures are introduced to approximate, using a single Gamma distribution, the weighted mixture formed by combining Gamma distributions corresponding to the same equipment type across multiple data sources. The proposed method is validated through an extensive simulation study and further illustrated by application to the accelerator systems of the IFMIF-DONES facility. Finally, the proposed algorithm has been implemented in the open-source software R to facilitate its use in the nuclear fusion community.
Distributed statistical learning has recently attracted significant attention owing to its computational efficiency and scalability in large-scale systems. However, the non-smooth quantile loss poses substantial challenges to the development of efficient decentralized algorithms for quantile regression. Existing methods may either suffer from slow convergence or rely on smoothing-based approximations that require additional stringent assumptions, particularly in high-dimensional settings. To address these issues, we propose a decentralized framework that integrates gradient tracking with iterative hard thresholding for sparse quantile estimation. The proposed method avoids smoothing approximations and relaxes the initialization requirement. Theoretically, we show that the proposed algorithm achieves linear convergence under stated regularity conditions and that the resulting estimator attains the same order of statistical accuracy as its centralized counterpart. We further develop a decentralized information criterion for selecting the working sparsity level and establish support screening guarantees. Numerical experiments and a real-data study illustrate the finite-sample performance of the proposed method.
After variable selection, uncertainty quantification often treats the selected model as fixed, although it may change under repeated sampling. We investigate whether this variability can be incorporated by repeating Cox-Lasso selection within a model-based bootstrap after unpenalised Cox refitting. We compare conventional post-selection Cox Wald intervals and a model-based bootstrap without re-selection, both of which keep the selected predictor set fixed, with two re-selection procedures: one retains the observed tuning-parameter value, whereas the other reapplies the cross-validation minimum (CV-min) tuning rule within each bootstrap sample. Against the data-generating coefficients, procedures holding the selected predictor set fixed frequently showed substantial undercoverage. Re-selection generally moved coverage closer to the nominal level, although performance varied across simulation regimes and coefficients. Sensitivity analyses evaluating coverage against the population Cox coefficients corresponding to each originally selected predictor set gave similar aggregate coverage for the primary interval comparison. In the NKI70 breast-cancer data, the four selected genes appeared together in only 0.6% of bootstrap samples under re-selection at the observed tuning-parameter value. These results suggest that selection instability can be an important component of post-selection coefficient uncertainty and that bootstrap re-selection can help incorporate this variability.
Reporting standardised mean differences of continuous psychometric measures is omnipresent in psychology. In a randomised controlled trial with a multilevel structure, one can define a cluster-adjusted version of these quantities naturally taking into account heterogeneous residual variances and accounting for the multilevel structure. In this article, we carefully define the underlying location-scale distributional mixed-effects models and derive cluster-adjusted standardised mean differences using partial pooling. Estimating these quantities is best handled using Bayesian statistics, and we perform an extensive Monte Carlo simulation study to test the performance of these metrics under different simulation configurations. We conclude that in order to measure the a-priori heterogeneous treatment effect across clusters in practically relevant scenarios, it is enough to consider random effects only on the location model and not on the scale model. Even under artificial severe heteroscedasticity, this reduced model performs well both in terms of parameter bias and predictive accuracy. The resulting cluster-adjusted standardised mean difference is an average value across clusters, where the heterogeneity is summarised via its standard deviation.
Randomized controlled trials are crucial for assessing new treatments, but they are often infeasible due to ethical or practical constraints. In these cases, single-arm trials combined with external control data constitute an alternative approach for estimating treatment effects. However, individual patient-level data on external controls are frequently inaccessible due to privacy and economic concerns. Therefore, investigators rely on data summaries published by previous studies of the control therapy. We propose a method for inference on the treatment effect of an experimental therapy using patient-level data from a single-arm trial and a few published data summaries from an external control group. We introduce and study the Interval of Potential Effects (IPE), described as the range of plausible levels of the treatment effect that are compatible with the patient-level data from the experimental treatment and the available summaries on the control group. The width of the IPE quantifies the uncertainty on the treatment effect attributable to the lack of individual patient data from the control group. We demonstrate, through simulations and a case study, how the IPE indicates whether the available summaries are sufficient for accurate estimation of the treatment effect. We also illustrate that IPEs can support decision-making in the early phases of development of new treatments.
Prioritizing causal variants in Genome-Wide Association Studies (GWAS) is challenging due to linkage disequilibrium and high dimensionality. We propose an iterative algorithm for the selection of genetic variables that allows the integration of biological priors. The proposed algorithm requires the specification of the power function of a test statistic and incorporates a priori information on each SNP. We derive the approximate and exact distributions of several commonly used test statistics. At each step of the proposed algorithm, we select the SNP with the largest a posteriori probability of being causal. We consider several stopping criteria for the algorithm. We conduct comprehensive simulation studies to assess the performance of several versions of the proposed algorithm and to make a comparison to some existing competitors. The use of the proposed prioritization algorithm is illustrated with a cis-eQTL application to a human gene-expression data, using European discovery and Yoruba cross-population replication.
The Ising model, a type of Markov Random Field, from statistical physics has a long history of modeling spatial and network data as well as particle interactions. In this work, we consider inserting the Ising model into the framework of a Hidden Markov Model, a type of Bayesian Network, to model temporal state changes in a sequence of Ising observations. While this is entirely a discrete probability model, a main challenge is the fitting of such a probabilistic graphical model due to the exponential intractability of normalizing constants within the Ising model. To circumvent this problem, we consider the usage of maximum pseudolikelihood estimators that have a long history in Markov Random Fields. Validity of the pseudolikelihood Baum-Welch algorithm, model identifiability under likelihood, and consistency for large networks is also proved. The proposed Hidden Markov Ising Model is further extended to incorporate exogenous variables resulting in a hybrid hidden Markov autologistic model with graphical dependencies integrated into the model. This model is applied to a five year sequence of weekly gun violence events across the 3109 counties of the continental United States to identify state changes in gun violence activity from 2021 to 2025.
Importance sampling is a building block of modern statistical computing, but the weighted estimate is prone to high or infinite variance. Moreover, its empirical standard error is often unreliable, leaving importance sampling without practical uncertainty quantification. Motivated by multilevel regression and poststratification (MRP) in survey sampling, we propose post-stratified importance sampling (PIS) for estimating target expectations from proposal draws. PIS reformulates arbitrarily high-dimensional weighting into two one-dimensional probabilistic modeling problems: regression for the conditional mean given the importance ratio and marginal density estimation for the importance ratio. Different modeling choices yield a family of PIS estimators, encompassing ordinary importance sampling as a special case. We develop a model-based implementation that combines one-dimensional Gaussian process regression for the conditional mean and score-informed heavy-tailed density estimation for the marginal law, using both the observed weights and the available density evaluations. We study the oracle gains of PIS and establish consistency under general and model-specific conditions. Simulations and applications to cross-validation demonstrate improvements in estimation accuracy and uncertainty quantification.
We introduce the Axiomatic Learning Structures (ALS) and Axiomatic Probabilistic Structure (APS) in order to study the interactions of Machine Learning, Optimization, with research frameworks- which, we take the view, can be characterized at the core by the corresponding axioms. The ALS explores the interplays between structures that one is interested in Learning Theory such as Loss functions and the axioms of the objects and transformations involved in these functions. The APS studies the interaction between structures of interest in Probability Theory such as Gibbs distributions, equations with noise, etc. and the axioms of the corresponding objects and transformations. We also introduce in Section 6 Transformations-Axiomatic Based Learning (TABLe) on Mathematical data-sets. One feature of this approach is the ranking of transformations from a table for a certain task based on scores which reflect which transformations work better than others, similarly to when human investigators try various methods on a task. We also study structure recovery from noisy data. In particular, to exemplify the ideas, we consider loss functions which contain derivations and derivations-like transformations. As a first application we study recovery in the context of polynomials with integer coefficients using standard derivatives along with other types of derivations which appear in Arithmetic: the p-derivations which satisfy \textit{different axioms} from the usual derivatives. As a first application of the APS we study a class of noisy (partial/mixed) derivations equations. In the last sections, to provide further potential applications of this framework, we give an overview of an initial simplified Opinion Dynamics model and its (simultaneous) study of rules of updates.
We establish minimax separation rates across sparsity levels for regression signal detection under $(\epsilon, \delta)$-central differential privacy, considering linear regression with an unknown error distribution and generalised linear models (GLMs) with an unknown intercept. Without sparsity restrictions, the separation rates for linear regression and GLMs with binary responses match the private Gaussian mean testing rate up to polylogarithmic factors. For GLMs with more than two possible response outcomes, additional information in the marginal response distribution permits faster detection of sufficiently small signals when the dimension exceeds the sample size $n$. Our upper bounds are based on a reduction to private mean testing, tests of dependence between the response and sparse covariate projections, and a response-only test that is adjusted for the unknown intercept. The resulting procedure adapts to unknown sparsity and matches our lower bounds up to polylogarithmic factors. We also construct polynomial-time tests that attain these rates when $\epsilon\gtrsim n^{-1/2}$, and compare their numerical performance with existing private procedures.
We study nonparametric estimation and forecasting for GARCH-type processes and propose an iterated neural estimator that alternates between estimating the recursion map and updating the unobserved conditional variances. Under mild conditions, we establish finite-sample bounds for the time-averaged in-sample variance-path error and show that the estimator is minimax rate-optimal up to logarithmic factors. Under stationarity and mixing conditions, we establish the same rate for recursion-map estimation and derive finite-sample bounds for out-of-sample variance forecasts and expected excess Gaussian quasi-likelihood (QLIKE) loss. The convergence rates depend on the true lag orders and remain unchanged when the fitted orders are overspecified. Simulation studies illustrate the effectiveness of our estimators under nonlinear and time-varying volatility dynamics. Real-data applications suggest that our neural methods perform competitively and can incorporate option-price and high-frequency information.
Many signal-processing objects are defined by periodicity, normalization, sign or global-phase ambiguity, basis invariance, or positive-semidefinite constraints and therefore cannot be modeled faithfully in ordinary Euclidean space. This paper develops a geometry driven framework for directional statistics in signal processing, connecting signal representations to their natural manifolds or quotient spaces and to probability models compatible with their invariances. After reviewing instances of directional statistics in the signal processing literature, we review circular, toroidal, spherical, projective, Stiefel, Grassmann, and covariance-valued distributions, emphasizing their statistical assumptions, computational requirements, and dimensional limitations relevant to signal processing. We describe various modeling and estimation strategies, including analytical and gradient-based maximum likelihood, structured low-rank parameterizations, mixtures and hidden Markov models (HMMs), score-based and contrastive estimation, and geometry-aware K-means. To exemplify the consequences of geometry-aware modeling, we study multivariate phase coherence in synthetic data and functional neuroimaging, comparing toroidal phase vectors to projective, Grassmann, and positive-semidefinite representations. We also study speaker separation in a car environment, showing that directional representations outperform a standard Euclidean model. These examples demonstrate that when the geometric representation, statistical distribution, inference method, and computational regime are considered together, they have the potential to lead to improved models across a broad range of signal processing applications.
Smart heat meter data provide an important basis for monitoring district-heating installations, but fault reports are delayed, incomplete and temporally coarse: a report follows the appearance of a fault, sometimes by weeks, so the onset is unknown. This blurs the distinction between statistical changes and operational faults and prevents the reliable labelling that fault classifiers require. Changepoint detectors such as the Bayesian Estimator of Abrupt Change, Seasonality, and Trend (BEAST) report only the probability that a change occurred, so a large jump and a much smaller one can receive the same probability, and probabilities alone cannot tell a utility which installation to inspect first. We therefore develop a decision score that applies BEAST to temperature difference and flow and scores each detected changepoint by its posterior probability together with its deviation from recent operation, jump size, local behaviour and physical direction. Cases with unfavourable changes in both channels, the pattern of insufficient cooling, receive the highest priority. Of 1,522 cases, 648 had a reported fault; these made up 96% of the 25 highest-priority cases and 90% of the highest 200, and Average Precision was 0.771 against a reported-fault prevalence of 0.426. As cases without a report are not verified to be fault-free, these figures are likely conservative. Without using the report date, the strongest unfavourable flow change in reported-fault cases fell a median of 11 days before the report. Level-based screening of temperature difference, as practised by utilities, ranked the cases comparably but cannot date the change, and the two kinds of evidence were complementary. The method can serve as the first step in a fault-detection workflow, providing a manageable subset of installations for inspection and candidate onset dates for training fault classifiers.
Variable selection is an important task for interpreting potentially complex systems in fields such as biology, business, and sociology. The importance of variable selection has attracted considerable attention in statistics and machine learning. In this paper, we are mainly interested in model-agnostic variable selection for identifying variables that have general nonlinear dependence on the response. We propose MambaPFN, short for Model-Agnostic Markov BlAnket selection with Prior-data Fitted Networks, which is built upon the Transformer, a building block of many popular large language models. Specifically, we consider a functional-mapping perspective on variable-selection learning, in which the Transformer is treated as a learnable mapping from a dataset to its relevant variables. To learn this mapping, MambaPFN introduces three new components: (i) a faithful, rich, and scalable scheme for generating pretraining data, (ii) an objective that targets each variable's posterior inclusion probability, and (iii) a modified Transformer architecture tailored to variable selection. Once pretrained, MambaPFN is applied to new datasets without parameter updates. Empirically, MambaPFN achieves desirable variable-selection performance in simulation studies and exhibits a scaling law---more pretraining datasets and larger models can continuously improve its performance. Applied to a Facebook comments dataset, MambaPFN selects the fewest variables and achieves the highest computational efficiency and predictive performance, improving upon the best baselines by about 20%.
End-to-end digital twins can be constructed by fitting separate stage models and connecting them after inference. We investigate the Frequentist prediction properties of a commonly used Monte Carlo simulation scheme that forwards complete predictive draws through such a chain when measured stage outputs contain independent normal errors, but the latent stage means are deterministic. Analytical two-stage examples show persistent asymptotic interval and coverage inflation and, with quadratic propagation, mean bias dependent on the parameters to be estimated. Simulations of four- and eight-stage linear-quadratic univariate chains with 7, 15 or 30 training runs demonstrate pronounced input-dependent distortions. A diagnostic using error-free latent training predictors shows that the effect is mainly not attributable to errors in variables, but to the MCSim approach even showing underestimation of coverage. We propose a model-based bootstrap: an initial chain forwards fitted means, bootstrap chains are refitted on fixed-center bootstrap data while forwarding their mean prediction, and a target-stage future error enters a studentized prediction root. Across the investigated scenarios, relative mean bias is substantially reduced, and coverage is closer to the nominal 90%. At sample size 30, the eight-stage bootstrap has mean width ratios of 1.16-1.42 relative to the known-parameter normal reference, compared with 4.11-4.97 for observed-to-observed MCSim. These results support improved endpoint prediction in the tested setting. Further investigations of the bootstrap scheme in the supplement for a bivariate 6-stage process show similar satisfactory results.