octave_packages/statistics-1.1.3/copularnd.m

   1 ## Copyright (C) 2012  Arno Onken
   2 ##
   3 ## This program is free software: you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation, either version 3 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program.  If not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {@var{x} =} copularnd (@var{family}, @var{theta}, @var{n})
  18 ## @deftypefnx {Function File} {} copularnd (@var{family}, @var{theta}, @var{n}, @var{d})
  19 ## @deftypefnx {Function File} {} copularnd ('t', @var{theta}, @var{nu}, @var{n})
  20 ## Generate random samples from a copula family.
  21 ##
  22 ## @subheading Arguments
  23 ##
  24 ## @itemize @bullet
  25 ## @item
  26 ## @var{family} is the copula family name. Currently, @var{family} can be
  27 ## @code{'Gaussian'} for the Gaussian family, @code{'t'} for the Student's t
  28 ## family, or @code{'Clayton'} for the Clayton family.
  29 ##
  30 ## @item
  31 ## @var{theta} is the parameter of the copula. For the Gaussian and Student's t
  32 ## copula, @var{theta} must be a correlation matrix. For bivariate copulas
  33 ## @var{theta} can also be a correlation coefficient. For the Clayton family,
  34 ## @var{theta} must be a vector with the same number of elements as samples to
  35 ## be generated or be scalar.
  36 ##
  37 ## @item
  38 ## @var{nu} is the degrees of freedom for the Student's t family. @var{nu} must
  39 ## be a vector with the same number of elements as samples to be generated or
  40 ## be scalar.
  41 ##
  42 ## @item
  43 ## @var{n} is the number of rows of the matrix to be generated. @var{n} must be
  44 ## a non-negative integer and corresponds to the number of samples to be
  45 ## generated.
  46 ##
  47 ## @item
  48 ## @var{d} is the number of columns of the matrix to be generated. @var{d} must
  49 ## be a positive integer and corresponds to the dimension of the copula.
  50 ## @end itemize
  51 ##
  52 ## @subheading Return values
  53 ##
  54 ## @itemize @bullet
  55 ## @item
  56 ## @var{x} is a matrix of random samples from the copula with @var{n} samples
  57 ## of distribution dimension @var{d}.
  58 ## @end itemize
  59 ##
  60 ## @subheading Examples
  61 ##
  62 ## @example
  63 ## @group
  64 ## theta = 0.5;
  65 ## x = copularnd ("Gaussian", theta);
  66 ## @end group
  67 ##
  68 ## @group
  69 ## theta = 0.5;
  70 ## nu = 2;
  71 ## x = copularnd ("t", theta, nu);
  72 ## @end group
  73 ##
  74 ## @group
  75 ## theta = 0.5;
  76 ## n = 2;
  77 ## x = copularnd ("Clayton", theta, n);
  78 ## @end group
  79 ## @end example
  80 ##
  81 ## @subheading References
  82 ##
  83 ## @enumerate
  84 ## @item
  85 ## Roger B. Nelsen. @cite{An Introduction to Copulas}. Springer, New York,
  86 ## second edition, 2006.
  87 ## @end enumerate
  88 ## @end deftypefn
  89
  90 ## Author: Arno Onken <asnelt@asnelt.org>
  91 ## Description: Random samples from a copula family
  92
  93 function x = copularnd (family, theta, nu, n)
  94
  95   # Check arguments
  96   if (nargin < 2)
  97     print_usage ();
  98   endif
  99
 100   if (! ischar (family))
 101     error ("copularnd: family must be one of 'Gaussian', 't', and 'Clayton'");
 102   endif
 103
 104   lower_family = lower (family);
 105
 106   # Check family and copula parameters
 107   switch (lower_family)
 108
 109     case {"gaussian"}
 110       # Gaussian family
 111       if (isscalar (theta))
 112         # Expand a scalar to a correlation matrix
 113         theta = [1, theta; theta, 1];
 114       endif
 115       if (! ismatrix (theta) || any (diag (theta) != 1) || any (any (theta != theta')) || min (eig (theta)) <= 0)
 116         error ("copularnd: theta must be a correlation matrix");
 117       endif
 118       if (nargin > 3)
 119         d = n;
 120         if (! isscalar (d) || d != size (theta, 1))
 121           error ("copularnd: d must correspond to dimension of theta");
 122         endif
 123       else
 124         d = size (theta, 1);
 125       endif
 126       if (nargin < 3)
 127         n = 1;
 128       else
 129         n = nu;
 130         if (! isscalar (n) || (n < 0) || round (n) != n)
 131           error ("copularnd: n must be a non-negative integer");
 132         endif
 133       endif
 134
 135     case {"t"}
 136       # Student's t family
 137       if (nargin < 3)
 138         print_usage ();
 139       endif
 140       if (isscalar (theta))
 141         # Expand a scalar to a correlation matrix
 142         theta = [1, theta; theta, 1];
 143       endif
 144       if (! ismatrix (theta) || any (diag (theta) != 1) || any (any (theta != theta')) || min (eig (theta)) <= 0)
 145         error ("copularnd: theta must be a correlation matrix");
 146       endif
 147       if (! isscalar (nu) && (! isvector (nu) || length (nu) != n))
 148         error ("copularnd: nu must be a vector with the same number of rows as x or be scalar");
 149       endif
 150       nu = nu(:);
 151       if (nargin < 4)
 152         n = 1;
 153       else
 154         if (! isscalar (n) || (n < 0) || round (n) != n)
 155           error ("copularnd: n must be a non-negative integer");
 156         endif
 157       endif
 158
 159     case {"clayton"}
 160       # Archimedian one parameter family
 161       if (nargin < 4)
 162         # Default is bivariate
 163         d = 2;
 164       else
 165         d = n;
 166         if (! isscalar (d) || (d < 2) || round (d) != d)
 167           error ("copularnd: d must be an integer greater than 1");
 168         endif
 169       endif
 170       if (nargin < 3)
 171         # Default is one sample
 172         n = 1;
 173       else
 174         n = nu;
 175         if (! isscalar (n) || (n < 0) || round (n) != n)
 176           error ("copularnd: n must be a non-negative integer");
 177         endif
 178       endif
 179       if (! isvector (theta) || (! isscalar (theta) && size (theta, 1) != n))
 180         error ("copularnd: theta must be a column vector with the number of rows equal to n or be scalar");
 181       endif
 182       if (n > 1 && isscalar (theta))
 183         theta = repmat (theta, n, 1);
 184       endif
 185
 186     otherwise
 187       error ("copularnd: unknown copula family '%s'", family);
 188
 189   endswitch
 190
 191   if (n == 0)
 192     # Input is empty
 193     x = zeros (0, d);
 194   else
 195
 196     # Draw random samples according to family
 197     switch (lower_family)
 198
 199       case {"gaussian"}
 200         # The Gaussian family
 201         x = normcdf (mvnrnd (zeros (1, d), theta, n), 0, 1);
 202         # No parameter bounds check
 203         k = [];
 204
 205       case {"t"}
 206         # The Student's t family
 207         x = tcdf (mvtrnd (theta, nu, n), nu);
 208         # No parameter bounds check
 209         k = [];
 210
 211       case {"clayton"}
 212         # The Clayton family
 213         u = rand (n, d);
 214         if (d == 2)
 215           x = zeros (n, 2);
 216           # Conditional distribution method for the bivariate case which also
 217           # works for theta < 0
 218           x(:, 1) = u(:, 1);
 219           x(:, 2) = (1 + u(:, 1) .^ (-theta) .* (u(:, 2) .^ (-theta ./ (1 + theta)) - 1)) .^ (-1 ./ theta);
 220         else
 221           # Apply the algorithm by Marshall and Olkin:
 222           # Frailty distribution for Clayton copula is gamma
 223           y = randg (1 ./ theta, n, 1);
 224           x = (1 - log (u) ./ repmat (y, 1, d)) .^ (-1 ./ repmat (theta, 1, d));
 225         endif
 226         k = find (theta == 0);
 227         if (any (k))
 228           # Produkt copula at columns k
 229           x(k, :) = u(k, :);
 230         endif
 231         # Continue argument check
 232         if (d == 2)
 233           k = find (! (theta >= -1) | ! (theta < inf));
 234         else
 235           k = find (! (theta >= 0) | ! (theta < inf));
 236         endif
 237
 238     endswitch
 239
 240     # Out of bounds parameters
 241     if (any (k))
 242       x(k, :) = NaN;
 243     endif
 244
 245   endif
 246
 247 endfunction
 248
 249 %!test
 250 %! theta = 0.5;
 251 %! x = copularnd ("Gaussian", theta);
 252 %! assert (size (x), [1, 2]);
 253 %! assert (all ((x >= 0) & (x <= 1)));
 254
 255 %!test
 256 %! theta = 0.5;
 257 %! nu = 2;
 258 %! x = copularnd ("t", theta, nu);
 259 %! assert (size (x), [1, 2]);
 260 %! assert (all ((x >= 0) & (x <= 1)));
 261
 262 %!test
 263 %! theta = 0.5;
 264 %! x = copularnd ("Clayton", theta);
 265 %! assert (size (x), [1, 2]);
 266 %! assert (all ((x >= 0) & (x <= 1)));
 267
 268 %!test
 269 %! theta = 0.5;
 270 %! n = 2;
 271 %! x = copularnd ("Clayton", theta, n);
 272 %! assert (size (x), [n, 2]);
 273 %! assert (all ((x >= 0) & (x <= 1)));
 274
 275 %!test
 276 %! theta = [1; 2];
 277 %! n = 2;
 278 %! d = 3;
 279 %! x = copularnd ("Clayton", theta, n, d);
 280 %! assert (size (x), [n, d]);
 281 %! assert (all ((x >= 0) & (x <= 1)));