MVN CDF -- have it coded, need help getting integrating into Distributions.jl

Summary
I have written a Julia function to compute the Multivariate Normal integral.
I need a partner who is an expert in Distributions.jl to properly integrate it into the code base.

Introduction

I have written a function in Julia to compute the central Multivariate Normal probabilities:

\Phi_k(\mathbf{a},\mathbf{b};\mathbf{\Sigma} ) = \frac{1}{\sqrt{\left | \mathbf{\Sigma} \right |{(2\pi )^k}}}\int_{a_1}^{b_1}\int_{a_2}^{b_2}\begin{align*} &...\end{align*} \int_{a_k}^{b_k}e^{^{-\frac{1}{2}}\mathbf{x}^t{\mathbf{\Sigma }}^{-1}\boldsymbol{\textbf{x}}}dx_k...dx_1

where \Sigma is a symmetric, positive definite k x k covariance matrix. Vector a is the lower integration limit, and b is the upper integration limit. As this is the “central” MVN, \mu = \vec{0}.

This is a translation of the MATLAB qsimvnv function written by Prof. Alan Genz.

Examples:

julia> r = [4 3 2 1;3 5 -1 1;2 -1 4 2;1 1 2 5]
4×4 Array{Int64,2}:
 4   3   2  1
 3   5  -1  1
 2  -1   4  2
 1   1   2  5
julia> a = [-Inf; -Inf; -Inf; -Inf]
4-element Array{Float64,1}:
 -Inf
 -Inf
 -Inf
 -Inf
julia> b = [1; 2; 3; 4 ]
4-element Array{Int64,1}:
 1
 2
 3
 4
julia> (p,e)=qsimvnv(r,a,b)
(0.6049240592769761, 0.0016812280451818128)

In this example, 0.604924… is the computed probability, and 0.00168122… is the estimated error of the calculation. The function uses quasi-random integration to compute the integral, so you get a slightly different result each time it’s run. The results of the function have been checked for correctness by comparing them to R, MATLAB, and published results.

julia> (p,e)=qsimvnv(r,a,b)
(0.6060880398660634, 0.001540997560720256)
julia> (p,e)=qsimvnv(r,a,b)
(0.6060880398660634, 0.001540997560720256)
julia> (p,e)=qsimvnv(r,a,b)
(0.6067495056529963, 0.0013883605551699985)

You can increase the iterations to reduce the error. The default is 1000*dimension of \Sigma. Decrease in error is roughly linearly with increase in iterations:

julia> (p,e)=qsimvnv(r,a,b;m=5000)
(0.6061177664363478, 0.0018974040888668662)

julia> (p,e)=qsimvnv(r,a,b;m=6000)
(0.605451604851963, 0.0009217826063167869)

julia> (p,e)=qsimvnv(r,a,b;m=7000)
(0.6057328187356109, 0.0009514371151440057)

julia> (p,e)=qsimvnv(r,a,b;m=80000)
(0.6056921305482306, 0.0002024445811897692)

julia> (p,e)=qsimvnv(r,a,b;m=9000)
(0.6060092106313155, 0.0009367451258961832)

julia> (p,e)=qsimvnv(r,a,b;m=10000)
(0.605821023964934, 0.0005674428432139305)

This could be the basis of a MVN cdf function for Distributions.jl which operates with the MvNormal objects.

julia> r=float(r)
4×4 Array{Float64,2}:
 4.0   3.0   2.0  1.0
 3.0   5.0  -1.0  1.0
 2.0  -1.0   4.0  2.0
 1.0   1.0   2.0  5.0
julia> mu=[0, 0, 0, 0]
4-element Array{Int64,1}:
 0
 0
 0
 0
julia> mvn1 = MvNormal(mu, r)
FullNormal(
dim: 4
μ: [0.0, 0.0, 0.0, 0.0]
Σ: [4.0 3.0 2.0 1.0; 3.0 5.0 -1.0 1.0; 2.0 -1.0 4.0 2.0; 1.0 1.0 2.0 5.0]
)
julia> mvn_cdf(D::AbstractMvNormal, a::AbstractVector, b::AbstractVector; m::Union{Integer, Nothing} = nothing) = qsimvnv(D.Σ.mat,a .- D.μ, b .- D.μ;m=m)
mvn_cdf (generic function with 1 method)

julia> mvn_cdf(mvn1,a,b;m=5000)
(0.6056845764526123, 0.0015238130294856398)

What I have

  • Julia implementation of qsimvnv()
  • Julia implementation of Prof. Genz’s function that 1) scales and re-arranges the integration vectors, and 2) computes the cholesky decomposition of \Sigma. This is called by qsimvnv().
  • Several test cases, consisting of matrices of various sizes, the expected probability and error estimate results, and tolerances for each. This could form the basis of test code.

What I Don’t Have

  • A good understanding of the structure of the Distributions.jl code and all of the different MvNormal objects. When I look at the code, I can’t even figure out where I would put the two functions I have.
  • Programming skills with Julia types, methods, multiple dispatch. The “mvn_cdf” example above is the pinnacle of my achievement, and it doesn’t work unless you specify m= in it (i.e., not specifying m to use the default results in an error).
  • Git experience

What I Need

  • I need a partner who is an expert in the structure of Distributions.jl and the MvNormal objects, who knows how to integrate it into the code base and write the methods for all of the different MvNormal objects. I would expect his partner would also identify improvements to the code I’ve written for qsimvnv() and _chlrdr (the cholesky decomposition function).
  • Note that MvNormal objects can have nonzero mean vectors; to use qsimvnv(), the \mu vector needs to be subtracted from the integration limits a and b.
  • Some of the MvNormal objects use Correlation matrices instead of Covariance matrices. These need to be converted to covariance matrices before calling the function.
  • Given a template to work from, I can probably contribute to writing the test code. I can also run tests, and I can compare the results to MATLAB as needed.

Next Steps
If someone wants to partner up on this, let me know what to do next. I could create a PR to dump the code in to the /src directory, but my partner would have to fix it afterwards. Or … let me know how to proceed.

Appendix - Statistics Notes

One way to define the Multivariate Normal CDF is for “simultaneous X” – that is, a for a given random vector \vec{x}, the probability that every component of \vec{x} is less than a reference vector, \vec{x_0}:

F(\vec{x_{0}}) = P(\vec{x} \leqslant \vec{x_{0}})

which can be generalized to:

F(\vec{a},\vec{b}) = P(\vec{a} \leqslant \vec{x} \leqslant \vec{b})

Fn(\vec{a},\vec{b};\mu,\Sigma) = P(ai < Xi < bi;i=1...n)

F_{n}(\vec{a},\vec{b};\vec{\mu} ,\Sigma ) = P(a_{i}<X_{i}<b_{i}; i=1...n) = \int_{a_{1}}^{b_{1}}\int_{a_{2}}^{b_{2}}\cdots \int_{a_{i}}^{b_{i}}f(x_{1},x_{2},...,x_{i})dx_{1}dx_{2}\cdots dx_{i}

We define:

F_{n}(\vec{a},\vec{b};\vec{\mu} ,\Sigma ) = \Phi_k(\mathbf{a},\mathbf{b},\mathbf{\Sigma} )

where:

\Phi_k(\mathbf{a},\mathbf{b};\mathbf{\Sigma} ) = \frac{1}{\sqrt{\left | \mathbf{\Sigma} \right |{(2\pi )^k}}}\int_{a_1}^{b_1}\int_{a_2}^{b_2}\begin{align*} &...\end{align*} \int_{a_k}^{b_k}e^{^{-\frac{1}{2}}\mathbf{x}^t{\mathbf{\Sigma }}^{-1}\boldsymbol{\textbf{x}}}dx_k...dx_1

Note that you’re only really computing the CDF if \vec{a} = -\infty.
If \vec{a}=\vec{0} and \vec{b} = +\vec{\infty}, this is called the orthant probability.

Non-Central MVN

The non-central MVN integral is given by:

\Phi_k(\mathbf{a},\mathbf{b};\mathbf{\Sigma},\delta ) = \frac{1}{\sqrt{\left | \mathbf{\Sigma} \right |{(2\pi )^k}}}\int_{a_1}^{b_1}\int_{a_2}^{b_2}\begin{align*} &...\end{align*} \int_{a_k}^{b_k}e^{^{-\frac{1}{2}}\mathbf{(x-\delta)}^t{\mathbf{\Sigma }}^{-1}\mathbf{(x-\delta)}}dx_k...dx_1

where \delta denotes the k x 1 non-centrality vector with -{\infty} < \delta_i < +{\infty}.

The non-central MVN integral can determined from the central MVN with shifted integration limits:

\Phi_k(\mathbf{a},\mathbf{b};\mathbf{\Sigma},\delta ) = \Phi_k(\mathbf{a-\delta},\mathbf{b-\delta };\mathbf{\Sigma})

Standardized MVN

The non-central MVN probability \Phi_k(\mathbf{a},\mathbf{b};\mathbf{\Sigma},\delta ) can also be transformed into the standardized MVN probability. The standardized MVN has \vec{\mu} = \vec{0} and \Sigma = I_n.

Let \mathbf{D} denote the diagonal matrix of the square roots of the diagonal entries of \Sigma. The correlation matrix \mathbf{R} is defined by \mathbf{\Sigma = DRD}. Then the transformation \mathbf{x = Dy + \delta} reduces the non-central MVN probability to

\Phi_k(\mathbf{a},\mathbf{b};\mathbf{\Sigma},\delta ) = \Phi_k(\mathbf{D^{-1}(a-\delta),D^{-1}(b-\delta)};\mathbf{R})

Mahalanobis Distance

The other way to think about the MVN probability is that of a sample vector \vec{x_0} falls within the ellipsoid defined by MVN(\vec{\mu},\Sigma). This has a closed form solution. If you’re interested in this, see Distribution of observation-level Mahalanobis distance.

5 Likes

First pass, Distributions.jl doesn’t implement the CDF for any MultivariateDistribution yet, so there’s no agreed upon standard there. StatsFuns.jl might be the appropriate home for this in the medium-term. However, you have to be careful about the licensing issues that might surround "a translation of the MATLAB qsimvnv".

2 Likes

Before you invest time in this, it’s worth making sure that you have the relevant IP rights from Alan Genz to include a translation of their code in a MIT-licensed project like Distributions.jl.

If that’s all set, I think you don’t need a full time partner to get started with this if you’re passionate and interested to experiment on your own. I would do the following:

  1. Locate your local copy of Distributions. You should be able to find this using something like
using Distributions
@which Bernoulli(0.5)
  1. Once you’ve found the directory, look for the mvnormal.jl file.

  2. Add a method like the following:

function cdf(d::MvNormal, x::AbstractArray)
    # YOUR CODE GOES HERE
end

Once you fill in code for cdf, you should be able to run it locally (since you’re editing your copy of Distributions) and test it.

If you’re confused by Git, you could wait to get your cdf function working and then submit a pull request using the steps from https://www.youtube.com/watch?v=0gshzIYqkcY&feature=youtu.be.

By the time you’ve done that, you’ll have learned a bunch and could finish the work with much less help from someone with commit rights to Distributions.jl. As @johnczito notes, the real blocker is getting people to agree what CDF for multivariate distributions does.

5 Likes

Prof. Genz allows redistribution of his code:

% Copyright © 2013, Alan Genz, All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided the following conditions are met:
% 1. Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
% 2. Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in
% the documentation and/or other materials provided with the
% distribution.
% 3. The contributor name(s) may not be used to endorse or promote
% products derived from this software without specific prior
% written permission.

1 Like

Great. And I remember now that a Genz port of the bivariate normal cdf already exists in StatsFuns.jl/src/tvpack.jl. It happens not to be exported. And in that case I think someone just emailed him for permission.

If the contribution is a perfectly numerical function that computes multivariate normal probabilities over hyperrectangles, then I’d say the best place to put that is StatsFuns.jl. And that process could be fairly uncomplicated. It sounds like you’re already pretty close. Observe the proper licensing issues, write some thorough tests, and then follow @johnmyleswhite’s advice above.

Like I said, trying to go direct to Distributions.jl could be a challenge because there’s no agreed-upon convention for how cdf(d::MultivariateDistribution, x::AbstractVector) ought to behave. So your PR could get caught up in lengthy discussions about that. But no matter what convention is ultimately settled on, this function will be relevant, so if you put it in StatsFuns.jl, it will be “ready and waiting.”

1 Like

If the contribution is a perfectly numerical function that computes multivariate normal probabilities over hyperrectangles, then I’d say the best place to put that is StatsFuns.jl .

This is how things are structured for a lot of the univariates anyway. For things like gamma, beta, etc. The function that actually computes the cdf is in StatsFuns.jl. Distributions.jl is just a front end.

It’s not so simple, I think.

For MvNormal, we have:
const IsoNormal = MvNormal{ScalMat, Vector{Float64}}
const DiagNormal = MvNormal{PDiagMat, Vector{Float64}}
const FullNormal = MvNormal{PDMat, Vector{Float64}}

const ZeroMeanIsoNormal{Axes} = MvNormal{ScalMat, Zeros{Float64,1,Axes}}
const ZeroMeanDiagNormal{Axes} = MvNormal{PDiagMat, Zeros{Float64,1,Axes}}
const ZeroMeanFullNormal{Axes} = MvNormal{PDMat, Zeros{Float64,1,Axes}}

The MvNormal object stores its data differently depending on how it’s defined. \Sigma may be stored as an array, a PDMat, a PDDiagMat, a ScalMat, a Vector, or a real-valued number. If the covariance matrix is stored in a compact form, a full k x k covariance matrix needs to be created before calling qsimvnv().

In addition, there are the MvNoralCannon types, which store the potential vector \mathbf{h} and precision matrix \mathbf{J} instead of \mu and \Sigma. The MvNormalCannon types has just as many variations as the MvNormal type.

If you look at function fit_mle() in MvNormal.jl, you’ll see all of the different cases that need to be covered. There’s a lot in there I don’t understand – objects and type system in general – but also objects like MvNormalKnownCov.

I thinking adding it to StatsFun.jl is a lot more manageable – more of the work of forming a \Sigma with the proper integration limits is put on the user.

Yeah, sorting out those front-end issues is secondary. It will all resolve to a function like what you have in the end.

It’s true that doing a good job isn’t that simple, but to unblock yourself, you could start by simplifying the problem to avoid dealing with types you don’t understand:

function cdf(d::MvNormal, x::AbstractArray)
    mu = Vector(d.μ)
    sigma = Matrix(d.Σ)
    # YOUR CODE HERE USES mu AND sigma

Eventually you’d likely need to handle all of the complexity, but you could make progress by simplifying the problem to a narrower space that you’re comfortable operating in and then gradually increasing complexity to provide the performance/accuracy you’d lose through the simplification.

I don’t contribute to Distributions.jl anymore and can’t speak to how much support you could rely on from the current maintainers, but I suspect blocking on support from a volunteer partner is likely to slow you down substantially more than trying to fumble your way through things and making mistakes that will eventually teach you how to do things through simple trial and error.

All that said, I think @johnczito’s approach is a good one – it’s another simplification that will help you get unblocked.

1 Like

Just stumbled on this issue. I’d be really good if your implementation was made widely available; I was looking for this at one point and I was disappointed not to find it. I can help out a little bit by looking at the code and suggesting how to make it more Julian and more ready for submission to StatsFun.jl. (Not that I’ve ever contibuted to StatsFun.jl, but, you know, general direction :slight_smile: ).
You’re welcome to ping me if/when you get the code onto github.

1 Like

Any update on this?

I’m very frustrated that cdf function in Distributions.jl does not have a method for multivariate distributions.

1 Like

Seven months have passed and I am wondering whether this function is ready for public release? If not, may I request for a private copy for my own non-commercial research use? Many thanks.

I think it would be best to make the code for the cdf available in its own mini-package (with tests), which Distributions.jl and other packages could just use.

IP issues can be avoided by proceeding directly from the book

“”"
qsimvnv(Σ,a,b;m=iterations)
Computes the Multivariate Normal probability integral using a quasi-random rule
with m points for positive definite covariance matrix Σ, mean [0,…], with lower
integration limit vector a and upper integration limit vector b.

\\Phi_k(\\mathbf{a},\\mathbf{b},\\mathbf{\\Sigma} ) = \\frac{1}{\\sqrt{\\left | \\mathbf{\\Sigma}  \\right |{(2\\pi )^k}}}\\int_{a_1}^{b_1}\\int_{a_2}^{b_2}\\begin{align*}
 &...\\end{align*} \\int_{a_k}^{b_k}e^{^{-\\frac{1}{2}}\\mathbf{x}^t{\\mathbf{\\Sigma }}^{-1}\\boldsymbol{\\mathbf{x}}}dx_k...dx_1

Probability p is output with error estimate e.

Arguments

  • Σ::AbstractArray: positive-definite covariance matrix of MVN distribution
  • a::Vector: lower integration limit column vector
  • b::Vector: upper integration limit column vector
  • m::Int64: number of integration points (default 1000*dimension)

Example

julia> r = [4 3 2 1;3 5 -1 1;2 -1 4 2;1 1 2 5]
julia> a=[-Inf; -Inf; -Inf; -Inf]
julia> b = [1; 2; 3; 4 ]
julia> m = 5000
julia> (p,e)=qsimvnv(r,a,b;m=m)
(0.605219554009911, 0.0015718064928452481)

Results will vary slightly from run-to-run due to the quasi-Monte Carlo
algorithm.

Non-central MVN distributions (with non-zero mean) can use this function by adjusting
the integration limits. Subtract the mean vector, μ, from each
integration vector.

Example

julia> #non-central MVN
julia> Σ=[4 2;2 3]
julia> μ = [1;2]
julia> a=[-Inf; -Inf]
julia> b=[2; 2]
julia> (p,e) = qsimvnv(Σ,a-μ,b-μ)
(0.4306346895870772, 0.00015776288569406053)

“”"
function qsimvnv(Σ,a,b;m=nothing)
#= rev 1.13

This function uses an algorithm given in the paper
"Numerical Computation of Multivariate Normal Probabilities", in
 J. of Computational and Graphical Stat., 1(1992), pp. 141-149, by
Alan Genz, WSU Math, PO Box 643113, Pullman, WA 99164-3113
Email : alangenz@wsu.edu
The primary references for the numerical integration are
"On a Number-Theoretical Integration Method"
H. Niederreiter, Aequationes Mathematicae, 8(1972), pp. 304-11, and
"Randomization of Number Theoretic Methods for Multiple Integration"
R. Cranley and T.N.L. Patterson, SIAM J Numer Anal, 13(1976), pp. 904-14.

Re-coded in Julia from the MATLAB function qsimvnv(m,r,a,b)

Alan Genz is the author the MATLAB qsimvnv() function.
Alan Genz software website: http://archive.is/jdeRh
Source code to MATLAB qsimvnv() function: http://archive.is/h5L37
% QSIMVNV(m,r,a,b) and _chlrdr(r,a,b)
%
% Copyright (C) 2013, Alan Genz,  All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided the following conditions are met:
%   1. Redistributions of source code must retain the above copyright
%      notice, this list of conditions and the following disclaimer.
%   2. Redistributions in binary form must reproduce the above copyright
%      notice, this list of conditions and the following disclaimer in
%      the documentation and/or other materials provided with the
%      distribution.
%   3. The contributor name(s) may not be used to endorse or promote
%      products derived from this software without specific prior
%      written permission.
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
% "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
% FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
% COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
% INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
% BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
% OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
% ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
% TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF USE
% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%

Julia dependencies
Distributions
PDMats
Primes
Random
LinearAlgebra

=#

if isnothing(m)
	m = 1000*size(Σ,1)  # default is 1000 * dimension
end

# check for proper dimensions
n=size(Σ,1)
nc=size(Σ,2) 	# assume square Cov matrix nxn
# check dimension > 1
n >= 2   || throw(ErrorException("dimension of Σ must be 2 or greater. Σ dimension: $(size(Σ))"))
n == nc  || throw(DimensionMismatch("Σ matrix must be square. Σ dimension: $(size(Σ))"))

# check dimensions of lower vector, upper vector, and cov matrix match
(n == size(a,1) == size(b,1)) || throw(DimensionMismatch("iconsistent argument dimensions. Sizes: Σ $(size(Σ))  a $(size(a))  b $(size(b))"))

# check that a and b are column vectors; if row vectors, fix it
if size(a,1) < size(a,2)
	a = transpose(a)
end
if size(b,1) < size(b,2)
	b = transpose(b)
end

# check that lower integration limit a < upper integration limit b for all elements
all(a .<= b) || throw(ArgumentError("lower integration limit a must be <= upper integration limit b"))

# check that Σ is positive definate; if not, print warning
isposdef(Σ) || @warn "covariance matrix Σ fails positive definite check"

# check if Σ, a, or b contains NaNs
if any(isnan.(Σ)) || any(isnan.(a)) || any(isnan.(b))
	p = NaN
	e = NaN
	return (p,e)
end

# check if a==b
if a == b
	p = 0.0
	e = 0.0
	return (p,e)
end

# check if a = -Inf & b = +Inf
if all(a .== -Inf) && all(b .== Inf)
	p = 1.0
	e = 0.0
	return (p,e)
end

# check input Σ, a, b are floats; otherwise, convert them
if eltype(Σ)<:Signed
	Σ = float(Σ)
end

if eltype(a)<:Signed
	a = float(a)
end

if eltype(b)<:Signed
	b = float(b)
end

##################################################################
#
# Special cases: positive Orthant probabilities for 2- and
# 3-dimesional Σ have exact solutions. Integration range [0,∞]
#
##################################################################

if all(a .== zero(eltype(a))) && all(b .== Inf) && n <= 3
	Σstd = sqrt.(diag(Σ))
	Rcorr = cov2cor(Σ,Σstd)

	if n == 2
		p = 1/4 + asin(Rcorr[1,2])/(2π)
		e = eps()
	elseif n == 3
		p = 1/8 + (asin(Rcorr[1,2]) + asin(Rcorr[2,3]) + asin(Rcorr[1,3]))/(4π)
		e = eps()
	end

	return (p,e)

end

##################################################################
#
# get lower cholesky matrix and (potentially) re-ordered integration vectors
#
##################################################################

(ch,as,bs)=_chlrdr(Σ,a,b) # ch =lower cholesky; as=lower vec; bs=upper vec

##################################################################
#
# quasi-Monte Carlo integration of MVN integral
#
##################################################################

### setup initial values
ai=as[1]
bi=bs[1]
ct=ch[1,1]

unitnorm = Normal() # unit normal distribution
rng=RandomDevice()

# if ai is -infinity, explicity set c=0
# implicitly, the algorith classifies anythign > 9 std. deviations as infinity
if ai > -9*ct
	if ai < 9*ct
		c1 = cdf.(unitnorm,ai/ct)
	else
		c1 = 1.0
	end
else
	c1 = 0.0
end

# if bi is +infinity, explicity set d=0
if bi > -9*ct
	if bi < 9*ct
		d1 = cdf(unitnorm,bi/ct)
	else
		d1 = 1.0
	end
else
	d1 = 0.0
end

#n=size(Σ,1) 	# assume square Cov matrix nxn
cxi=c1			# initial cxi; genz uses ci but it conflicts with Lin. Alg. ci variable
dci=d1-cxi		# initial dcxi
p=0.0			# probablity = 0
e=0.0			# error = 0

# Richtmyer generators
ps=sqrt.(primes(Int(floor(5*n*log(n+1)/4)))) # Richtmyer generators
q=ps[1:n-1,1]
ns=12
nv=Int(max(floor(m/ns),1))

Jnv    = ones(1,nv)
cfill  = transpose(fill(cxi,nv)) 	# evaulate at nv quasirandom points row vec
dpfill = transpose(fill(dci,nv))
y      = zeros(n-1,nv)				# n-1 rows, nv columns, preset to zero

#=Randomization loop for ns samples
 j is the number of samples to integrate over,
     but each with a vector nv in length
 i is the number of dimensions, or integrals to comptue =#

for j in 1:ns					# loop for ns samples
	c  = copy(cfill)
	dc = copy(dpfill)
	pv = copy(dpfill)
		for i in 2:n
			x=transpose(abs.(2.0 .* mod.((1:nv) .* q[i-1] .+ rand(rng),1) .- 1))	 # periodizing transformation
			# note: the rand() is not broadcast -- it's a single random uniform value added to all elements
			y[i-1,:] = quantile.(unitnorm,c .+ x.*dc)
			s = transpose(ch[i,1:i-1]) * y[1:i-1,:]
			ct=ch[i,i]										# ch is cholesky matrix
			ai=as[i] .- s
			bi=bs[i] .- s
			c=copy(Jnv)										# preset to 1 (>9 sd, +∞)
			d=copy(Jnv)										# preset to 1 (>9 sd, +∞)

			c[findall(x -> isless(x,-9*ct),ai)] .= 0.0		# If < -9 sd (-∞), set to zero
			d[findall(x -> isless(x,-9*ct),bi)] .= 0.0		# if < -9 sd (-∞), set to zero
			tstl = findall(x -> isless(abs(x),9*ct),ai)		# find cols inbetween -9 and +9 sd (-∞ to +∞)
			c[tstl] .= cdf.(unitnorm,ai[tstl]/ct)			# for those, compute Normal CDF
			tstl = (findall(x -> isless(abs(x),9*ct),bi))	# find cols inbetween -9 and +9 sd (-∞ to +∞)
			d[tstl] .= cdf.(unitnorm,bi[tstl]/ct)
			@. dc=d-c
			@. pv=pv * dc
		end # for i=
	d=(mean(pv)-p)/j
	p += d
	e=(j-2)*e/j+d^2
end # for j=

e=3*sqrt(e) 	# error estimate is 3 times standard error with ns samples

return (p,e)  	# return probability value and error estimate

end # function qsimvnv

“”"
Computes permuted lower Cholesky factor c for R which may be singular,
also permuting integration limit vectors a and b.

Arguments

r		matrix			Matrix for which to compute lower Cholesky matrix
						when called by mvn_cdf(), this is a covariance matrix

a		vector			column vector for the lower integration limit
						algorithm may permutate this vector to improve integration
						accuracy for mvn_cdf()

b		vector			column vector for the upper integration limit
						algorithm may pertmutate this vector to improve integration
						accuracy for mvn_cdf()

Output
tuple An a tuple with 3 returned arrays:
1 - lower Cholesky root of r
2 - lower integration limit (perhaps permutated)
3 - upper integration limit (perhaps permutated)

Examples

r = [1 0.25 0.2; 0.25 1 0.333333333; 0.2 0.333333333 1]
a = [-1; -4; -2]
b = [1; 4; 2]

(c, ap, bp) = _chlrdr(r,a,b)

result:
Lower cholesky root:
c = [ 1.00 0.0000 0.0000,
0.20 0.9798 0.0000,
0.25 0.2892 0.9241 ]
Permutated upper input vector:
ap = [-1, -2, -4]
Permutated lower input vector:
bp = [1, 2, 4]

Related Functions

mvn_cdf - multivariate Normal CDF function makes use of this function

“”"
function _chlrdr(Σ,a,b)

# Rev 1.13

# define constants
# 64 bit machien error 1.0842021724855e-19 ???
# 32 bit machine error 2.220446049250313e-16 ???
ep = 1e-10 # singularity tolerance
if Sys.WORD_SIZE == 64
    fpsize=Float64
    ϵ = eps(0.0) # 64-bit machine error
else
    fpsize=Float32
    ϵ = eps(0.0f0) # 32-bit machine error
end

if !@isdefined sqrt2π
    sqrt2π = √(2π)
end

# unit normal distribution
unitnorm = Normal()

n = size(Σ,1) # covariance matrix n x n square

ckk = 0.0
dem = 0.0
am = 0.0
bm = 0.0
ik = 0.0

if eltype(Σ)<:Signed
    c = copy(float(Σ))
else
    c = copy(Σ)
end

if eltype(a)<:Signed
    ap = copy(float(a))
else
    ap = copy(a)
end

if eltype(b)<:Signed
    bp = copy(float(b))
else
    bp = copy(b)
end

d=sqrt.(diag(c))
for i in 1:n
    if d[i] > 0.0
        c[:,i] /= d[i]
        c[i,:] /= d[i]
        ap[i]=ap[i]/d[i]     # ap n x 1 vector
        bp[i]=bp[i]/d[i]     # bp n x 1 vector
    end
end

y=zeros(fpsize,n) # n x 1 zero vector to start

for k in 1:n
    ik = k
    ckk = 0.0
    dem = 1.0
    s = 0.0
    #pprinta(c)
    for i in k:n
        if c[i,i] > ϵ  # machine error
            cii = sqrt(max(c[i,i],0))

            if i>1 && k>1
                s=(c[i,1:(k-1)].*y[1:(k-1)])[1]
            end

            ai=(ap[i]-s)/cii
            bi=(bp[i]-s)/cii
            de = cdf(unitnorm,bi) - cdf(unitnorm,ai)

            if de <= dem
                ckk = cii
                dem = de
                am = ai
                bm = bi
                ik = i
            end
        end # if c[i,i]> ϵ
    end # for i=
    i = n

    if ik>k
        ap[ik] , ap[k] = ap[k] , ap[ik]
        bp[ik] , bp[k] = bp[k] , bp[ik]

        c[ik,ik] = c[k,k]

        if k > 1
            c[ik,1:(k-1)] , c[k,1:(k-1)] = c[k,1:(k-1)] , c[ik,1:(k-1)]
        end

        if ik<n
            c[(ik+1):n,ik] , c[(ik+1):n,k] = c[(ik+1):n,k] , c[(ik+1):n,ik]
        end

        if k<=(n-1) && ik<=n
            c[(k+1):(ik-1),k] , c[ik,(k+1):(ik-1)] = transpose(c[ik,(k+1):(ik-1)]) , transpose(c[(k+1):(ik-1),k])
        end
    end # if ik>k

    if ckk > k*ep
        c[k,k]=ckk
        if k < n
            c[k:k,(k+1):n] .= 0.0
        end

        for i in (k+1):n
            c[i,k] /= ckk
            c[i:i,(k+1):i] -= c[i,k]*transpose(c[(k+1):i,k])
        end

        if abs(dem)>ep
            y[k] = (exp(-am^2/2)-exp(-bm^2/2))/(sqrt2π*dem)
        else
            if am<-10
                y[k] = bm
            elseif bm>10
                y[k]=am
            else
                y[k]=(am+bm)/2
            end
        end # if abs
    else
        c[k:n,k] .== 0.0
        y[k] = 0.0
    end # if ckk>ep*k
end # for k=

return (c, ap, bp)

end # function _chlrdr

function testmvn(;m=nothing)

#Typical Usage/Example Code
#Example multivariate Normal CDF for various vectors
println()

# from MATLAB documentation
 r = [4 3 2 1;3 5 -1 1;2 -1 4 2;1 1 2 5]
 	a=[-Inf; -Inf; -Inf; -Inf] # -inf for each
	b = [1; 2; 3; 4 ]
	m = 5000
	#m=4000 # rule of thumb: 1000*(number of variables)
	(myp,mye)=qsimvnv(r,a,b;m=m)
	println("Answer should about 0.6044 to 0.6062")
	println(myp)
	println("The Error should be <= 0.001 - 0.0014");
	println(mye)

	r=[1 0.6 0.333333;
	   0.6 1 0.733333;
	   0.333333 0.733333 1]
	r  = [  1   3/5   1/3;
		  3/5    1    11/15;
		  1/3  11/15    1]
	a=[-Inf;-Inf;-Inf]
	b=[1;4;2]
	#m=3000;
	(myp,mye)=qsimvnv(r,a,b;m=4000)
	println()
	println("Answer shoudl be about 0.82798")
	# answer from Genz paper 0.827984897456834
	println(myp)
	println("The Error should be <= 2.5-05")
	println(mye)

	r=[1 0.25 0.2; 0.25 1 0.333333333; 0.2 0.333333333 1]
	a=[-1;-4;-2]
	b=[1;4;2]
	#m=3000;
	(myp,mye)=qsimvnv(r,a,b;m=4000);
	println()
	println("Answer should be about 0.6537")
	println(myp)
	println("The Error should be <= 2.5-05")
	println(mye)

	# Genz problem 1.5 p. 4-5  & p. 63
	# correct answer F(a,b) = 0.82798
    r = [1/3 3/5 1/3;
         3/5 1.0 11/15;
         1/3 11/15 1.0]
    a = [-Inf; -Inf; -Inf]
    b = [1; 4; 2]
	(myp,mye)=qsimvnv(r,a,b;m=4000)
	println()
	#println("Answer shoudl be about 0.82798")
	println("Answer should be 0.9432")
	println(myp)
	println("The error should be < 6e-05")
	println(mye)

	# Genz page 63 uses a different r Matrix
	r = [1 0 0;
		3/5 1 0;
		1/3 11/15 1]
		a = [-Inf; -Inf; -Inf]
		b = [1; 4; 2]
	(myp,mye)=qsimvnv(r,a,b;m=4000)
	println()
	println("Answer shoudl be about 0.82798")
	println(myp)
	println("The error should be < 6e-05")
	println(mye)
	# mystery solved - he used the wrong sigma Matrix
	# when computing the results on page 6


	# singular cov Example
	r = [1 1 1; 1 1 1; 1 1 1]
	a = [-Inf, -Inf, -Inf]
	b = [1, 1, 1]
	(myp,mye)=qsimvnv(r,a,b)
	println()
	println("Answer should be 0.841344746068543")
	println(myp)
	println("The error should be 0.0")
	println(mye)
	println("Cov matrix is singular")
	println("Problem reduces to a univariate problem with")
	println("p = cdf.(Normal(),1) = ",cdf.(Normal(),1))

	# 5 dimensional Example
	#c = LinearAlgebra.tri(5)
	r = [1 1 1 1 1;
		 1 2 2 2 2;
		 1 2 3 3 3;
		 1 2 3 4 4;
		 1 2 3 4 5]
	a = [-1,-2,-3,-4,-5]
	b = [2,3,4,5,6]
	(myp,mye)=qsimvnv(r,a,b)
	println()
	println("Answer should be ~ 0.7613")
	# genz gives 0.4741284  p. 5 of book
	# Julia, MATLAB, and R all give ~0.7613 !
	println(myp)
	println("The error should be < 0.001")
	println(mye)

	# genz reversed the integration limits when computing
	# see p. 63
	a = sort!(a)
	b = 1 .- a
	(myp,mye)=qsimvnv(r,a,b)
	println()
	println("Answer should be ~ 0.4741284")
	# genz gives 0.4741284  p. 5 of book
	# now matches with reversed integration limits
	println(myp)
	println("The error should be < 0.001")
	println(mye)

	# positive orthant of above
	a = [0,0,0,0,0]
	(myp,mye)=qsimvnv(r,a,b)
	println()
	println("Answer should be ~  0.11353418")
	# genz gives 0.11353418   p. 6 of book
	println(myp)
	println("The error should be < 0.001")
	println(mye)

	# now with a = -inf
	a = [-Inf,-Inf,-Inf,-Inf,-Inf]
	(myp,mye)=qsimvnv(r,a,b)
	println()
	println("Answer should be ~ 0.81031455")
	# genz gives 0.81031455  p. 6 of book
	println(myp)
	println("The error should be < 0.001")
	println(mye)

	# eight dimensional Example
	r = [1 1 1 1 1 1 1 1;
		 1 2 2 2 2 2 2 2;
		 1 2 3 3 3 3 3 3;
		 1 2 3 4 4 4 4 4;
		 1 2 3 4 5 5 5 5;
		 1 2 3 4 5 6 6 6;
		 1 2 3 4 5 6 7 7;
		 1 2 3 4 5 6 7 8]
	a = -1*[1,2,3,4,5,6,7,8]
	b = [2,3,4,5,6,7,8,9]
	(myp,mye)=qsimvnv(r,a,b)
	println()
	println("Answer should be ~ 0.7594")
	# genz gives 0.32395   p. 6 of book
	# MATLAB gives 0.7594362
	println(myp)
	println("The error should be < 0.001")
	println(mye)


	# eight dim orthant
	a=[0,0,0,0,0,0,0,0]
	b=[Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf]
	(myp,mye)=qsimvnv(r,a,b)
	println()
	println("Answer should be ~ 0.196396")
	# genz gives 0.076586  p. 6 of book
	# MATLB gives 0.196383
	println(myp)
	println("The error should be < 0.001")
	println(mye)

	# eight dim with a = -inf
	a=-Inf*[1,1,1,1,1,1,1,1]
	b = [2,3,4,5,6,7,8,9]
	#b = [0,0,0,0,0,0,0,0]
	(myp,mye)=qsimvnv(r,a,b)
	println()
	println("Answer should be ~ 0.9587")
	# genz gives 0.69675    p. 6 of book
	# MATLAB gives 0.9587
	println(myp)
	println("The error should be < 0.001")
	println(mye)

end # testmvn

3 Likes

@blackeneth Very nice! Thank you for this, I’ve been looking for this for a long time.

Do you have this anywhere available on GitHub, so we could use it in the meantime before it’s integrated to Distributions or StatsFuns?

1 Like

Thank you @blackeneth for making it available to us!