octave_packages/statistics-1.1.3/mvtcdf.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} =} mvtcdf (@var{x}, @var{sigma}, @var{nu})
  18 ## @deftypefnx {Function File} {} mvtcdf (@var{a}, @var{x}, @var{sigma}, @var{nu})
  19 ## @deftypefnx {Function File} {[@var{p}, @var{err}] =} mvtcdf (@dots{})
  20 ## Compute the cumulative distribution function of the multivariate
  21 ## Student's t 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{sigma} is the correlation matrix.
  32 ##
  33 ## @item
  34 ## @var{nu} is the degrees of freedom.
  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 ## sigma = [1.0 0.5; 0.5 1.0];
  58 ## nu = 4;
  59 ## p = mvtcdf (x, sigma, nu)
  60 ## @end group
  61 ##
  62 ## @group
  63 ## a = [-inf 0];
  64 ## p = mvtcdf (a, x, sigma, nu)
  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 Student's t distribution
  81
  82 function [p, err] = mvtcdf (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) == 3)
  90     x = varargin{1};
  91     sigma = varargin{2};
  92     nu = varargin{3};
  93     a = -Inf .* ones (size (x));
  94   elseif (length (varargin) == 4)
  95     a = varargin{1};
  96     x = varargin{2};
  97     sigma = varargin{3};
  98     nu = varargin{4};
  99   else
 100     print_usage ();
 101   endif
 102
 103   # Dimension
 104   q = size (sigma, 1);
 105   cases = size (x, 1);
 106
 107   # Check parameters
 108   if (size (x, 2) != q)
 109     error ("mvtcdf: x must have the same number of columns as sigma");
 110   endif
 111
 112   if (any (size (x) != size (a)))
 113     error ("mvtcdf: a must have the same size as x");
 114   endif
 115
 116   if (! isscalar (nu) && (! isvector (nu) || length (nu) != cases))
 117     error ("mvtcdf: nu must be a scalar or a vector with the same number of rows as x");
 118   endif
 119
 120   # Convert to correlation matrix if necessary
 121   if (any (diag (sigma) != 1))
 122     svar = repmat (diag (sigma), 1, q);
 123     sigma = sigma ./ sqrt (svar .* svar');
 124   endif
 125   if (q < 1 || size (sigma, 2) != q || any (any (sigma != sigma')) || min (eig (sigma)) <= 0)
 126     error ("mvtcdf: sigma must be nonempty symmetric positive definite");
 127   endif
 128
 129   nu = nu(:);
 130   c = chol (sigma)';
 131
 132   # Number of integral transformations
 133   n = 1;
 134
 135   p = zeros (cases, 1);
 136   varsum = zeros (cases, 1);
 137
 138   err = ones (cases, 1) .* err_eps;
 139   # Apply crude Monte-Carlo estimation
 140   while any (err >= err_eps)
 141     # Sample from q-1 dimensional unit hypercube
 142     w = rand (cases, q - 1);
 143
 144     # Transformation of the multivariate t-integral
 145     dvev = tcdf ([a(:, 1) / c(1, 1), x(:, 1) / c(1, 1)], nu);
 146     dv = dvev(:, 1);
 147     ev = dvev(:, 2);
 148     fv = ev - dv;
 149     y = zeros (cases, q - 1);
 150     for i = 1:(q - 1)
 151       y(:, i) = tinv (dv + w(:, i) .* (ev - dv), nu + i - 1) .* sqrt ((nu + sum (y(:, 1:(i-1)) .^ 2, 2)) ./ (nu + i - 1));
 152       tf = (sqrt ((nu + i) ./ (nu + sum (y(:, 1:i) .^ 2, 2)))) ./ c(i + 1, i + 1);
 153       dvev = tcdf ([(a(:, i + 1) - c(i + 1, 1:i) .* y(:, 1:i)) .* tf, (x(:, i + 1) - c(i + 1, 1:i) .* y(:, 1:i)) .* tf], nu + i);
 154       dv = dvev(:, 1);
 155       ev = dvev(:, 2);
 156       fv = (ev - dv) .* fv;
 157     endfor
 158
 159     n++;
 160     # Estimate standard error
 161     varsum += (n - 1) .* ((fv - p) .^ 2) ./ n;
 162     err = gamma .* sqrt (varsum ./ (n .* (n - 1)));
 163     p += (fv - p) ./ n;
 164   endwhile
 165
 166 endfunction