We analyse the probability distribution functions of the energy and magnetisation of the two-dimensional Blume–Capel (BC) and Baxter–Wu (BW) models with spin values
S
∈
{
1
/
2
,
1
,
3
/
2
}
in the presence of a crystal field Δ. By employing extensive single-spin flip Monte Carlo simulations and a recently developed method of studying the zeros of the energy probability distribution we are able to probe, with a good numerical accuracy, several critical characteristics of the transitions. Additionally, the universal aspects of these transitions are scrutinised by computing the corresponding probability distribution functions. The energy distribution has been underutilised in the literature when compared to that of the magnetisation. Somewhat surprisingly, however, the former appears to be more robust in characterising the universality class for both models upon varying the crystal field Δ than the latter. Finally, our analysis suggests that in contrast to the BC ferromagnet, the BW model appears to suffer from strong finite-size effects, especially upon increasing Δ and S, that obscure the application of traditional finite-size scaling approaches.