A common technique to quantitatively estimate P speciation in soil samples is to apply linear combination fitting (LCF) to normalized P K-edge X-ray absorption near-edge structure (XANES) spectra. Despite the rapid growth of such applications, the uncertainties of the fitted weights are still poorly known. Further, there are few reports to what extent the LCF standards represent unique end-members. Here, the co-variance between 34 standards was determined and their significance for LCF was discussed. We present a probabilistic approach for refining the calculation of LCF weights based on Latin hypercube sampling of normalized XANES spectra, where the contributions of energy calibration and normalization to fit uncertainty were considered. Many of the LCF standards, particularly within the same standard groups, were strongly correlated. This supports an approach in which the LCF standards are grouped. Moreover, adsorbed phytates and monetite were well described by other standards, which puts into question their use as end-members in LCF. Use of the probabilistic method resulted in uncertainties ranging from 2 to 11 percentage units. Uncertainties in the calibrated energy were important for the LCF weights, particularly for organic P, which changed with up to 2.7 percentage units per 0.01 eV error in energy. These results highlight the necessity of careful energy calibration and the use of frequent calibration checks. The probabilistic approach, in which at least 100 spectral variants are analyzed, improves our ability to identify the most likely P compounds present in a soil sample, and a procedure for this is suggested in the paper.