1 ## Copyright (C) 2008 Arno Onken <asnelt@asnelt.org>
3 ## This program is free software; you can redistribute it and/or modify it under
4 ## the terms of the GNU General Public License as published by the Free Software
5 ## Foundation; either version 3 of the License, or (at your option) any later
8 ## This program is distributed in the hope that it will be useful, but WITHOUT
9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 ## You should have received a copy of the GNU General Public License along with
14 ## this program; if not, see <http://www.gnu.org/licenses/>.
17 ## @deftypefn {Function File} {@var{p} =} mvncdf (@var{x}, @var{mu}, @var{sigma})
18 ## @deftypefnx {Function File} {} mvncdf (@var{a}, @var{x}, @var{mu}, @var{sigma})
19 ## @deftypefnx {Function File} {[@var{p}, @var{err}] =} mvncdf (@dots{})
20 ## Compute the cumulative distribution function of the multivariate
21 ## normal distribution.
23 ## @subheading Arguments
27 ## @var{x} is the upper limit for integration where each row corresponds
31 ## @var{mu} is the mean.
34 ## @var{sigma} is the correlation matrix.
37 ## @var{a} is the lower limit for integration where each row corresponds
38 ## to an observation. @var{a} must have the same size as @var{x}.
41 ## @subheading Return values
45 ## @var{p} is the cumulative distribution at each row of @var{x} and
49 ## @var{err} is the estimated error.
52 ## @subheading Examples
58 ## sigma = [1.0 0.5; 0.5 1.0];
59 ## p = mvncdf (x, mu, sigma)
64 ## p = mvncdf (a, x, mu, sigma)
68 ## @subheading References
72 ## Alan Genz and Frank Bretz. Numerical Computation of Multivariate
73 ## t-Probabilities with Application to Power Calculation of Multiple
74 ## Constrasts. @cite{Journal of Statistical Computation and Simulation},
75 ## 63, pages 361-378, 1999.
79 ## Author: Arno Onken <asnelt@asnelt.org>
80 ## Description: CDF of the multivariate normal distribution
82 function [p, err] = mvncdf (varargin)
84 # Monte-Carlo confidence factor for the standard error: 99 %
89 if (length (varargin) == 1)
92 sigma = eye (size (x, 2));
93 a = -Inf .* ones (size (x));
94 elseif (length (varargin) == 3)
98 a = -Inf .* ones (size (x));
99 elseif (length (varargin) == 4)
112 # Default value for mu
118 if (size (x, 2) != q)
119 error ("mvncdf: x must have the same number of columns as sigma");
122 if (any (size (x) != size (a)))
123 error ("mvncdf: a must have the same size as x");
127 mu = ones (1, q) .* mu;
128 elseif (! isvector (mu) || size (mu, 2) != q)
129 error ("mvncdf: mu must be a scalar or a vector with the same number of columns as x");
132 x = x - repmat (mu, cases, 1);
134 if (q < 1 || size (sigma, 2) != q || any (any (sigma != sigma')) || min (eig (sigma)) <= 0)
135 error ("mvncdf: sigma must be nonempty symmetric positive definite");
140 # Number of integral transformations
143 p = zeros (cases, 1);
144 varsum = zeros (cases, 1);
146 err = ones (cases, 1) .* err_eps;
147 # Apply crude Monte-Carlo estimation
148 while any (err >= err_eps)
149 # Sample from q-1 dimensional unit hypercube
150 w = rand (cases, q - 1);
152 # Transformation of the multivariate normal integral
153 dvev = normcdf ([a(:, 1) / c(1, 1), x(:, 1) / c(1, 1)]);
157 y = zeros (cases, q - 1);
159 y(:, i) = norminv (dv + w(:, i) .* (ev - dv));
160 dvev = normcdf ([(a(:, i + 1) - c(i + 1, 1:i) .* y(:, 1:i)) ./ c(i + 1, i + 1), (x(:, i + 1) - c(i + 1, 1:i) .* y(:, 1:i)) ./ c(i + 1, i + 1)]);
163 fv = (ev - dv) .* fv;
167 # Estimate standard error
168 varsum += (n - 1) .* ((fv - p) .^ 2) ./ n;
169 err = gamma .* sqrt (varsum ./ (n .* (n - 1)));