function [y,err] = mvncdf(varargin)
%MVNCDF Multivariate normal cumulative distribution function (cdf).
% Y = MVNCDF(X) returns the cumulative probability of the multivariate
% normal distribution with zero mean and identity covariance matrix,
% evaluated at each row of X. Rows of the N-by-D matrix X correspond to
% observations or points, and columns correspond to variables or
% coordinates. Y is an N-by-1 vector.
%
% Y = MVNCDF(X,MU,SIGMA) returns the cumulative probability of the
% multivariate normal distribution with mean MU and covariance SIGMA,
% evaluated at each row of X. MU is a 1-by-D vector, and SIGMA is a
% D-by-D symmetric, positive definite matrix. MU can also be a scalar
% value, which MVNCDF replicates to match the size of X. If the
% covariance matrix is diagonal, containing variances along the diagonal
% and zero covariances off the diagonal, SIGMA may also be specified as a
% 1-by-D matrix containing just the diagonal. Pass in the empty
% matrix for MU to use its default value when you want to only specify
% SIGMA.
%
% The multivariate normal cumulative probability at X is defined as the
% probability that a random vector V, distributed as multivariate normal,
% will fall within the semi-infinite rectangle with upper limits defined by
% X, i.e., Pr{V(1)<=X(1), V(2)<=X(2), ... V(D)<=X(D)}.
%
% Y = MVNCDF(XL,XU,MU,SIGMA) returns the multivariate normal cumulative
% probability evaluated over the rectangle (hyper-rectangle for D>2)
% with lower and upper limits defined by XL and XU, respectively.
%
% [Y,ERR] = MVNCDF(...) returns an estimate of the error in Y. For
% bivariate and trivariate distributions, MVNCDF uses adaptive quadrature on
% a transformation of the t density, based on methods developed by Drezner
% and Wesolowsky, and by Genz, as described in the references. The default
% absolute error tolerance for these cases is 1e-8. For four or more
% dimensions, MVNCDF uses a quasi-Monte Carlo integration algorithm based on
% methods developed by Genz and Bretz, as described in the references. The
% default absolute error tolerance for these cases is 1e-4. For
% univariate distributions and when SIGMA is specified as a diagonal,
% numerical integration is not used and ERR returns NaN.
%
% [...] = MVNCDF(...,OPTIONS) specifies control parameters for the numerical
% integration, used to compute Y when D > 1 and is not specified as a
% diagonal. This argument can be created by a call to STATSET. Choices
% of STATSET parameters are:
%
% 'TolFun' - Maximum absolute error tolerance. Default is 1e-8