Genz-Bretz multivariate normal in Python

I’ve been fighting for some time to try and get Genz-Bretz’s method for calculating orthant probabilities in multivariate normal distributions imported into Python. I downloaded the fortran code from Alan Genz’s site and was unsuccessful in using f2py to link it with Python. However, I discovered the usefulness of the Python ctypes module in linking with shared libraries (see here). So, I compiled the fortran code using

gfortran mvtdstpack.f -shared -o libmvt.so

and then, within Python, did

from ctypes import *
libmvn = cdll.LoadLibrary('./libmvt.so')
pmvn = libmvn.mvtdst_  # the underscore is added while compiling

This allows us access to the function.  The inputs have to be formatted properly, with the use of c_int, c_double and numpy.ctypeslib.ndpointer, and placed in pmvn.argtypes to prototype the function. The variables can then be initialized and passed through the function or subroutine.

For my purposes this took a bit of a learning curve, but I found a nice example online to make the formatting easier, and the results are really quite fast.  I may actually create a larger Fortran library for this project to speed things up, once I profile my program.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s