<p>Two critical challenges, among others, in unsupervised damage detection of bridges are the lack of uncertainty quantification for damage indices, such as the widely used Mahalanobis squared distance (MSD), and issues related to limited data, which arise when the number of damage-sensitive features is close to or exceeds the number of available time series samples. This paper introduces a novel method that addresses both challenges by computing the variability of the MSD as a damage index through explicit propagation of uncertainties in autoregressive (AR) coefficients and their associated covariance matrices. The method also addresses the limited data problem by taking advantage of the least squares estimation of AR coefficients to derive estimator covariance matrices that characterize the uncertainty of the features. The method is validated using experimental data from a full-scale, densely instrumented bridge in both undamaged and damaged conditions. The results show that the proposed method allows for a rigorous non-parametric statistical analysis of damage index distributions, achieves excellent damage detection performance, and supplements the damage detection results with a quantifiable measure of uncertainty. As such, the proposed method enables more robust and efficient data-driven decision-making for existing bridges.</p>