We consider the probability distributions of the subsystem (staggered) magnetization in ordered and disordered models of quantum magnets in D dimensions. We focus on Heisenberg antiferromagnets and long-range transverse-field Ising models as particular examples. By employing a range of self-consistent time-dependent mean-field approximations in conjunction with Holstein-Primakoff, Dyson-Maleev, Schwinger boson and modified spin-wave theory representations we obtain results in thermal equilibrium as well as during non-equilibrium evolution after quantum quenches. To extract probability distributions we derive a simple formula for the characteristic function of generic quadratic observables in any Gaussian theory of bosons.
