octave_packages/statistics-1.1.3/mvncdf.m

   1 ## Copyright (C) 2008 Arno Onken <asnelt@asnelt.org>
   2 ##
   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
   6 ## version.
   7 ##
   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
  11 ## details.
  12 ##
  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/>.
  15
  16 ## -*- texinfo -*-
  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.
  22 ##
  23 ## @subheading Arguments
  24 ##
  25 ## @itemize @bullet
  26 ## @item
  27 ## @var{x} is the upper limit for integration where each row corresponds
  28 ## to an observation.
  29 ##
  30 ## @item
  31 ## @var{mu} is the mean.
  32 ##
  33 ## @item
  34 ## @var{sigma} is the correlation matrix.
  35 ##
  36 ## @item
  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}.
  39 ## @end itemize
  40 ##
  41 ## @subheading Return values
  42 ##
  43 ## @itemize @bullet
  44 ## @item
  45 ## @var{p} is the cumulative distribution at each row of @var{x} and
  46 ## @var{a}.
  47 ##
  48 ## @item
  49 ## @var{err} is the estimated error.
  50 ## @end itemize
  51 ##
  52 ## @subheading Examples
  53 ##
  54 ## @example
  55 ## @group
  56 ## x = [1 2];
  57 ## mu = [0.5 1.5];
  58 ## sigma = [1.0 0.5; 0.5 1.0];
  59 ## p = mvncdf (x, mu, sigma)
  60 ## @end group
  61 ##
  62 ## @group
  63 ## a = [-inf 0];
  64 ## p = mvncdf (a, x, mu, sigma)
  65 ## @end group
  66 ## @end example
  67 ##
  68 ## @subheading References
  69 ##
  70 ## @enumerate
  71 ## @item
  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.
  76 ## @end enumerate
  77 ## @end deftypefn
  78
  79 ## Author: Arno Onken <asnelt@asnelt.org>
  80 ## Description: CDF of the multivariate normal distribution
  81
  82 function [p, err] = mvncdf (varargin)
  83
  84   # Monte-Carlo confidence factor for the standard error: 99 %
  85   gamma = 2.5;
  86   # Tolerance
  87   err_eps = 1e-3;
  88
  89   if (length (varargin) == 1)
  90     x = varargin{1};
  91     mu = [];
  92     sigma = eye (size (x, 2));
  93     a = -Inf .* ones (size (x));
  94   elseif (length (varargin) == 3)
  95     x = varargin{1};
  96     mu = varargin{2};
  97     sigma = varargin{3};
  98     a = -Inf .* ones (size (x));
  99   elseif (length (varargin) == 4)
 100     a = varargin{1};
 101     x = varargin{2};
 102     mu = varargin{3};
 103     sigma = varargin{4};
 104   else
 105     print_usage ();
 106   endif
 107
 108   # Dimension
 109   q = size (sigma, 1);
 110   cases = size (x, 1);
 111
 112   # Default value for mu
 113   if (isempty (mu))
 114     mu = zeros (1, q);
 115   endif
 116
 117   # Check parameters
 118   if (size (x, 2) != q)
 119     error ("mvncdf: x must have the same number of columns as sigma");
 120   endif
 121
 122   if (any (size (x) != size (a)))
 123     error ("mvncdf: a must have the same size as x");
 124   endif
 125
 126   if (isscalar (mu))
 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");
 130   endif
 131
 132   x = x - repmat (mu, cases, 1);
 133
 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");
 136   endif
 137
 138   c = chol (sigma)';
 139
 140   # Number of integral transformations
 141   n = 1;
 142
 143   p = zeros (cases, 1);
 144   varsum = zeros (cases, 1);
 145
 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);
 151
 152     # Transformation of the multivariate normal integral
 153     dvev = normcdf ([a(:, 1) / c(1, 1), x(:, 1) / c(1, 1)]);
 154     dv = dvev(:, 1);
 155     ev = dvev(:, 2);
 156     fv = ev - dv;
 157     y = zeros (cases, q - 1);
 158     for i = 1:(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)]);
 161       dv = dvev(:, 1);
 162       ev = dvev(:, 2);
 163       fv = (ev - dv) .* fv;
 164     endfor
 165
 166     n++;
 167     # Estimate standard error
 168     varsum += (n - 1) .* ((fv - p) .^ 2) ./ n;
 169     err = gamma .* sqrt (varsum ./ (n .* (n - 1)));
 170     p += (fv - p) ./ n;
 171   endwhile
 172
 173 endfunction