A Bayesian fitting algorithm was combined with analytical approximations of the Bloch-McConnell (BM) equations with the aim to considerably reduce processing time. The accuracy of the algorithm was assessed with simulated data and data from phantom experiments and compared to fit results obtained with the numerical solution of the BM equations. Continuous-wave and pulsed saturation was considered. The results showed agreement between estimates and ground truth as well as between the approximate analytical and numerical model implementations of the Bayesian algorithm. A considerable reduction of processing time was achieved.