# Cdf of multivariate normal in distributions.jl

**URL:** <https://discourse.julialang.org/t/cdf-of-multivariate-normal-in-distributions-jl/4773>\
**Category:** Statistics\
**Created:** [July 10, 2017, 9:21pm UTC](https://discourse.julialang.org/t/cdf-of-multivariate-normal-in-distributions-jl/4773 "2017-07-10T21:21:42Z")\
**Posts on this page:** 5\
**Page:** 2

<div class="post-metadata">

**Author:** ![Fuzeq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fuzeq/32/49148_2.png) [@Fuzeq](https://discourse.julialang.org/u/Fuzeq)\
**Post date:** [April 20, 2023, 2:22pm UTC](https://discourse.julialang.org/t/cdf-of-multivariate-normal-in-distributions-jl/4773/22 "2023-04-20T14:22:34Z")

</div>

> [@Dan](#):
>
> Can you add a plot of the difference also?

![image](https://global.discourse-cdn.com/julialang/original/3X/2/c/2ce74f254792f7aae5494bf47fc0842f0b0a0ba7.png)

---

<div class="post-metadata">

**Author:** ![Fuzeq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fuzeq/32/49148_2.png) [@Fuzeq](https://discourse.julialang.org/u/Fuzeq)\
**Post date:** [April 20, 2023, 5:07pm UTC](https://discourse.julialang.org/t/cdf-of-multivariate-normal-in-distributions-jl/4773/23 "2023-04-20T17:07:42Z")

</div>

> [@Dan](#):
>
> Can you add a plot of the difference also?

Also here is function that calculates mean value of absolute difference from n random selected points:

```julia
function diff_avg(n,μ, Σ)
    x_low = μ[1]-3*Σ[1,1]
    y_low = μ[2]-3*Σ[2,2]
    x_high = μ[1]+3*Σ[1,1]
    y_high = μ[2]+3*Σ[2,2]
    dist = MultivariateNormal(μ, Σ)
    sample = rand(dist, n) 
    xs = rand(Uniform(x_low,x_high),n)
    ys = rand(Uniform(y_low,y_high),n)
    ecdf(x,y) = sum((sample[1,i]<x && sample[2,i]<y) for i in 1:n)/n
    _mv_norm_cdf(x,y) = hcubature(x -> pdf(dist,x), [x_low, y_low], [x, y])[1]
    return mean(abs.(_mv_norm_cdf.(xs,ys)-ecdf.(xs,ys)))
end
diff_avg(10^5,μ,Σ)

```

It’s less than 0.005 for 10^5 points.

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [April 20, 2023, 8:20pm UTC](https://discourse.julialang.org/t/cdf-of-multivariate-normal-in-distributions-jl/4773/24 "2023-04-20T20:20:20Z")

</div>

Hi! Just try to calculate more than 3 dimensional problem… that’s why QMC used for qsimvnv in Matlab.

For 2 and 3 dimensional case you can explore:

Drezner, Z. “Computation of the Trivariate Normal Integral.” _Mathematics of Computation_. Vol. 63, 1994, pp. 289–294.  
Drezner, Z., and G. O. Wesolowsky. “On the Computation of the Bivariate Normal Integral.” _Journal of Statistical Computation and Simulation_. Vol. 35, 1989, pp. 101–107.

I didn’t check, but maybe that methods is slightly better than QMC. If you make compaction - PR welcome to [MvNormalCDF.jl](https://github.com/PharmCat/MvNormalCDF.jl)

---

<div class="post-metadata">

**Author:** ![Fuzeq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fuzeq/32/49148_2.png) [@Fuzeq](https://discourse.julialang.org/u/Fuzeq)\
**Post date:** [April 21, 2023, 9:30am UTC](https://discourse.julialang.org/t/cdf-of-multivariate-normal-in-distributions-jl/4773/25 "2023-04-21T09:30:27Z")

</div>

> [@PharmCat](#):
>
> Hi! Just try to calculate more than 3 dimensional problem…

I think this should work, but I’m not 100% sure:

```julia
using HCubature, Distributions,Statistics, StatsBase, LinearAlgebra

params(d::MvNormal) = (d.μ, d.Σ)
function mv_norm_cdf(dist::MvNormal, coords::Vector)
    dims = length(dist)
    μ, Σ = params(dist)
    lows = [μ[i] - 5*Σ[i,i] for i in 1:dims]
    return hcubature(x -> pdf(dist,x), lows,coords)[1]
end

n = 4
μ = [0.0 for i in 1:n]
Σ = diagm(ones(n))
dist = MvNormal(μ,Σ)

mv_norm_cdf(dist,[0,0,0,0])

```

Of course there is some error one directly from HCubature and also one, because we cut some mass to shorten calculation time.

---

<div class="post-metadata">

**Author:** ![Fuzeq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fuzeq/32/49148_2.png) [@Fuzeq](https://discourse.julialang.org/u/Fuzeq)\
**Post date:** [April 24, 2023, 8:27am UTC](https://discourse.julialang.org/t/cdf-of-multivariate-normal-in-distributions-jl/4773/26 "2023-04-24T08:27:22Z")

</div>

So I compared it with function from [MvNormalCDF.jl](https://github.com/PharmCat/MvNormalCDF.jl). It’s slightly slower, but it’s always returns the same value, because I don’t use QMC.

Here is (probably) final version:

```julia
using HCubature, Distributions,Statistics, StatsBase, LinearAlgebra

params(d::MvNormal) = (d.μ, d.Σ)

function mvnorm_cdf(dist::MvNormal, coords::Vector)
    dims = length(dist)
    μ, Σ = params(dist)
    lows = [μ[i] - 6*Σ[i,i] for i in 1:dims]
    return hcubature(x -> pdf(dist,x), lows,coords, maxevals=10^5*dims)[1]
end

```

There is still posibility to manipulate with lower integration limits, lower they are longer calculation time, but more accurate result.

[Previous page](https://discourse.julialang.org/t/cdf-of-multivariate-normal-in-distributions-jl/4773.md?page=1)