Paleoliquefaction features can be used to estimate lower bounds on the magnitude and ground motion associated with the earthquake that caused the liquefaction feature. The engineering back-analysis of paleoliquefaction features is usually conducted using state of the practice liquefaction-triggering analysis methodologies. Recent studies have shown that these methodologies are associated with variable probabilities of liquefaction depending on the soil parameters. This would imply that estimates of magnitude and ground motion intensity obtained from these methodologies would not be consistent for all soil sites. Moreover, these estimates could be unconservative. In this paper, the use of a probabilistic methodology for the backanalysis of paleoliquefaction features is proposed. The proposed methodology permits the incorporation of model and parameter uncertainty into the analysis and results in more robust estimates of past magnitude and a measure of the uncertainty associated with these predictions. Previously published paleoliquefaction data are used to demonstrate the applicability of the proposed method. Magnitude estimates obtained with the proposed method do not differ significantly from those obtained using deterministic methodologies, but the proposed methodology permits a quantification of the uncertainty associated with magnitude estimates.