Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

From |
"Lacy,Michael" <Michael.Lacy@colostate.edu> |

To |
"statalist@hsphsun2.harvard.edu" <statalist@hsphsun2.harvard.edu> |

Subject |
st: Why is Mata's binomialp() applied to a matrix slower than DIY version with loops? |

Date |
Sun, 23 Oct 2011 22:30:04 +0000 |

In the context of needing to calculate a sizeable matrix of binomial probabilities using Mata, I was very surprised that doing this by using loops and my own code to calculate the probability appears to be about twice as fast as doing this using the built-in function -binomialp() applied to a matrix. This is completely counter to my experience that using matrix operations in Mata is typically 4 or 5 times faster than looping. This matters to me because this particular operation is a large part of total execution time in the program in which I need to do this. Can someone make sense of this? And, any other guidance on gaining efficiency in this task? Some illustrative code follows. mata: mata clear // DIY vs. built-in binomial probability; N trials, X successes, probability PI void test(real scalar N, real scalar nsize, real scalar reps) { real scalar i, j, r real matrix X, PI, Prob1, Prob2 // Make data for N, X, and PI X = J(nsize, nsize, .) Prob1 = J(nsize, nsize, .) Prob2 = J(nsize, nsize, .) PI = runiform(nsize, nsize) X = 1 :+ floor(N * runiform(nsize, nsize)) // Calculate binomial probabilities DIY vs. built-in, for many reps. for (r =1; r <= reps; r++) { timer_on(1) //DIY for (i = 1; i <= nsize ; i++ ) { for (j = 1; j <= nsize; j++) { Prob1[i,j] = comb(N,X[i,j]) * PI[i,j]^X[i,j] * (1-PI[i,j])^(N-X[i,j]) } } timer_off(1) // timer_on(2) //built-in Prob2 = binomialp(N,X,PI) timer_off(2) } mreldif(Prob1, Prob2) // check that I'm actually calculating the same thing in both instances } end // local N = 20 // number of trials local nsize = 30 // size of matrix local reps = 1000 // mata: timer_clear() mata: test(`N', `nsize', `reps') mata:timer() ---------------------------------------- Regards, Mike Lacy Dept. of Sociology Colorado State University Fort Collins CO 80523-1784 * * For searches and help try: * http://www.stata.com/help.cgi?search * http://www.stata.com/support/statalist/faq * http://www.ats.ucla.edu/stat/stata/

- Prev by Date:
**st: New command -des2- on SSC** - Next by Date:
**st: My axis(2)-title is ignored when using by() in twoway scatter?** - Previous by thread:
**st: New command -des2- on SSC** - Next by thread:
**st: My axis(2)-title is ignored when using by() in twoway scatter?** - Index(es):