1
$\begingroup$

I have a custom probability distribution function and I want to generate numbers according to this PDF. Interestingly enough, this works without trouble on my MacBook with Mathematica 10.0, but when I try to run the exact same notebook on a Linux machine having later Mathematica versions, the execution takes orders of magnitude longer and even tends to crash the kernel.

An example of a PDF that shows up in my computations is the following:

plot of probability distribution function

It blows up for t=0,1,2. However, on my Macbook this does not cause any problems. If I cut off the pdf at a certain maximum value, the linux implementation also seems to work, so I suspect these singularities to cause the issue. What could be the reason for this discrepancy between different operating systems?

I'm running Mathematica 10.0.2 on the MacBook, and tried several different versions on the Linux machine, including 9.0, 10.3 and 11.0, so I don't expect this to cause the problem. The problem was confirmed to be due to the version difference, not the operating system difference (see comments).

The reason I didn't provide code initially is that the pdf is generated through a complicated process that will not illustrate the problem. However, to provide a little more detail, here is an example of a typical pdf:

pdf = 1. (1/(2 \[Pi]) \[Sqrt](0.0332814 (12. Cos[
            0.5 (1.45946 - 2.91891 (-1 + t))] - 
          1. Sqrt[-80. + 
            144. Cos[0.5 (1.45946 - 2.91891 (-1 + t))]^2])^2 + 
       0.015625 (17.5135 Sin[0.5 (1.45946 - 2.91891 (-1 + t))] - (
          210.162 Cos[0.5 (1.45946 - 2.91891 (-1 + t))] Sin[
            0.5 (1.45946 - 2.91891 (-1 + t))])/
          Sqrt[-80. + 
           144. Cos[0.5 (1.45946 - 2.91891 (-1 + t))]^2])^2) UnitStep[
     2 - t] UnitStep[-1. + t] + (1/(
   2 \[Pi]))\[Sqrt](0.0332814 (12. Cos[0.5 (-1.45946 + 2.91891 t)] + 
          Sqrt[-80. + 144. Cos[0.5 (-1.45946 + 2.91891 t)]^2])^2 + 
       0.015625 (-17.5135 Sin[0.5 (-1.45946 + 2.91891 t)] - (
          210.162 Cos[0.5 (-1.45946 + 2.91891 t)] Sin[
            0.5 (-1.45946 + 2.91891 t)])/
          Sqrt[-80. + 
           144. Cos[0.5 (-1.45946 + 2.91891 t)]^2])^2) UnitStep[
     1 - t] UnitStep[t])

The pdfs are always normalized to 1. Sometimes it acquires a small imaginary part due to numerical problems, but I define the ProbabilityDistribution as

cdist = ProbabilityDistribution[Re[pdf], {t, 0, 2}];
clist = RandomVariate[cdist, n];

so the imaginary part is taken care of. The second line is where the problem starts. As I said, cutting off pdf at some maximum value seems to help somewhat, but it's still orders of magnitude slower than on macOS.

I'm currently considering to just generate random numbers according to a uniform distribution, and subsequently mapping them to this distribution. However, I'm still interested in finding out what may cause this problem.

$\endgroup$
  • 1
    $\begingroup$ Is your version of Mathematica the same on both computers? Have you tried running it in the cloud, since the cloud runs the Linux version? $\endgroup$ – Carl Lange May 8 at 18:42
  • $\begingroup$ Sorry, forgot to mention that. I have Mathematica 10.0.2 on my Macbook, and several different versions on the linux machine, including 10.3 and 11.0. However, all the versions on the linux machine show the same problem. I have not tried running it in the cloud - let me give that a shot too. $\endgroup$ – JorenB May 8 at 18:44
  • 6
    $\begingroup$ Can you please provide your code? Otherwise it's difficult to make any statements. $\endgroup$ – Roman May 8 at 18:56
  • $\begingroup$ Your pdf looks like two separate beta distributions but does the pdf integrate to 1 over 0 to 2? It looks like it could integrate something larger than 1 (like maybe 2). $\endgroup$ – JimB May 9 at 0:55
  • $\begingroup$ I added some code, but I'm not sure this is actually helpful. And @JimB, the pdfs are always normalized to 1. $\endgroup$ – JorenB May 9 at 9:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.