Principal component analysis (PCA) is a well-established technique in paleomagnetism and provides a means to estimate magnetic remanence directions from univectorial segments of stepwise demagnetization data. Derived directions constrain past geomagnetic field behavior and form the foundation of chronological and tectonic reconstructions. PCA of isolated remanence segments relies on estimates of the segment mean and covariance matrix, which can carry large uncertainties given the relatively small number of demagnetization data points used to characterize individual specimens. Traditional PCA does not, however, lend itself to quantification of these uncertainties, and inferences drawn from paleomagnetic reconstructions suffer from an inability to propagate uncertainties from individual specimens to higher levels, such as in calculations of paleomagnetic site mean directions and pole positions. In this study, we employ a probabilistic reformulation of PCA that represents the unknowns involved in the data fitting process as probability density functions. Such probability density functions represent our state of knowledge about the unknowns in the fitting process and provide a tractable framework with which to rigorously quantify uncertainties associated with remanence directions estimated from demagnetization data. These uncertainties can be propagated readily through each step of a paleomagnetic reconstruction to enable quantification of uncertainties for all stages of the data interpretation sequence, removing the need for arbitrary selection/rejection criteria at the specimen level. Rigorous uncertainty determination helps to protect against spurious inferences being drawn from uncertain data.