The standard quantum mechanical expressions for the singlet and triplet survival probabilities and product yields of a radical pair recombination reaction involve a trace over the states in a combined electronic and nuclear spin Hilbert space. If this trace is evaluated deterministically, by performing a separate time-dependent wavepacket calculation for each initial state in the Hilbert space, the computational effort scales as O(Z2log⁡Z), where Z is the total number of nuclear spin states. Here we show that the trace can also be evaluated stochastically, by exploiting the properties of spin coherent states. This results in a computational effort of O(MZlog⁡Z), where M is the number of Monte Carlo samples needed for convergence. Example calculations on a strongly coupled radical pair with Z>106 show that the singlet yield can be converged to graphical accuracy using just M=200 samples, resulting in a speed up by a factor of >5000 over a standard deterministic calculation. We expect that this factor will greatly facilitate future quantum mechanical simulations of a wide variety of radical pairs of interest in chemistry and biology.