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.