A developer gave me this suggestion,
= = =
He is attempting to use ProbabilityDistribution to model a mix of
non-overlapping distributions.
His distribution is equivalent to
spd = SplicedDistribution[{0.03673118299607324`, 0.8606546277105262`,
0.10261418929340049`}, {0, 151, 12846,
Infinity}, {PowerDistribution[1/151, 1.8923],
BetaPrimeDistribution[1, 2701/10000, 5000/3333,
495.9127480196388`], ParetoDistribution[12846, 3253/2500]}];
Now the sampling is faster and is more precise:
In[317]:= AbsoluteTiming[mpf = RandomVariate[spd, 10^6];]
Out[317]= {1.190119, Null}
= = =