One can use many simultaneous random chains. It considerably speeds up the simulation. Also we can implement own distribution (see Defining Distributional Generators [here](http://reference.wolfram.com/mathematica/tutorial/RandomNumberGeneration.html)) ClearAll[lifetime] lifetime /: Random`DistributionVector[lifetime[size_, chains_ : 200], n_, prec_] := Module[{m, t, res, pos}, Flatten[Table[ m = ConstantArray[-1, size chains]; t = 0; res = ConstantArray[0, chains]; While[Min[res] == 0, t++; pos = Range@chains + chains RandomInteger[size - 1, chains]; m[[pos]] *= -1; res += UnitStep[-res] UnitStep[Total@Partition[m, chains] - size] t; ]; res, {⌈n/chains⌉}]][[1 ;; n]] ]; Here `size` is the size of the list, `chains` is the number of simultaneous chains (the default value 200 is optimal) and `n` is the number of generated numbers. Now we can use it as a regular distribution with `RandomVariate` N@Mean@RandomVariate[lifetime[6], 100000] > 83.4748 Histogram[RandomVariate[lifetime[6], 100000], {0, 150, 2}, AxesOrigin -> 0] ![enter image description here][1] The speed is several times faster then heropup's uncompiled version because I use [packed arrays](https://mathematica.stackexchange.com/q/3496/4678). One can put the code above to `Compile`: lifetimeCFun = Compile[{{size, _Integer}, {chains, _Integer}, {n, _Integer}}, Module[{m, t, res, pos}, Flatten[Table[ m = ConstantArray[-1, size chains]; t = 0; res = ConstantArray[0, chains]; While[Min[res] == 0, t++; pos = Range@chains + chains RandomInteger[size - 1, chains]; m[[pos]] *= -1; res += UnitStep[-res] UnitStep[Total@Partition[m, chains] - size] t; ]; res, {⌈n/chains⌉}]][[1 ;; n]] ], CompilationTarget -> "C", RuntimeOptions -> "Speed"]; lifetimeC /: Random`DistributionVector[lifetimeC[size_, chains_ : 200], n_, prec_] := lifetimeCFun[size, chains, n]; The speed of my compiled solution is comparable to heropup's compiled solution. [1]: https://i.sstatic.net/WVMHF.png