Perpendicular magnetic tunnel junction (pMTJ)-based true-random number generators (RNGs) can consume orders of magnitude less energy per bit than CMOS pseudo-RNGs. Here, we numerically investigate with a macrospin Landau–Lifshitz-Gilbert equation solver the use of pMTJs driven by spin–orbit torque to directly sample numbers from arbitrary probability distributions with the help of a tunable probability tree. The tree operates by dynamically biasing sequences of pMTJ relaxation events, called ‘coinflips’, via an additional applied spin-transfer-torque current. Specifically, using a single, ideal pMTJ device we successfully draw integer samples on the interval [0, 255] from an exponential distribution based on p-value distribution analysis. In order to investigate device-to-device variations, the thermal stability of the pMTJs are varied based on manufactured device data. It is found that while repeatedly using a varied device inhibits ability to recover the probability distribution, the device variations average out when considering the entire set of devices as a ‘bucket’ to agnostically draw random numbers from. Further, it is noted that the device variations most significantly impact the highest level of the probability tree, with diminishing errors at lower levels. The devices are then used to draw both uniformly and exponentially distributed numbers for the Monte Carlo computation of a problem from particle transport, showing excellent data fit with the analytical solution. Finally, the devices are benchmarked against CMOS and memristor RNGs, showing faster bit generation and significantly lower energy use.