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

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

 From "Lacy,Michael" To "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