octave_packages/statistics-1.1.3/mnpdf.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{y} =} mnpdf (@var{x}, @var{p})
  18 ## Compute the probability density function of the multinomial distribution.
  19 ##
  20 ## @subheading Arguments
  21 ##
  22 ## @itemize @bullet
  23 ## @item
  24 ## @var{x} is vector with a single sample of a multinomial distribution with
  25 ## parameter @var{p} or a matrix of random samples from multinomial
  26 ## distributions. In the latter case, each row of @var{x} is a sample from a
  27 ## multinomial distribution with the corresponding row of @var{p} being its
  28 ## parameter.
  29 ##
  30 ## @item
  31 ## @var{p} is a vector with the probabilities of the categories or a matrix
  32 ## with each row containing the probabilities of a multinomial sample.
  33 ## @end itemize
  34 ##
  35 ## @subheading Return values
  36 ##
  37 ## @itemize @bullet
  38 ## @item
  39 ## @var{y} is a vector of probabilites of the random samples @var{x} from the
  40 ## multinomial distribution with corresponding parameter @var{p}. The parameter
  41 ## @var{n} of the multinomial distribution is the sum of the elements of each
  42 ## row of @var{x}. The length of @var{y} is the number of columns of @var{x}.
  43 ## If a row of @var{p} does not sum to @code{1}, then the corresponding element
  44 ## of @var{y} will be @code{NaN}.
  45 ## @end itemize
  46 ##
  47 ## @subheading Examples
  48 ##
  49 ## @example
  50 ## @group
  51 ## x = [1, 4, 2];
  52 ## p = [0.2, 0.5, 0.3];
  53 ## y = mnpdf (x, p);
  54 ## @end group
  55 ##
  56 ## @group
  57 ## x = [1, 4, 2; 1, 0, 9];
  58 ## p = [0.2, 0.5, 0.3; 0.1, 0.1, 0.8];
  59 ## y = mnpdf (x, p);
  60 ## @end group
  61 ## @end example
  62 ##
  63 ## @subheading References
  64 ##
  65 ## @enumerate
  66 ## @item
  67 ## Wendy L. Martinez and Angel R. Martinez. @cite{Computational Statistics
  68 ## Handbook with MATLAB}. Appendix E, pages 547-557, Chapman & Hall/CRC, 2001.
  69 ##
  70 ## @item
  71 ## Merran Evans, Nicholas Hastings and Brian Peacock. @cite{Statistical
  72 ## Distributions}. pages 134-136, Wiley, New York, third edition, 2000.
  73 ## @end enumerate
  74 ## @end deftypefn
  75
  76 ## Author: Arno Onken <asnelt@asnelt.org>
  77 ## Description: PDF of the multinomial distribution
  78
  79 function y = mnpdf (x, p)
  80
  81   # Check arguments
  82   if (nargin != 2)
  83     print_usage ();
  84   endif
  85
  86   if (! ismatrix (x) || any (x(:) < 0 | round (x(:) != x(:))))
  87     error ("mnpdf: x must be a matrix of non-negative integer values");
  88   endif
  89   if (! ismatrix (p) || any (p(:) < 0))
  90     error ("mnpdf: p must be a non-empty matrix with rows of probabilities");
  91   endif
  92
  93   # Adjust input sizes
  94   if (! isvector (x) || ! isvector (p))
  95     if (isvector (x))
  96       x = x(:)';
  97     endif
  98     if (isvector (p))
  99       p = p(:)';
 100     endif
 101     if (size (x, 1) == 1 && size (p, 1) > 1)
 102       x = repmat (x, size (p, 1), 1);
 103     elseif (size (x, 1) > 1 && size (p, 1) == 1)
 104       p = repmat (p, size (x, 1), 1);
 105     endif
 106   endif
 107   # Continue argument check
 108   if (any (size (x) != size (p)))
 109     error ("mnpdf: x and p must have compatible sizes");
 110   endif
 111
 112   # Count total number of elements of each multinomial sample
 113   n = sum (x, 2);
 114   # Compute probability density function of the multinomial distribution
 115   t = x .* log (p);
 116   t(x == 0) = 0;
 117   y = exp (gammaln (n+1) - sum (gammaln (x+1), 2) + sum (t, 2));
 118   # Set invalid rows to NaN
 119   k = (abs (sum (p, 2) - 1) > 1e-6);
 120   y(k) = NaN;
 121
 122 endfunction
 123
 124 %!test
 125 %! x = [1, 4, 2];
 126 %! p = [0.2, 0.5, 0.3];
 127 %! y = mnpdf (x, p);
 128 %! assert (y, 0.11812, 0.001);
 129
 130 %!test
 131 %! x = [1, 4, 2; 1, 0, 9];
 132 %! p = [0.2, 0.5, 0.3; 0.1, 0.1, 0.8];
 133 %! y = mnpdf (x, p);
 134 %! assert (y, [0.11812; 0.13422], 0.001);