In order to calculate the CDF of a multivariate normal, I followed this example (for the univariate case) but cannot interpret the output produced by scipy:
from scipy.stats import norm
import numpy as np
mean = np.array([1,5])
covariance = np.matrix([[1, 0.3 ],[0.3, 1]])
distribution = norm(loc=mean,scale = covariance)
[[ 8.41344746e-01 4.29060333e-04]
[ 9.99570940e-01 1.58655254e-01]]
P (X1 ≤ x1, . . . ,Xn ≤ xn)
After searching a lot, I think this blog entry by Noah H. Silbert describes the only readymade code from a standard library that can be used for computing the cdf for a multivariate normal in Python. Scipy has a way to do it but as mentioned in the blog, it is difficult to find. The approach is based on a paper by Alan Genz’s.
From the blog, this is how it works.
from scipy.stats import mvn import numpy as np low = np.array([-10, -10]) upp = np.array([.1, -.2]) mu = np.array([-.3, .17]) S = np.array([[1.2,.35],[.35,2.1]]) p,i = mvn.mvnun(low,upp,mu,S) print p 0.2881578675080012