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} =} 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.
23 ## @subheading Arguments
27 ## @var{x} is the upper limit for integration where each row corresponds
31 ## @var{sigma} is the correlation matrix.
34 ## @var{nu} is the degrees of freedom.
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
57 ## sigma = [1.0 0.5; 0.5 1.0];
59 ## p = mvtcdf (x, sigma, nu)
64 ## p = mvtcdf (a, x, sigma, nu)
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 Student's t distribution
82 function [p, err] = mvtcdf (varargin)
84 # Monte-Carlo confidence factor for the standard error: 99 %
89 if (length (varargin) == 3)
93 a = -Inf .* ones (size (x));
94 elseif (length (varargin) == 4)
108 if (size (x, 2) != q)
109 error ("mvtcdf: x must have the same number of columns as sigma");
112 if (any (size (x) != size (a)))
113 error ("mvtcdf: a must have the same size as x");
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");
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');
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");
132 # Number of integral transformations
135 p = zeros (cases, 1);
136 varsum = zeros (cases, 1);
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);
144 # Transformation of the multivariate t-integral
145 dvev = tcdf ([a(:, 1) / c(1, 1), x(:, 1) / c(1, 1)], nu);
149 y = zeros (cases, 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);
156 fv = (ev - dv) .* fv;
160 # Estimate standard error
161 varsum += (n - 1) .* ((fv - p) .^ 2) ./ n;
162 err = gamma .* sqrt (varsum ./ (n .* (n - 1)));