octave_packages/statistics-1.1.3/copulapdf.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} =} copulapdf (@var{family}, @var{x}, @var{theta})
  18 ## Compute the probability density function of a copula family.
  19 ##
  20 ## @subheading Arguments
  21 ##
  22 ## @itemize @bullet
  23 ## @item
  24 ## @var{family} is the copula family name. Currently, @var{family} can
  25 ## be @code{'Clayton'} for the Clayton family, @code{'Gumbel'} for the
  26 ## Gumbel-Hougaard family, @code{'Frank'} for the Frank family, or
  27 ## @code{'AMH'} for the Ali-Mikhail-Haq family.
  28 ##
  29 ## @item
  30 ## @var{x} is the support where each row corresponds to an observation.
  31 ##
  32 ## @item
  33 ## @var{theta} is the parameter of the copula. The elements of
  34 ## @var{theta} must be greater than or equal to @code{-1} for the
  35 ## Clayton family, greater than or equal to @code{1} for the
  36 ## Gumbel-Hougaard family, arbitrary for the Frank family, and greater
  37 ## than or equal to @code{-1} and lower than @code{1} for the
  38 ## Ali-Mikhail-Haq family. Moreover, @var{theta} must be non-negative
  39 ## for dimensions greater than @code{2}. @var{theta} must be a column
  40 ## vector with the same number of rows as @var{x} or be scalar.
  41 ## @end itemize
  42 ##
  43 ## @subheading Return values
  44 ##
  45 ## @itemize @bullet
  46 ## @item
  47 ## @var{p} is the probability density of the copula at each row of
  48 ## @var{x} and corresponding parameter @var{theta}.
  49 ## @end itemize
  50 ##
  51 ## @subheading Examples
  52 ##
  53 ## @example
  54 ## @group
  55 ## x = [0.2:0.2:0.6; 0.2:0.2:0.6];
  56 ## theta = [1; 2];
  57 ## p = copulapdf ("Clayton", x, theta)
  58 ## @end group
  59 ##
  60 ## @group
  61 ## p = copulapdf ("Gumbel", x, 2)
  62 ## @end group
  63 ## @end example
  64 ##
  65 ## @subheading References
  66 ##
  67 ## @enumerate
  68 ## @item
  69 ## Roger B. Nelsen. @cite{An Introduction to Copulas}. Springer,
  70 ## New York, second edition, 2006.
  71 ## @end enumerate
  72 ## @end deftypefn
  73
  74 ## Author: Arno Onken <asnelt@asnelt.org>
  75 ## Description: PDF of a copula family
  76
  77 function p = copulapdf (family, x, theta)
  78
  79   # Check arguments
  80   if (nargin != 3)
  81     print_usage ();
  82   endif
  83
  84   if (! ischar (family))
  85     error ("copulapdf: family must be one of 'Clayton', 'Gumbel', 'Frank', and 'AMH'");
  86   endif
  87
  88   if (! isempty (x) && ! ismatrix (x))
  89     error ("copulapdf: x must be a numeric matrix");
  90   endif
  91
  92   [n, d] = size (x);
  93
  94   if (! isvector (theta) || (! isscalar (theta) && size (theta, 1) != n))
  95     error ("copulapdf: theta must be a column vector with the same number of rows as x or be scalar");
  96   endif
  97
  98   if (n == 0)
  99     # Input is empty
 100     p = zeros (0, 1);
 101   else
 102     if (n > 1 && isscalar (theta))
 103       theta = repmat (theta, n, 1);
 104     endif
 105
 106     # Truncate input to unit hypercube
 107     x(x < 0) = 0;
 108     x(x > 1) = 1;
 109
 110     # Compute the cumulative distribution function according to family
 111     lowerarg = lower (family);
 112
 113     if (strcmp (lowerarg, "clayton"))
 114       # The Clayton family
 115       log_cdf = -log (max (sum (x .^ (repmat (-theta, 1, d)), 2) - d + 1, 0)) ./ theta;
 116       p = prod (repmat (theta, 1, d) .* repmat (0:(d - 1), n, 1) + 1, 2) .* exp ((1 + theta .* d) .* log_cdf - (theta + 1) .* sum (log (x), 2));
 117       # Product copula at columns where theta == 0
 118       k = find (theta == 0);
 119       if (any (k))
 120         p(k) = 1;
 121       endif
 122       # Check theta
 123       if (d > 2)
 124         k = find (! (theta >= 0) | ! (theta < inf));
 125       else
 126         k = find (! (theta >= -1) | ! (theta < inf));
 127       endif
 128     elseif (strcmp (lowerarg, "gumbel"))
 129       # The Gumbel-Hougaard family
 130       g = sum ((-log (x)) .^ repmat (theta, 1, d), 2);
 131       c = exp (-g .^ (1 ./ theta));
 132       p = ((prod (-log (x), 2)) .^ (theta - 1)) ./ prod (x, 2) .* c .* (g .^ (2 ./ theta - 2) + (theta - 1) .* g .^ (1 ./ theta - 2));
 133       # Check theta
 134       k = find (! (theta >= 1) | ! (theta < inf));
 135     elseif (strcmp (lowerarg, "frank"))
 136       # The Frank family
 137       if (d != 2)
 138         error ("copulapdf: Frank copula PDF implemented as bivariate only");
 139       endif
 140       p = (theta .* exp (theta .* (1 + sum (x, 2))) .* (exp (theta) - 1))./ (exp (theta) - exp (theta + theta .* x(:, 1)) + exp (theta .* sum (x, 2)) - exp (theta + theta .* x(:, 2))) .^ 2;
 141       # Product copula at columns where theta == 0
 142       k = find (theta == 0);
 143       if (any (k))
 144         p(k) = 1;
 145       endif
 146       # Check theta
 147       k = find (! (theta > -inf) | ! (theta < inf));
 148     elseif (strcmp (lowerarg, "amh"))
 149       # The Ali-Mikhail-Haq family
 150       if (d != 2)
 151         error ("copulapdf: Ali-Mikhail-Haq copula PDF implemented as bivariate only");
 152       endif
 153       z = theta .* prod (x - 1, 2) - 1;
 154       p = (theta .* (1 - sum (x, 2) - prod (x, 2) - z) - 1) ./ (z .^ 3);
 155       # Check theta
 156       k = find (! (theta >= -1) | ! (theta < 1));
 157     else
 158       error ("copulapdf: unknown copula family '%s'", family);
 159     endif
 160
 161     if (any (k))
 162       p(k) = NaN;
 163     endif
 164
 165   endif
 166
 167 endfunction
 168
 169 %!test
 170 %! x = [0.2:0.2:0.6; 0.2:0.2:0.6];
 171 %! theta = [1; 2];
 172 %! p = copulapdf ("Clayton", x, theta);
 173 %! expected_p = [0.9872; 0.7295];
 174 %! assert (p, expected_p, 0.001);
 175
 176 %!test
 177 %! x = [0.2:0.2:0.6; 0.2:0.2:0.6];
 178 %! p = copulapdf ("Gumbel", x, 2);
 179 %! expected_p = [0.9468; 0.9468];
 180 %! assert (p, expected_p, 0.001);
 181
 182 %!test
 183 %! x = [0.2, 0.6; 0.2, 0.6];
 184 %! theta = [1; 2];
 185 %! p = copulapdf ("Frank", x, theta);
 186 %! expected_p = [0.9378; 0.8678];
 187 %! assert (p, expected_p, 0.001);
 188
 189 %!test
 190 %! x = [0.2, 0.6; 0.2, 0.6];
 191 %! theta = [0.3; 0.7];
 192 %! p = copulapdf ("AMH", x, theta);
 193 %! expected_p = [0.9540; 0.8577];
 194 %! assert (p, expected_p, 0.001);