Problem with RandomVariate in Mathematica 10 (and 9) Is it a bug?

Questions
,
1

Hello,

I have a problema with RandomVariate in Mathematica 10 and 9, in the sense that
for a specific distribution constructed by myself, the RandomVariate sample does not follow
the distribution the sample comes from.

I would like to discard that it be a specific problem of my machine or my installation.

Thus I would be grateful if anyone of you could run the attached notebook with Mathematica 10
and send it back to me. I have scanned the notebook and it is free from viruses according
to my antivirus.

Thanks in advance,
Arturo Ramos

2

Here you go! Same result as you, it appears.

In[1]:= $Version

Out[1]= "10.0 for Mac OS X x86 (64-bit) (June 29, 2014)"
3

You can always make a function that is hard to numerically integrate. In such a case, it’s sometimes hard to know whether you should consider it to be a bug or not. Similarly, it is always possible to make a PDF that doesn’t get sampled well enough at the tails with RandomVariate for numerical reasons. And it’s hard to know whether it should be considered a bug or not. There isn’t a always a clear yes or no answer.

It looks like it’s not RandomVariate isn’t sampling well at some very extreme tails. If you ask it to make the ProbabilityDistribution of f[Exp[x]] Exp[x] instead of f to make the tails more reasonable, everything seems to work fine.

There’s a couple of things you can try here to help this to work. The first I would try is to rewrite as a Piecewise function. The function f seems numerically complicated and this should help simplify things down a bit. Can you reproduce the problem with simpler values?

The function f doesn’t integrate easily with NIntegrate. This might be better fixed by using a Piecewise function If something is difficult to numerically integrate (say with NIntegrate) it stands to reason that you’d have a hard time making a very good random number generator from it.

4

Thanks a lot Jesse, you have been of great help!.

Best,
Arturo

5

Thanks a lot Sean,

You have given me a pair of good ideas and
now I could write my code more efficiently adapted to the mathematics in play.

In fact, with the Exp transformation anything works well and the Piecewise command
is something I will try.

Thanks again!

Best,
Arturo

6

Dear Sean,

Now I have tried the sampling of f[Exp[y]]Exp[y] and works very well.

I have also tried using Piecewise[] instead of the HeavisideTheta[]
and the result is that the computations are quite slower with Piecewise[]
and it does not solve the initial problem. So I keep the first idea for the future!

Thanks again and all the best.

7

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}

= = =

8

Thanks Bruce,

It is another good idea to take into account.