octave_packages/statistics-1.1.3/mvtrnd.m

   1 ## Copyright (C) 2012  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{x} =} mvtrnd (@var{sigma}, @var{nu})
  18 ## @deftypefnx {Function File} {@var{x} =} mvtrnd (@var{sigma}, @var{nu}, @var{n})
  19 ## Generate random samples from the multivariate t-distribution.
  20 ##
  21 ## @subheading Arguments
  22 ##
  23 ## @itemize @bullet
  24 ## @item
  25 ## @var{sigma} is the matrix of correlation coefficients. If there are any
  26 ## non-unit diagonal elements then @var{sigma} will be normalized.
  27 ##
  28 ## @item
  29 ## @var{nu} is the degrees of freedom for the multivariate t-distribution.
  30 ## @var{nu} must be a vector with the same number of elements as samples to be
  31 ## generated or be scalar.
  32 ##
  33 ## @item
  34 ## @var{n} is the number of rows of the matrix to be generated. @var{n} must be
  35 ## a non-negative integer and corresponds to the number of samples to be
  36 ## generated.
  37 ## @end itemize
  38 ##
  39 ## @subheading Return values
  40 ##
  41 ## @itemize @bullet
  42 ## @item
  43 ## @var{x} is a matrix of random samples from the multivariate t-distribution
  44 ## with @var{n} row samples.
  45 ## @end itemize
  46 ##
  47 ## @subheading Examples
  48 ##
  49 ## @example
  50 ## @group
  51 ## sigma = [1, 0.5; 0.5, 1];
  52 ## nu = 3;
  53 ## n = 10;
  54 ## x = mvtrnd (sigma, nu, n);
  55 ## @end group
  56 ##
  57 ## @group
  58 ## sigma = [1, 0.5; 0.5, 1];
  59 ## nu = [2; 3];
  60 ## n = 2;
  61 ## x = mvtrnd (sigma, nu, 2);
  62 ## @end group
  63 ## @end example
  64 ##
  65 ## @subheading References
  66 ##
  67 ## @enumerate
  68 ## @item
  69 ## Wendy L. Martinez and Angel R. Martinez. @cite{Computational Statistics
  70 ## Handbook with MATLAB}. Appendix E, pages 547-557, Chapman & Hall/CRC, 2001.
  71 ##
  72 ## @item
  73 ## Samuel Kotz and Saralees Nadarajah. @cite{Multivariate t Distributions and
  74 ## Their Applications}. Cambridge University Press, Cambridge, 2004.
  75 ## @end enumerate
  76 ## @end deftypefn
  77
  78 ## Author: Arno Onken <asnelt@asnelt.org>
  79 ## Description: Random samples from the multivariate t-distribution
  80
  81 function x = mvtrnd (sigma, nu, n)
  82
  83   # Check arguments
  84   if (nargin < 2)
  85     print_usage ();
  86   endif
  87
  88   if (! ismatrix (sigma) || any (any (sigma != sigma')) || min (eig (sigma)) <= 0)
  89     error ("mvtrnd: sigma must be a positive definite matrix");
  90   endif
  91
  92   if (!isvector (nu) || any (nu <= 0))
  93     error ("mvtrnd: nu must be a positive scalar or vector");
  94   endif
  95   nu = nu(:);
  96
  97   if (nargin > 2)
  98     if (! isscalar (n) || n < 0 | round (n) != n)
  99       error ("mvtrnd: n must be a non-negative integer")
 100     endif
 101     if (isscalar (nu))
 102       nu = nu * ones (n, 1);
 103     else
 104       if (length (nu) != n)
 105         error ("mvtrnd: n must match the length of nu")
 106       endif
 107     endif
 108   else
 109     n = length (nu);
 110   endif
 111
 112   # Normalize sigma
 113   if (any (diag (sigma) != 1))
 114     sigma = sigma ./ sqrt (diag (sigma) * diag (sigma)');
 115   endif
 116
 117   # Dimension
 118   d = size (sigma, 1);
 119   # Draw samples
 120   y = mvnrnd (zeros (1, d), sigma, n);
 121   u = repmat (chi2rnd (nu), 1, d);
 122   x = y .* sqrt (repmat (nu, 1, d) ./ u);
 123 endfunction
 124
 125 %!test
 126 %! sigma = [1, 0.5; 0.5, 1];
 127 %! nu = 3;
 128 %! n = 10;
 129 %! x = mvtrnd (sigma, nu, n);
 130 %! assert (size (x), [10, 2]);
 131
 132 %!test
 133 %! sigma = [1, 0.5; 0.5, 1];
 134 %! nu = [2; 3];
 135 %! n = 2;
 136 %! x = mvtrnd (sigma, nu, 2);
 137 %! assert (size (x), [2, 2]);