<body><script type="text/javascript"> function setAttributeOnload(object, attribute, val) { if(window.addEventListener) { window.addEventListener('load', function(){ object[attribute] = val; }, false); } else { window.attachEvent('onload', function(){ object[attribute] = val; }); } } </script> <div id="navbar-iframe-container"></div> <script type="text/javascript" src="https://apis.google.com/js/platform.js"></script> <script type="text/javascript"> gapi.load("gapi.iframes:gapi.iframes.style.bubble", function() { if (gapi.iframes && gapi.iframes.getContext) { gapi.iframes.getContext().openChild({ url: 'https://www.blogger.com/navbar.g?targetBlogID\x3d6813476980165976394\x26blogName\x3dThe+Bayesian+Conspiracy\x26publishMode\x3dPUBLISH_MODE_BLOGSPOT\x26navbarType\x3dBLUE\x26layoutType\x3dCLASSIC\x26searchRoot\x3dhttps://bayesianconspiracy.blogspot.com/search\x26blogLocale\x3den\x26v\x3d2\x26homepageUrl\x3dhttp://bayesianconspiracy.blogspot.com/\x26vt\x3d-948109813487012623', where: document.getElementById("navbar-iframe-container"), id: "navbar-iframe", messageHandlersFilter: gapi.iframes.CROSS_ORIGIN_IFRAMES_FILTER, messageHandlers: { 'blogger-ping': function() {} } }); } }); </script>

The Bayesian Conspiracy

The Bayesian Conspiracy is a multinational, interdisciplinary, and shadowy group of scientists. It is rumored that at the upper levels of the Bayesian Conspiracy exist nine silent figures known only as the Bayes Council.

This blog is produced by James Durbin, a Bayesian Initiate.

Recent

Gibbs sampler in Groovy

Stamped: 5:32 PM |

I recently read a couple of nice articles by Darren Wilkinson about implementing MCMC in various languages. The posts are here and here. Wilkinson apparently uses, or is considering using, Python for a lot of prototyping and C for a lot of his actual MCMC runs. However, since he feels that Java is in some ways nicer than C, and almost as fast, he has been using Java some also for the final MCMC runs. I thought I'd see how Groovy performed on this task.
A Groovy version of Wilkinson's Gibbs sampler is given here: 

#!/usr/bin/env groovy 

import static java.lang.Math.*
import cern.jet.random.tdouble.engine.*
import cern.jet.random.tdouble.*

N=20000
thin=500
rngEngine= new DoubleMersenneTwister(new Date())
rngN= new Normal(0.0,1.0,rngEngine)
rngG= new Gamma(1.0,1.0,rngEngine)
x=y=0.0
println("Iter x y")
for(i in 1..N){
  for(j in 1..thin){
    x=rngG.nextDouble(3.0,y*y+4)
    y=rngN.nextDouble(1.0/(x+1),1.0/sqrt(x+1))
  }
  println("$i $x $y")
}

You'd run it like:
./Gibbs.gv > data.tab

On a Macbook Pro 2.2 GHz Intel Core 2 Duo, comparing Wilkinson's Java (1.6.1) and Python (2.7.2) versions with my Groovy (1.8.0) version, I get the following times:
  • java: 5.23 sec
  • python: 1m49s
  • groovy: 16.5s
Not bad for a dynamic scripting language! Now, there are lots of things one can do to speed up the Python code, from using pypy to calling native functions and so on, so this is nowhere close to optimal for Python but, if you're going to ultimately translate your prototype into Java or C anyway, why bother with Python at all when you can rapid-prototype in Groovy then do the simple translation to Java, using the same libraries, later? Moreover, at one third the speed of Java, many times you won't even need to bother with the final speed-up step.

Labels: , ,

« Home
Using Groovlets in jQuery tutorial. »
gcsvsql »
Executing arbitrary SQL on CSV files. »
Groovy++ »
Create/open OmniOutliner notebooks from command line. »
Weka: getting predictions from cross validation »
Groovy in Google App Engine »
Invalid duplicate class definition »
Groovy compared to Perl »
Groovy for bioinformatics. »

You can leave your response or bookmark this post to del.icio.us by using the links below.

Post a Comment | Bookmark | Go to end |
  • Anonymous Anonymous:  

    Did you try compiling the Groovy to see what difference there was between the compile Java and compiled Groovy execution times?

    There's also supposed to be some Java 7-related work which will speed up dynamic languages. I think there's some work that the dynamic languages will have to do to take advantage of it though. (July 25, 2011 at 5:43 PM) top

  • Blogger James Durbin:  

    I have never seen compiling Groovy help much. I gave it a whirl here and, as usual, it didn't make any noticeable difference. I think for processes that run more than a few seconds, JIT compiling pretty much smooths out any differences there would be. I do compile lots of Groovy code, to call it from Java, but I gave up on that as a way to improve performance.

    Groovy++ really does have a big performance gain, but you lose a lot of Groovy goodness in the process. It might be useful in this case, though... maybe I'll give it a try and update the post. (July 25, 2011 at 9:23 PM) top

  • Blogger James Durbin:  

    Actually, I tried Groovy++ on this example and it only had a small impact, maybe 10%. I've seen examples where the improvement was large. (July 27, 2011 at 2:00 PM) top

  • Blogger James Durbin:  

    Hold on... I'm going to try something:

    The probability of getting \(k\) heads when flipping \(n\) coins is

    \[P(E) = {n \choose k} p^k (1-p)^{ n-k}\] (July 29, 2011 at 11:23 PM) top