During magnetic resonance (MR)-guided thermal therapies, water proton resonance frequency shift (PRFS)-based MR temperature imaging can quantitatively monitor tissue temperature changes. It is widely known that the PRFS technique is easily perturbed by tissue motion, tissue susceptibility changes, magnetic field drift, and modality-dependent applicator-induced artifacts. Here, a referenceless Gaussian process modeling (GPM)-based estimation of the PRFS is investigated as a methodology to mitigate unwanted background field changes. The GPM offers a complementary trade-off between data fitting and smoothing and allows prior information to be used. The end result being the GPM provides a full probabilistic prediction and an estimate of the uncertainty. GPM was employed to estimate the covariance between the spatial position and MR phase measurements. The mean and variance provided by the statistical model extrapolated background phase values from nonheated neighboring voxels used to train the model. MR phase predictions in the heating ROI are computed using the spatial coordinates as the test input. The method is demonstrated in ex vivo rabbit liver tissue during focused ultrasound heating with manually introduced perturbations (n = 6) and in vivo during laser-induced interstitial thermal therapy to treat the human brain (n = 1) and liver (n = 1). Temperature maps estimated using the GPM referenceless method demonstrated a RMS error of <0.8°C with artifact-induced reference-based MR thermometry during ex vivo heating using focused ultrasound. Nonheated surrounding areas were <0.5°C from the artifact-free MR measurements. The GPM referenceless MR temperature values and thermally damaged regions were within the 95% confidence interval during in vivo laser ablations. A new approach to estimation for referenceless PRFS temperature imaging is introduced that allows for an accurate probabilistic extrapolation of the background phase. The technique demonstrated reliable temperature estimates in the presence of the background phase changes and was demonstrated useful in the in vivo brain and liver ablation scenarios presented.