====================================================================
John P. Cunningham, Columbia University
Copyright (C) 2011 John P. Cunningham, Philipp Hennig, and Simon Lacoste-Julien

Expectation propagation for multivariate gaussian probabilities
(code for the basic routines from Cunningham, Hennig, and Lacoste-Julien, ArXiv, 2011. http://arxiv.org/abs/1111.6832)

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
====================================================================

Please cite http://arxiv.org/abs/1111.6832 in any relevant publications.  You may also find this work useful, from my collaborators at Columbia: arxiv.org/pdf/1208.4118 . 

(0) All code requires truncNormMoments.  This is a very numerically stable way to use matlab to calculate the moments of truncated univariate normals.  I recommend you use it.  Otherwise, you will probably run into underflow issues in tail probabilities, which are important test cases for Gaussian integrals.  It is one of my 'utils' functions that gets shared reasonably often.

(1) Based on your application, you may only want axis-aligned hyperrectangles... if so, use axisepmgp.m.  It's a straightforward EP implementation, and quicker than the arbitrary polyhedron version.

(2) Calculate some arbitrary polyhedral integrals...  For that, use epmgp.m, which is "EP for Multivariate Gaussian Probabilities".  It's pretty well commented, so you should be able to use it no sweat.  (you can ignore maxSteps and alphaCorrection, which are some ways we are tweaking  the algorithm... see the latter sections of our ArXiv paper).

(3) Also I include genzmgp.m, which is the numerical guy that we compare against.  It is the state of the art and the best proxy to ground truth.  See our ArXiv paper for more details of that, if you are not clear on that accuracy point already.  genzmgp.m is a wrapper for qsclatmvn.m, the latter of which is horribly written.  But for any accuracy test, I would use genz also.

(4) I include a quick sample of how I might run these things.  The 'run_example.m' file will get you started.

(5) Note: I have extended this EP code to work for arbitrary rank one likelihoods (given a function of that form); contact me if you would like to use that code.

Important: if you are using only axis aligned hyperrectangles, epmgp and axisepmgp give precisely the same result.   If they don't, it's a likely user error.  A common mistake is to specify halfspaces instead of box functions.  Use box functions.