Cdf of a multivariate normal in Julia

I’m trying to use the cdf of a multivariate normal distribution in Julia using Distributions.jl.

d = Normal(0,1)
d2 = MvNormal([0.;0.], [1. 0.;0. 1.])
x = [1. ; 2.]
e1 = cdf(d,x[1])
e2 = cdf(d,x[2])

works (as long as the mean vector and the var-cov matrix for the multivariate normal are Floats ), I get an error whenever I write
e3 = cdf(d2,x)
Here is the error message:

ERROR: MethodError: no method matching cdf(::Distributions.MvNormal{Float64,PDMats.PDMat{Float64,Array{Float64,2}},Array{Float64,1}}, ::Array{Float64,1})
Closest candidates are:
  cdf(::Distributions.Distribution{Distributions.Univariate,S<:Distributions.ValueSupport}, ::AbstractArray{T,N}) at /Applications/

It seems like the cdf function does not support multivariate distributions. I couldn’t find anything on this in the documentation of Distributions.jl. Is there any other way to find a cdf of a multivariate normal in Julia?

Yes, it seems that evaluating the multivariate normal cdf is not supported yet. I didn’t realize that.

There’s some info on the algorithm here, if you’re willing to code it up:

Or perhaps you could combine PyCall and SciPy…?

For reference there is an open issue

Though it doesn’t look there has been much traction…

Now that I look into PyCall and SciPy, it turns out that scipy.stats have multivariate_normal but doesn’t have cdf defined as a method. Same issue with Julia.

The only suggested solution I’ve found so far is this: but for some reason I can’t find the scipy.stats.mvn that they’re talking about.

Hmm. That’s disappointing. At a pinch you could use Monte Carlo, replacing the probability of your random vector X being less than a fixed vector x with the empirical probability. If speed is not vital it will give you reasonable results.

The latest version of SciPy — version 1.0 — has a cdf method for the multivariate normal. See here.

That’s great, thanks!