Hi, yet another useless benchmarking exercise: I was playing around a little comparing the speed of various matrix languages with respect to doing stochastic simulations. The codes below draw repeated random samples, compute their means, and find the 95% percent quantile of the empirical distribution of the means. I distinguish between a "naive" loop programming variant and one working only with matrix functions. Here are the results on Ubuntu Linux 9.10 (I ran the scripts/programs several times and they are pretty robust): Octave 3.0.5: loop ca. 3.3 seconds vectorized ca. 0.8 seconds Python 2.6.4/Numpy 1.3.0/matplotlib 0.99(?): loop ca. 2.1 seconds vectorized ca. 1.5 seconds gretl 1.8.6cvs: loop ca. 2.6 seconds vectorized ca. 2.0 seconds Ox 5.1: loop ca. 1.1 seconds vectorized ca. 0.7 seconds Ox confirms its reputation about speed, but for me it was surprising that Octave is almost as fast with the vectorized code. (But Octave has the slowest loops.) Gretl seems to have the smallest relative gain going from loops to vectorized code. But considering that gretl is quite new in the field of matrix programming I think it is doing fine. cheers, sven <Octave-code> length = 1000; iterations = 10000; whichp = 0.95; tic; # looping variant averages = []; for i=1:iterations; mrandom = randn(length,1); averages = [averages ; mean(mrandom)]; endfor; myquant = quantile(averages,whichp); toc; printf ("Result: %f\n", myquant); tic; # vectorized variant mrandom2 = randn(length,iterations); averages2 = mean(mrandom2); # row vector myquant2 = quantile(averages2',whichp); toc; printf ("Result: %f\n", myquant2); </Octave-code> <Python-code> from numpy import matlib as nm from matplotlib import mlab as ml import time length = 1000 iterations = 10000 whichp = 0.95 starttime = time.time() # looping variant averages = nm.empty((0,1)) for i in range(iterations): mrandom = nm.randn(length,1) averages = nm.vstack((averages, nm.mean(mrandom))) myquant = ml.prctile(averages,whichp*100) # per cents endtime1 = time.time() print "Result: " + str(myquant) print "Execution time loop variant, length " + str(length) + ", iterations " + str(iterations) + ": " + str(endtime1-starttime) # vectorized variant starttime2 = time.time() mrandom2 = nm.randn(length,iterations) averages2 = nm.mean(mrandom2, axis=0) myquant2 = ml.prctile(averages2,whichp*100) endtime2 = time.time() print "Result: " + str(myquant2) print "Execution time vectorized variant, length " + str(length) + ", iterations " + str(iterations) + ": " + str(endtime2-starttime2) </Python-code> <gretl-code> length = 1000 iterations = 10000 whichp = 0.95 set echo off set messages off set stopwatch # looping variant (non-vectorized) matrix averages = {} # will be col vector loop iterations matrix mrandom = mnormal(length,1) matrix averages = averages | meanc(mrandom) end loop matrix myquant = quantile(averages, whichp) time = $stopwatch printf "Result: %f\n", myquant printf "Execution time loop variant, length %d, iterations %d: %f\n", length, iterations, time time = $stopwatch # don't count printing overhead # vectorized variant matrix mrandom2 = mnormal(length,iterations) matrix averages2 = transp(meanc(mrandom2)) matrix myquant2 = quantile(averages2, whichp) time = $stopwatch printf "Result: %f\n", myquant2 printf "Execution time vectorized variant, length %d, iterations %d: %f\n", length, iterations, time </gretl-code> </Ox-code> #include <oxstd.h> const decl length = 1000; const decl iterations = 10000; const decl whichp = 0.95; main(){ decl averages, averages2; decl myquant, myquant2; decl mrandom, mrandom2; decl starttime, endtime; decl i; starttime = timer(); // looping variant averages = <>; for (i=0; i<iterations; ++i) { mrandom = rann(length,1); averages = averages | meanc(mrandom); } myquant = quantilec(averages,whichp); endtime = timer(); print("Result: ", myquant, "\n"); print("elapsed time: ", timespan(starttime,endtime)); starttime = timer(); // vectorized variant mrandom2 = rann(length,iterations); averages2 = meanc(mrandom2); // row vector myquant2 = quantiler(averages2,whichp); endtime = timer(); print("\nResult: ", myquant2,"\n"); print("elapsed time: ", timespan(starttime,endtime),"\n"); } </Ox-code>