Tuesday, July 14, 2009

Multivariate normal distribution in Python

I could not find a Python function to evaluate the multivariate normal distribution in Python. Here's one that gives equivalent results to the dmvnorm function in the mvtnorm package for R. It's something that works. I've not had time or need yet to fix it up.

b: A vector
mean: The mean of the elements in b (same dimensions as b)
cov: The covariance matrix of the multivariate normal distribution

License: GPL version 2
http://www.gnu.org/licenses/gpl-2.0.html
k = b.shape[0]
part1 = numpy.exp(-0.5*k*numpy.log(2*numpy.pi))
part2 = numpy.power(numpy.linalg.det(cov),-0.5)
dev = b-mean
part3 = numpy.exp(-0.5*numpy.dot(numpy.dot(dev.transpose(),numpy.linalg.inv(cov)),dev))
dmvnorm = part1*part2*part3

9 comments:

Anonymous said...

Amiable post and this fill someone in on helped me alot in my college assignement. Say thank you you as your information.

Anonymous said...

Interesting blog you got here. I'd like to read more concerning this topic. Thanx for sharing this information.

Anonymous said...

Genial post and this fill someone in on helped me alot in my college assignement. Say thank you you as your information.

Anonymous said...

i honestly adore all your posting choice, very helpful.
don't quit as well as keep writing since it just that is worth to read it.
excited to look into a lot more of your own posts, goodbye ;)

Anonymous said...

Well your article helped me very much in my college assignment. Hats afar to you post, will look ahead for more interdependent articles soon as its sole of my pet issue to read.

mrhoward said...

This has actually been very helpful, as I also haven't been able to find a proper implementation of the PDF of multivariate normal distributions. I was just wondering if there was any particular reason for you using

exp(-0.5*log(2*pi))

instead of

(2*pi)**-0.5.

(which gives the same result)

mrhoward said...

Of course, I meant

(2*pi)**-k/2.

lmf said...

@mrhoward: I'm not sure to be honest. It's been a while since I wrote this post. My recollection is that I found an implementation in R and translated it to Python.

buy erection pills said...

Thank you for the good writeup. It in fact was a amusement account it. Look advanced to more added agreeable from you! By the way, how could we communicate?