Bookmark and Share

Notice: On March 31, it was announced that Statalist is moving from an email list to a forum. The old list will shut down on April 23, and its replacement, is already up and running.

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

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

From   "Lacy,Michael" <>
To   "" <>
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 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_on(2)   //built-in
	Prob2 = binomialp(N,X,PI)
 mreldif(Prob1, Prob2)  // check that I'm actually calculating the same thing in both instances
local N = 20  // number of trials
local nsize = 30  // size of matrix
local reps = 1000   
mata: timer_clear()
mata: test(`N', `nsize', `reps')



Mike Lacy
Dept. of Sociology
Colorado State University
Fort Collins CO 80523-1784

*   For searches and help try:

© Copyright 1996–2015 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   Site index