Cdf of multivariate normal in distributions.jl

Specific Domains Statistics
22

image

1 Like
23

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

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.

1 Like
24

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

2 Likes
25

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

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.

26

So I compared it with function from 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:

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.