Drawing from a multinomial distribution
Claire Runser
Hi all,
In the model I am currently working on, I need to be able to draw from a multinomial distribution. Unfortunately, Eidos does not possess a function such as rmultinom that would allow me to do that. I have tried to code one in Eidos in multiple ways, but they all run very slowly. See underneath the two ways I have tried to do it that run the fastest.
This one uses a series of binomial draws :
function (integer) multinom (integer$ size, float probs){
s=size;
rho=1;
k=length(probs);
counts=c(1:k);
for (i in 0:(k-2)){
if (rho != 0){
p=probs[i]/rho;
if(p>1.0){p=1.0;}
counts[i] = rbinom(1,s,p);}
else
{counts[i] = 0;}
s=s-counts[i];
rho=rho-probs[i];
}
counts[k-1]=s;
return counts;}
This one draws size samples from k categories, with probabilities probs, then counts how many times each category appears :
function(integer)multinom2(integer$ k, integer$ size, float probs) {
draws=sample(x=0:(k-1),size=size,replace=T,weights=probs);
counts=tabulate(draws,maxbin=99);
return counts;}
No matter how i try to code it, it runs very slowly. It really is a problem because I call this function very often. Has anybody encountered this problem before ? Any way to code this in a way that runs fast ? Thanks!
Claire Runser.
Ben Haller
I've just added a rmultinom() function to Eidos for you. The only reason it wasn't there already is that nobody has ever asked for it before. :->
It is indeed much faster. With your multinom2() function:
> executeLambda("multinom2(5, 100000000, c(0.1, 0.5, 0.1, 5.0, 0.25));", T)
// ********** executeLambda() elapsed time: 4.03563
1679701 8403316 1678966 84035998 4202019 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
With the new rmultinom() function:
> executeLambda("rmultinom(1, 100000000, c(0.1, 0.5, 0.1, 5.0, 0.25));", T);
// ********** executeLambda() elapsed time: 8.6e-05
[,0]
[0,] 1678721
[1,] 8403512
[2,] 1679716
[3,] 84036912
[4,] 4201139
Hopefully that solves your speed problem. :-> Note that the new rmultinom() matches rmultinom() in R pretty much exactly, and returns a matrix of results for (if you wish) multiple draws. So the semantics of it are a bit different from your multinom2() function, and your code will need to be adjusted accordingly.
This has just been pushed to the `master` branch of the SLiM repo. Instructions for building from source are in chapter 2 of the SLiM manual. It will be released in the next release of SLiM, but given that I just released SLiM 5.2 a couple of days ago, that will probably be a while from now – probably months. So if you want this new function now, you'll need to build from sources; but that's pretty straightforward to do.
Let me know if you notice any problems with it. Good luck and happy modeling!
Cheers, -B. Benjamin C. Haller Messer Lab Cornell University
--
SLiM forward genetic simulation: http://messerlab.org/slim/
---
You received this message because you are subscribed to the Google Groups "slim-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to slim-discuss...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/slim-discuss/b45712ea-653c-4e57-a706-3d9e38254888n%40googlegroups.com.