The Binomial and Multinomial random number generators in Mathematica are fast if multiple draws are needed from the same distribution, i.e., when the distribution parameters do not change across the draws. This can be seen by generating, for example, 500,000 draws from the Binomial distribution.

In[30]:= AbsoluteTiming[ RandomVariate[BinomialDistribution[100, 0.6], 500000];]
Out[30]= {0.017365, Null}

However, the speed is slower, relative to some other systems, when the parameters change across draws. Such a situation arises often when performing certain Monte Carlo simulations.

For example, if we have a vector nvec that contains the number of binomial 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 can be addressed via an implementation of fast random number generators for the Binomial as 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 use, as in the following function.

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 another function from the Communications of ACM paper referred above.

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 for different values of the parameters 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}
];

The above functions result in considerable improvement in speed of random number generation for the Binomial.

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 by leveraging the listability property of the functions as follows.

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

Out[33]= {0.156881, 32.5337}

The above functions can then be used for generating multinomial draws. The multinomial distribution is specified via two parameters. n is the number of trials, and p is a vector of probabilities that sum to 1. Multinomial variates can be generated using the following function:

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}