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

["Lacy,Michael" <Michael.Lacy@colostate.edu>]

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/