2026-09-07 | | Total: 44
We study finite-sample valid tests of the one-sided mean hypothesis $H_0:μ\leq 1$ against $H_1:μ>1$ for nonnegative random variables. To do so, we develop a leave-one-out dual certificate framework, where certain pointwise inequalities imply p-value validity under the conditional mean null $\mathbb{E}[X_i\mid\mathbf{X}_{-i}]\leq 1$, and which also gives conditions that allow combining dual certificates for p-values to show that their pointwise minimum is also a valid p-value. The framework proves finite-sample validity of Wang and Zhao's nonparametric likelihood-ratio statistic $T_{\mathrm{nplr}}$, yields a new p-value $T_{\mathrm{bin}+}$ extending the Clopper--Pearson binomial test to general nonnegative random variables, and shows that the pointwise minimum $\min\{T_{\mathrm{nplr}},T_{\mathrm{bin}+}\}$ is itself a valid and more powerful p-value. We establish sharp optimality results for such testing problems in two regimes: both $T_{\mathrm{nplr}}$ and $T_{\mathrm{bin}+}$ attain a universal detectability boundary for the null $H_0$ without moment or tail assumptions, and $T_{\mathrm{bin}+}$ attains a nonparametric power lower bound under $n^{-1/2}$-local alternatives to $H_0$. Efficient algorithms and numerical experiments demonstrate substantial finite-sample power gains over existing valid methods.
Consider $n$ independent random vectors sampled from uniform distribution on $(p-1)$-dimensional unit sphere. This paper investigates the limiting joint distribution for the minimum and the maximum values of their pairwise angles. It proves that the minimum and the maximum angles are asymptotically independent when both $n$ and $p$ tend to infinity, which solves an open problem raised in Cai, Fan and Jiang (2013) [\emph{Journal of Machine Learning Research} 14, 1837-1864]. Cai, Fan and Jiang (2013) obtained the limiting marginal distributions for both the minimum and the maximum angles under assumption $\lim_{n\to\infty}{\ln n}/p=β$ according to whether $β=0$, $β\in (0,\infty)$, or $β=\infty$. This paper presents unified limits for both joint distributions and marginal distributions regardless of the relative divergence rate of $n$ and $p$. The paper also derives the limiting distributions for some statistics based on the minimum and the maximum angles,
Bayesian pseudo-posteriors based on moment conditions, such as Bayesian empirical likelihood (EL) and Bayesian exponentially tilted empirical likelihood (ETEL), provide a robust route to Bayesian inference when the model is specified only through moment restrictions, but their computation is often prohibitive. In this work, we address two barriers to Bayesian pseudo-inference under moment restrictions. First, scalable sampling is challenging for tall data because the pseudo-likelihoods are not log-additive across observations, so standard likelihood-subsampling methods do not apply. Second, in simulation-based and latent-variable problems, the moment function itself may be defined through an intractable expectation, making full-data pseudo-posterior inference infeasible. We develop a new family of subsampled pseudo-posteriors that replaces the full-data pseudo-posterior with an aggregate of mini-batch pseudo-posteriors. We show that naive mini-batching induces a shifted-mixture distortion, yielding aggregates that are correctly centered but overly diffuse. To remove this distortion, we introduce moment-function-level control variates that align mini-batch moment equations around a reference estimator and recover the full-data posterior shape. For Bayesian ETEL, we establish finite-sample total variation bounds to the full-data pseudo-posterior, allowing discontinuous and intractable moment functions. Numerical experiments show that the method applies broadly across moment-condition-based Bayesian pseudo-inference, closely approximates full-data pseudo-posteriors, and substantially reduces computational cost.
Suppose that the eventual use of data is not known when the data are reduced or collected. This note considers two simple boundary cases. In a finite statistical experiment, a statistic preserves the Bayes risk for every finite later decision problem if and only if it is sufficient. Hence, when the minimal sufficient statistic is one-to-one, exact preservation of all later decision problems permits no nontrivial reduction. We then consider adaptive sampling from $m$ independent Gaussian streams when an external query specifies the coordinate to be classified only after sampling stops. Under coordinatewise error control, the optimal symmetric average sample size is exactly $m$ times the one-coordinate optimum. A change-of-measure argument gives the corresponding pointwise lower bound in terms of binary relative entropy.
Long COVID is a condition usually defined by persisting symptoms following infection by the SARS-CoV-2 virus beyond the acute phase of infection. The condition has a significant impact on healthcare systems, the economy, and the individuals living with it. Due to the diverse and extensive symptomatology of long COVID, symptom co-occurrence tracking can be used to capture patient experiences and improve understanding, diagnosis, and management. Here, we leverage the UK Office for National Statistics (ONS) COVID-19 Infection Survey (CIS). The CIS was run between April 2020 and March 2023. On February 3, 2021, the ONS launched a new CIS question evaluating self-reported symptom persistence of 23 long COVID symptoms post self-reported COVID-19 infection. We use Jaccard symptom-by symptom distance matrices, derived from the binary survey responses of the presence of each symptom. Three methods are then used to visualise symptom co-occurrence: heatmaps, non-metric multidimensional scaling, and agglomerative hierarchical clustering with complete linkage. We split our analysis into two parts, first looking at all long COVID survey responses (n = 207,319) and then looking at responses stratified by time-since-onset of long COVID up to 24 months (n = 30,224 at zero months-since-onset). We find higher symptom co-occurrence in prolonged long COVID. We also find clusters of neurological/systemic symptoms, gastrointestinal symptoms, and respiratory symptoms, with separability in symptom co-occurrence by organ system becoming less pronounced as time-since-onset of long COVID increases. Shortness of breath, weakness/tiredness, and muscle ache present as core symptoms long COVID. We demonstrate that the symptom experience of long COVID evolves from early stages to late stages with higher symptom burden and multi-systemic presentation. We also find three possible phenotypes of early long COVID.
Forecasting time series accurately is critical for applications with complex data ranging from energy systems to healthcare and finance. Among current state of the art models, generative latent variable models are increasingly implemented; yet principled generalisation guarantees for modern latent variable models remain limited. In particular, while Variational AutoEncoders are widely used for sequential data, their theoretical analysis is largely restricted to i.i.d. settings. In this work, we develop a PAC-Bayesian framework for latent variables models applied to time series. Building on reconstruction-based bounds, we extend PAC-Bayesian guarantees to Markovian latent structures, capturing temporal dependencies through a sequential generative process. These guarantees do not grow with the length of the trajectory. Our bounds depend on assumptions which are common in the literature; we provide an example framework where they would be verified to show that they are not as restrictive as they may seem.
Dynamical symbolic regression methods identify governing differential equations from noisy data, balancing interpretability and predictive accuracy. However, standard methods often produce expressions that violate known physical laws. To address this, we propose FluxDisco, a physics-informed framework tailored for flux-based, stoichiometric ODE systems. By leveraging a known stoichiometry, we reduce the expression search space and ensure physical adherence. Our framework adapts the Monte Carlo Graph Search algorithm for the unique challenges associated with joint flux discovery of stoichiometric systems. We evaluate our method across a range of physical and biological systems, demonstrating its ability to accurately recover governing dynamics through interpretable equations.
We develop a finite-horizon calibration method for raw-score martingale posteriors, with von Mises--Fisher models as the main worked example. Starting from the maximum likelihood estimator, predictive paths are generated by simulating future observations from the current fitted model and updating the natural parameter by unpreconditioned score increments. The main methodological step is to separate predictive simulation from covariance calibration. Raw-score increments have Fisher-information covariance, whereas Bernstein--von Mises calibration requires inverse-information covariance. We therefore apply a terminal linear correction based on a local information estimate. For more efficient implementation, we also introduce a hybrid version that replaces the omitted tail of the infinite predictive continuation by a Gaussian approximation with matching leading-order quadratic variation. We prove fixed-n convergence, a finite-horizon approximation bound for the Gaussian tail, and a Bernstein--von Mises limit for the hybrid post-corrected sampler under local regularity and consistent terminal calibration. Simulations show that tail correction reduces truncation-induced underdispersion, and an OSCAR ocean-current example illustrates local directional uncertainty summaries.
We present an approach for determining sample sizes required to detect underestimations of the expected shortfall with a prescribed power when applying the recently proposed e-backtesting procedure. We consider scenarios in which the value-at-risk at level $p$ is always estimated correctly, while the difference between the true expected shortfall and the value-at-risk is underestimated by a given factor $r$. We show that exploiting the structure of the backtest e-statistic proposed for backtesting the expected shortfall at level $p$ enables the derivation of approximate lower bounds for the required sample sizes by considering a sequence of independent and identically distributed Bernoulli random variables. We also discuss potential limitations of this approximation and compare the resulting sample size requirements with those obtained in practical applications using Monte Carlo simulations. Furthermore, we present generalizations of the e-backtesting procedure, in particular to risk measures which constitute Bayes pairs.
We develop a framework for simultaneous change-point inference of high-dimensional functional time series. The observations are modeled as temporally dependent vectors whose coordinates take values in possibly different separable Hilbert spaces, thereby covering a broad class of functional data. Heterogeneous mean changes may occur at coordinate-specific locations, and the contemporaneous dependence across coordinates is left unrestricted. Our procedure is based on coordinatewise cumulative-sum statistics and a residual block multiplier bootstrap that provides a common critical value for the global test and all coordinatewise decisions. We establish a nonasymptotic Gaussian approximation, quantitative nonasymptotic bounds for strong family-wise error control under arbitrary mixtures of changed and unchanged coordinates, covariance-adaptive detection guarantees, and simultaneous high-probability bounds for change-point localization. The bounds accommodate high-dimensional regimes in which the number of functional coordinates grows exponentially in a power of the sample size. We investigate finite-sample performance in simulations and illustrate the method using river discharge curves and high-frequency financial log returns.
We study how the response of a Bayesian posterior statistic to future observations changes as information accumulates. For a non-decreasing function $T$, define $Π_n^T=\E[T(Θ)\vert \mathcal F_n]$, where $Θ$ has an arbitrary prior and the observations come from a one-parameter exponential family. Conditioning on the same current value of $Π^T$, we show that the posterior statistic after additional observations is larger in convex order when the current posterior is based on fewer observations. We also prove preservation of convexity: the expected value of a convex function of the future posterior statistic is convex in the current posterior statistic. Together, these two properties provide structural tools for establishing time-monotonicity results in dynamic Bayesian decision and optimal stopping problems. If the exponential family contains an infinitely divisible distribution, the results extend to a continuous-time observation model through a family of Lévy processes.
We study regression under bounded responses in terms of excess mean squared error. When the comparator class is finite, this setting is known as model selection aggregation, and achieving minimax excess risk requires improper learning algorithms. Contrary to this, in the universal learning framework no improperness is needed, as simple empirical risk minimization achieves the best-possible exponential learning rate. Hence, the two frameworks suggest different optimal algorithmic principles. This poses the question of best-of-both-worlds guarantees: Are minimax and universal exponential rates achievable by the same algorithm? For finite hypothesis classes, we answer this question in the affirmative by showing that the $Q$-aggregation estimator - which is known to achieve minimax optimal tails - achieves exponential universal rates. A wide range of other estimators and algorithmic principles (ERM, sequential averaging, pruning, and star estimation) do not achieve both. For countably infinite hypothesis classes, we answer the question in the negative by showing that there is an inherent trade-off between achieving exponential universal and minimax uniform rates. This trade-off is exactly traced by combining optimal algorithms from each world using $Q$-aggregation. Besides these results, we prove several additional structural results about universal rates in learning with squared loss.
We consider density estimation under the relaxed local differential privacy condition that the privatized distributions are $α$-close in total variation distance. We show that adding independent noise with a convenient symmetrized Gamma distribution to each sensitive observation attains the $α$-TV-LDP. We prove that the deconvolution estimator of $r$-Sobolev smooth functions attains the pointwise rate $(nα)^{-\frac{2r-1}{2r}}$ up to log factors which is faster than $(nα^2)^{-\frac{2r-1}{2r+1}}$ under the classical $α$-LDP and closer to the nonprivate minimax rate $n^{-\frac{2r-1}{2r}}$. Next, we use a Goldenshluger-Lepski procedure to build a free of the smoothness adaptive procedure and show optimality of our rates in the convolution model of our privatisation scheme. We illustrate the benefits of this simple privacy mechanism by implementing a neural network estimator which does not need to add more noise in the optimization steps. Numerical results show significant improvement of the estimation rate over the Laplace and the private-SGD mechanisms.
Self-supervised learning relies on so-called data augmentations $φ(x)$ of unlabeled datapoints $x$ --- for example, masking random pixels in an image $x$ --- that should leave the label of $x$ invariant and are often used to learn a lower-complexity invariant subspace $\cal V$ for downstream tasks. In practice, such augmentations $\{ φ_l(x_i) \}$ are pooled together to learn $\cal V$, despite obvious inter-dependencies between different augmentations $φ_l(x), φ_k(x)$ of the same datapoint $x$. However, theoretical works on the subject typically consider procedures that avoid such dependencies, and are therefore limited to operate on smaller subsets of independent data. We show in this work that pooling augmentations together, despite inter-dependencies, is a better alternative than the baseline of partitioning the data into subsets of independent data. More precisely, in the context of estimating $\cal V$, the statistical estimation error bounds for pooling are never worse than the partitioning baseline, and in some cases --- such as masking or noise injection-based augmentations over a shallow neural network --- naive pooling leads to faster rates in terms of the number of augmentations. The benefits of pooling are particularly prominent when the correlations between different augmentations $φ_l(x), φ_k(x)$ have mild effects on estimation or help decrease the estimation variance. The analysis, therefore, yields new insights into the success of pooling augmented samples in self-supervised pre-training, and provides an intuition behind the practical preference towards using many augmentations.
Bridging the gap between animal and human experiments remains a major challenge in translational medicine, particularly in early drug development. Progress is constrained by financial cost, the difficulty of integrating heterogeneous in vitro and in vivo data, and the desire to reduce the use of animal testing balanced against minimising the risk to human participants. We present a statistical machine learning framework using multi-fidelity Gaussian processes, in which animal studies are considered as lower fidelity but informative approximations to human experiments. This allows cross-species similarities and nonlinear exposure-response relationships to be learned simultaneously, enabling principled extrapolation between species while quantifying uncertainty. By leveraging information from multiple experimental fidelities, our method improves estimation of clinically relevant quantities of interest and supports the replacement, reduction, and refinement of in vivo testing. We first illustrate this approach in a simulated scenario, before validating it on real clinical data. We simulate data for drug-induced QT-interval prolongation, a key cardiac safety assessment required for regulatory approval. This framework provides a probabilistic surrogate capable of integrating in vitro pharmacology, animal experiments and human data within a unified statistical model. Crucially, it achieves this at no additional experimental cost while also enabling transfer learning across compounds. As a result, predictions and uncertainty quantification for new drugs can be generated from in vitro findings alone, providing additional efficiency gains and accelerating decision making. For validation, we use a clinical dataset measuring change in heart rate under autonomic blockade, which represents some of the challenges commonly found in multi-species datasets.
Problem definition: Fixed parameters and inputs in a stochastic simulator induce a distribution over complete trajectories, not one trajectory. Calibration must assess this distribution, including variability and temporal dependence, against observations. Yet stochastic car-following models are commonly calibrated with trajectory-error objectives inherited from deterministic modelling. Methodology/results: We establish a scoring-rule theory of stochastic calibration. Strict propriety requires the data-generating distribution to uniquely minimise expected score. MRMean-I, the average run-wise error, drives separable stochastic spread to zero; MRMean-II, the error of the ensemble-mean trajectory, cannot identify a parameter that changes only spread; and MRMin, the error of the closest simulated run, has a population target that changes with ensemble size. These results are confirmed for stochastic Intelligent Driver Model extensions with additive acceleration noise and random desired headway. We recommend exact maximum likelihood when the correct transition density is available; otherwise, an unbiased simulation-based estimator of a strictly proper score. The energy score meets this requirement and gives the best held-out distributional prediction among the evaluated simulation-based objectives, although both models retain too-narrow bands and miss persistent disturbances. Implications:Strict propriety separates a valid calibration target from parameter identifiability and model adequacy. The theory applies to vector-valued outputs from stochastic transportation simulators; the car-following experiments illustrate its scope.
Multivariate correlated outcomes occur across disciplines, including ecology, social sciences, and psychometrics. This paper focuses on clustering these outcomes across observational units, specifically, finding groups of units with the same ``outcome profile". Our motivation comes from bioregionalization in ecology, which aims to cluster sites into bioregions with the same species profiles, where site membership can depend on environmental or habitat covariates. To accomplish this, we propose finite mixtures of generalized estimating equations (MixGEE). Unlike existing approaches to model-based bioregionalization, MixGEE partitions sites into regions while accounting for between-species correlations through a region-specific working correlation structure. Thus, each region is characterized by a marginal species mean vector and a between-species correlation matrix. Unlike likelihood-based finite mixture models, MixGEE does not require a full joint distribution of the multivariate outcomes. Instead, we construct a pseudo-posterior probability for region membership motivated by the large-sample distribution of the estimating equation. This leads to an iterative algorithm alternating between updating these probabilities and solving weighted estimating equations. We determine the number of regions using cross-validation based on predictive performance for held-out sites and species components, and use a clustered Dirichlet random-weight bootstrap for uncertainty quantification. Simulations demonstrate reliable estimation and inference under various correlation structures and more stable selection of the number of groups than methods that ignore dependence. Applying MixGEE to presence--absence records of fish species around the Kerguelen Plateau reveals three distinct fish assemblage profiles with heterogeneous occurrence and within-site correlation patterns.
Longitudinal clinical data are increasingly available, offering opportunities to study treatment effects over time but also raising challenges related to time-varying confounding and evolving treatment decisions. We investigate how longitudinal treatment dynamics affect causal effect estimation when baseline and longitudinal methods target different causal estimands. Using a structural causal model (SCM), we simulate data with time-varying confounding, binary treatment, and an absorbing binary outcome under 9 scenarios combining functional complexity and treatment persistence. We compare three baseline and two longitudinal estimators against Monte Carlo ground-truth risk ratios (RRs) under sustained treatment regimes. Results show that treatment persistence is the main driver of estimator behaviour. High persistence reduces the discrepancy between baseline and sustained-regime estimands, making baseline estimators closer to the sustained-regime ground truth (average relative deviation of baseline IPTW/TMLE decreasing from 61% under low persistence to 10% under high persistence), whereas low persistence induces practical positivity challenges and increases the variability of longitudinal estimators, with empirical 95% interval widths increasing from 0.22 to 0.46 for longitudinal IPTW and from 0.20 to 0.36 for LTMLE when moving from high to low persistence. These findings emphasise that estimator performance should be interpreted jointly with the target intervention and the treatment process generating the observed data.
Testing whether two independent samples arise from the same underlying distribution is a fundamental statistical problem. We propose a new class of two-sample distribution tests based on a family of probability metrics $Δ_{n,p}$, constructed from repeatedly integrated quantile functions. On their respective domains, these metrics are proved to be genuine distributional distances. The case $n=1$ recovers the $p$-Wasserstein distance, which requires finite $p$-th moments; for $n\geq2$, the proposed metrics are well defined and require only finite first moments. The asymptotic properties of the plug-in statistic are established, including strong consistency and limiting distributions under the null and fixed alternatives. A permutation calibration for finite-sample inference is also proposed. We further derive an asymptotic power function under local alternatives. Finally, the finite-sample performance of the proposed tests is examined through simulation studies, and their reduced sensitivity to extreme upper-tail observations is illustrated through a real data application.
The statistical analysis of heavy-tailed data has received considerable attention because extreme observations frequently arise in many practical applications. The Pareto type-I distribution is a fundamental heavy-tailed model used in economics, finance, actuarial science, insurance, reliability, and extreme value analysis. In this paper, we propose novel goodness-of-fit tests for the Pareto distribution using a mean residual life characterization. The test statistic is constructed using U-statistic theory, and its asymptotic behaviour is established under both the null and alternative hypotheses. Its finite-sample performance is evaluated through Monte Carlo simulations using maximum-likelihood and method-of-moments estimation and compared with existing tests. The results show that the proposed test controls the nominal significance level and performs competitively in terms of power across a broad range of alternatives. Finally, the proposed methodology is illustrated using the Danish fire insurance loss and pollution datasets.
We propose a Bayesian model-based approach for clustering functional data with an unknown number of clusters and temporally correlated observations. Cluster-specific mean functions are represented using B-spline basis expansions, while within-curve dependence is modeled through an Ornstein--Uhlenbeck covariance structure. A truncated Dirichlet process mixture is used to infer the effective number of clusters, and a variational EM algorithm is developed for efficient posterior approximation. Simulation studies show that the proposed method performs well under both correctly specified and misspecified mean-function settings and compares favorably with several existing functional clustering methods. Comparisons with MCMC indicate that the variational approximation yields consistent clustering and parameter estimates but at substantially lower computational cost. An application to Canadian daily temperature curves further demonstrates the practical usefulness of the method in identifying interpretable functional clusters while accounting for temporal dependence.
While diffusion-based methods have recently emerged as effective tools for probing the intrinsic geometry of high-dimensional data, their statistical difficulty remains largely unexplored. We study estimation of the finite-scale population functional underlying FLIPD (Kamkari et al., 2024; arXiv:2406.03537), a diffusion-based local intrinsic dimension (LID) quantity defined through the logarithmic scale derivative of a Gaussian-smoothed density. Intuitively, Gaussian smoothing turns local dimension into a scale law: near a $d$-dimensional manifold, the kernel mass grows like $σ^d$, so differentiating with respect to the noise scale reveals the intrinsic exponent. Under a regular manifold model, we show uniformly over the model class that the finite-scale field differs from the manifold dimension $d$ by at most $O(σ^2)$. We then establish a minimax lower bound of order $(nσ^d)^{-1}$ for estimating this finite-scale field from $n$ observations, for $n^{-1/(2α+d)}\lesssimσ\leσ_0$. At the smallest scale covered by our lower-bound construction, the bound becomes the nonparametric rate $n^{-2α/(2α+d)}$.
This study proposes two novel bivariate distributions for jointly modeling temperature and rainfall by integrating Kumaraswamy-Teissier marginals with Clayton and Gumbel copula structures. To capture a wide range of dependence patterns, including both positive and negative associations, rotated copula variants (90°, 180°, and 270°) are incorporated along with their corresponding tail dependence characteristics. Model parameters are estimated using maximum likelihood and the inference functions for margins (IFM) approach, and their finite-sample performance is assessed through a comprehensive Monte Carlo simulation study. The proposed models are applied to monthly gridded temperature and rainfall data from the Northwest Himalayas, a region characterized by complex hydro-climatic variability. Comparative analysis demonstrates that the proposed framework outperforms several existing bivariate models and effectively captures lower-tail, upper-tail, and asymmetric dependence structures across summer and winter seasons. Based on the selected best-fitting copula models, univariate, joint, and conditional return periods are derived to quantify the risk of compound extremes. The results highlight the capability of the proposed approach to provide a more realistic representation of hydro-climatic dependence and offer a robust framework for assessing the risk of extreme temperature and rainfall events in mountainous regions.
This paper introduces a quantile-based Kumaraswamy-Teissier autoregressive moving average (KTARMA) model for positive-valued time series. Leveraging the flexibility of the extended Kumaraswamy-Teissier distribution within an observation-driven framework, the random component of the distribution is conditioned on the historical process and time-varying covariates, and is parameterized explicitly via its $ρ$-th conditional quantile, where $ρ\in (0,1)$. To capture temporal dependence, the systematic component maps an ARMA-type structure to this conditional quantile via an appropriate link function. For inference, we implement a conditional maximum likelihood framework and derive explicit analytical expressions for the resulting score vector and conditional information matrix, followed by the development of model diagnostic and forecasting procedures. The finite-sample performance of the developed estimators is evaluated through a Monte Carlo simulation study across various parameter configurations and quantile levels. Finally, the practical utility of the study is demonstrated by modeling monthly rainfall data over the Northwest Himalayas (2001-2025), where 525 grids are grouped into four homogeneous zones using a Self-Organizing Map and relevant atmospheric variables and large-scale climate indices are incorporated as predictive regressors. Out-of-sample forecasting evaluations reveal that the KTARMA model delivers highly competitive predictive performance, achieving consistently lower mean squared errors across all identified zones compared to KARMA and $β$ARMA models.
In medical research, it is often of interest to evaluate the predictive performance of a biomarker. Statistical approaches based on the Receiver Operating Characteristic (ROC) curve and its summary measures, such as the area under the curve (AUC) and the Youden index, are widely used to evaluate the prognostic performance of these biomarkers. In time-to-event studies, ROC analysis poses additional challenges due to change in disease status over time and the presence of censored individuals. To address these issues, time-dependent ROC curves were introduced. In this paper, we propose a non-parametric estimator of the time-dependent two-way partial AUC for right-censored data. We also discuss the partial Youden index and the associated optimal biomarker cutoff estimator for the right-censored data. We conduct an extensive simulation study to investigate the finite sample performance of the proposed estimators. The simulation study indicates that the proposed non-parametric estimators efficiently account for right censoring. Finally, we illustrate the proposed methods using two real data sets, one from the Primary Biliary Cirrhosis study and the other from the Molecular Taxonomy of Breast Cancer International Consortium trial.