4
$\begingroup$

A number of users have discussed the speed of Random number generation in Mathematica.

The Binomial and Multinomial random number generators in Mathematica are fast if multiple draws are needed for the same distribution parameters. For example, generating 500000 draws from the Binomial distribution is very quick

In[30]:= AbsoluteTiming[
 RandomVariate[BinomialDistribution[100, 0.6], 500000];]

Out[30]= {0.017365, Null}

However, the speed is slow compared to that in R and Julia when the parameters change across draws, as may be required when performing certain Monte Carlo simulations.

For example, if we have a vector nvec that contains the number of trials for each draw and a vector pvec that contains the corresponding probabilities of success.

nvec = RandomInteger[{5, 125}, 500000];

pvec = RandomReal[{0, 1}, 500000];

Then we have

In[28]:= AbsoluteTiming[
 Mean[Table[
    RandomVariate[BinomialDistribution[nvec[[i]], pvec[[i]]]], {i, 1, 
     Length@nvec}]] // N
 ]

Out[28]= {36.2144, 32.5283}

This hit in speed most probably stems from how these are implemented internally in Mathematica.

Are there alternate methods that are fast for the case when the parameter distribution changes across the draws?

$\endgroup$
  • $\begingroup$ I thought using ParameterMixtureDistribution might help. It didn't. dist = ParameterMixtureDistribution[BinomialDistribution[n, p], {n \[Distributed] DiscreteUniformDistribution[{5, 125}], p \[Distributed] UniformDistribution[{0, 1}]}] is just as slow. $\endgroup$ – Chris K Apr 24 at 12:22
5
$\begingroup$

Fast random number generators for the Binomial are described in

Kachitvichyanukul, V.; Schmeiser, B.W. "Binomial random variate generation." Comm.ACM 31 (1988), no .2, 216 - 222.

A Mathematica implementation based on the above paper involves the following three functions.

When the number of trials is small, the Geometric method from

Devroye. L. "Generating the maximum of independent identically distributed random variables." Computers and Mathematics with Applications 6, 1960, 305-315. can be used.

ablRanBinGeom = Compile[{{n, _Integer}, {p, _Real}}, 
   Module[
        {y = 0, x = 0, comp = 0, er, v, scale = 0.0},
        If[p >= 1.0, Return[n]];
        If[p <= 0.5, comp = 0; scale = -1.0/Internal`Log1p[-p], comp = 1;
      scale = -1.0/Log[p]];
        While[True, 
            er = -Log[RandomReal[]];
            v = er*scale;
            If[v > n, Break[]];
            y = y + Ceiling[v];
            If[y > n, Break[]];
            x = x + 1;
        ];
        If[comp == 1, n - x, x]
        ],
   CompilationTarget -> "C", RuntimeAttributes -> {Listable}
   ];

For larger n, we can use

    ablRanBinBtpe = Compile[{{n, _Integer}, {p, _Real}},
   Module[
        {comp = 0, r, q, nrq, fM, M, Mi, p1, xM, xL, xR, c, a, lamL, 
     lamR, p2, p3, p4,
        y, u, v, x, k, S, F, t, A, x1, f1, z, w, rho},

        If[p >= 1.0, Return[n]];

        If[p <= 0.5, 
            comp = 0; r = p; q = 1.0 - p, 
            comp = 1; r = 1.0 - p; q = p
        ];

        nrq = n*r*q;
        fM = (n + 1)*r;
        M = Floor[fM];
        Mi = Floor[M];
        p1 = Floor[2.195*Sqrt[nrq] - 4.6*q] + 0.5;
        xM = M + 0.5;
        xL = xM - p1;
        xR = xM + p1;
        c = 0.134 + 20.5/(15.3 + M);
        a = (fM - xL)/(fM - xL*r);
        lamL = a*(1.0 + 0.5*a);
        a = (xR - fM)/(xR*q);
        lamR = a*(1.0 + 0.5*a);
        p2 = p1*(1.0 + 2.0*c);
        p3 = p2 + c/lamL;
        p4 = p3 + c/lamR;

        y = 0;

        While[True, (* Step 1 *)

            u = p4*RandomReal[];
            v = RandomReal[];
            Which[
                u <=  p1, 
                y = Floor[xM - p1*v + u];
                Break[],

                u <= p2, (* Step 2 *)
                x = xL + (u - p1)/c;
                v = v*c + 1.0 - Abs[M - x + 0.5]/p1;
                If[v > 1, Continue[]];
                y = Floor[x],

                u <= p3 ,(* Step 3 *)
                y = Floor[xL + Log[v]/lamL];
                If[y < 0, Continue[]];
                v = v*(u - p2)*lamL,

                True, (* Step 4 *)

                y = Floor[xR - Log[v]/lamR];
                If[y > n, Continue[]];
                v = v*(u - p3)*lamR
            ];
            A = Log[v]; 
     If[A > (LogGamma[Mi + 1] + LogGamma[n - Mi + 1] + 
         LogGamma[y + 1] + LogGamma[n - y + 1] + (y - Mi)*Log[r/q]), 
      Continue[]
      ];
     Break[];
     ];
        If[comp == 1, n - y, y]
    ], 
   CompilationTarget -> "C", 
   CompilationOptions -> {"InlineExternalDefinitions" -> True}, 
   RuntimeAttributes -> {Listable}
   ];

We can automatically choose the appropriate method using the function

ablRanBinomial = Compile[{{n, _Integer}, {p, _Real}},
   Module[
        {q = 1.0 - p},
        If[n*Min[p, q] > 20, ablRanBinBtpe[n, p], ablRanBinGeom[n, p]]
    ],
   CompilationTarget -> "C", 
   CompilationOptions -> {"InlineExternalDefinitions" -> True}, 
   RuntimeAttributes -> {Listable}
   ];

Here is the performance improvement.

In[32]:= AbsoluteTiming[
  Mean@Table[
    ablRanBinomial[nvec[[i]], pvec[[i]]], {i, 1, Length@nvec}]] // N

Out[32]= {0.413019, 32.5307}

This can be further improved using the Listability attribute for the function.

In[33]:= AbsoluteTiming[Mean@ablRanBinomial[nvec, pvec] // N]

Out[33]= {0.156881, 32.5337}

Once we have the function for Binomial draws, this can be used in a function for fast generation of multinomial draws with n being the number of trails and p being a vector of probabilities that sum to 1.

ablRanMultinomial=Compile[{{n, _Integer},{p, _Real, 1}},
Module[
    {k=Length[p], rp,i, km1,cn,pi,xi,x},
    rp=1.0;cn=n;
    i=0;
    km1=k-1;
    x=Table[0, {k}];
    While[i<km1 && cn >0,
        i+=1;
        pi=p[[i]];
        If[pi < rp,
            xi=ablRanBinomial[cn, pi/rp];
            x[[i]]=xi;
            cn-=xi;
            rp-=pi;,
            x[[i]]=cn;
            cn=0
        ];
    ];
    If[i==km1, x[[k]]=cn,
        Do[x[[j]]=0,{j,i+1, k}]
    ];
    x
],
CompilationTarget->"C", 
CompilationOptions->{"InlineExternalDefinitions" -> True}, 
RuntimeAttributes->{Listable}
];

This can be used as follows:

In[36]:= ablRanMultinomial[20, {0.5, 0.3, 0.2}]

Out[36]= {12, 7, 1}
$\endgroup$
  • 1
    $\begingroup$ Performance improvement looks impressive - it would be helpful though to know if there was any quality impediment in adopting the above. $\endgroup$ – wolfies Apr 24 at 6:05

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.