Improving speed of Binomial and Multinomial Random Draws

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?