unitproc/octave/invgaurnd.m

   1 ## Copyright (C) 2012 Rik Wehbring
   2 ## Copyright (C) 1995-2016 Kurt Hornik
   3 ## Copyright (C) 2021 Steven Baltakatei Sandoval
   4 ##
   5 ## This program is free software: you can redistribute it and/or
   6 ## modify it under the terms of the GNU General Public License as
   7 ## published by the Free Software Foundation, either version 3 of the
   8 ## License, or (at your option) any later version.
   9 ##
  10 ## This program is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with this program; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn  {} {} invgaurnd (@var{mu}, @var{lambda})
  21 ## @deftypefnx {} {} invgaurnd (@var{mu}, @var{lambda}, @var{r})
  22 ## @deftypefnx {} {} invgaurnd (@var{mu}, @var{lambda}, @var{r}, @var{c}, @dots{})
  23 ## @deftypefnx {} {} invgaurnd (@var{mu}, @var{lambda}, [@var{sz}])
  24 ## Return a matrix of random samples from the inverse gaussian distribution with
  25 ## parameters mean @var{mu} and shape parameter @var{lambda}.
  26 ##
  27 ## When called with a single size argument, return a square matrix with
  28 ## the dimension specified.  When called with more than one scalar argument the
  29 ## first two arguments are taken as the number of rows and columns and any
  30 ## further arguments specify additional matrix dimensions.  The size may also
  31 ## be specified with a vector of dimensions @var{sz}.
  32 ##
  33 ## If no size arguments are given then the result matrix is the common size of
  34 ## @var{mu} and @var{lambda}.
  35 ## @end deftypefn
  36
  37 ## Author: Steven Sandoval <baltakatei@gmail.com>
  38 ## Description: Random variates from the inverse gaussian distribution
  39
  40 function rnd = invgaurnd (mu, lambda, varargin)
  41
  42   if (nargin < 2)
  43     print_usage ();
  44   endif
  45
  46   if (! isscalar (mu) || ! isscalar (lambda))
  47     [retval, mu, lambda] = common_size (mu, lambda);
  48     if (retval > 0)
  49       error ("invgaurnd: MU and LAMBDA must be of common size or scalars");
  50     endif
  51   endif
  52
  53   if (iscomplex (mu) || iscomplex (lambda))
  54     error ("invgaurnd: MU and LAMBDA must not be complex");
  55   endif
  56
  57   if (nargin == 2)
  58     sz = size (mu);
  59   elseif (nargin == 3)
  60     if (isscalar (varargin{1}) && varargin{1} >= 0)
  61       sz = [varargin{1}, varargin{1}];
  62     elseif (isrow (varargin{1}) && all (varargin{1} >= 0))
  63       sz = varargin{1};
  64     else
  65       error ("invgaurnd: dimension vector must be row vector of non-negative integers");
  66     endif
  67   elseif (nargin > 3)
  68     if (any (cellfun (@(x) (! isscalar (x) || x < 0), varargin)))
  69       error ("invgaurnd: dimensions must be non-negative integers");
  70     endif
  71     sz = [varargin{:}];
  72   endif
  73
  74   if (! isscalar (mu) && ! isequal (size (mu), sz))
  75     error ("invgaurnd: MU and LAMBDA must be scalar or of size SZ");
  76   endif
  77
  78   if (isa (mu, "single") || isa (lambda, "single"))
  79     cls = "single";
  80   else
  81     cls = "double";
  82   endif;
  83
  84   # Convert mu and lambda from scalars into matrices since scalar multiplication used
  85   if isscalar (mu)
  86     mu = mu * ones(sz);
  87   endif;
  88   if isscalar (lambda)
  89     lambda = lambda * ones(sz);
  90   endif;
  91
  92   # Generate random variates
  93   # Ref/Attrib: Michael, John R., Generating Random Variates Using
  94   #   Transformations with Multiple Roots. The American Statistician, May 1976,
  95   #   Vol. 30, No. 2. https://doi.org/10.2307/2683801
  96   nu = randn(sz,cls);
  97   y = nu .** 2;
  98   x1 = mu;
  99   x2 = mu .** 2 .* y ./ (2 .* lambda);
 100   x3 = (- mu ./ (2 .* lambda)) .* sqrt(4 .* mu .* lambda .* y + mu .** 2 .* y .** 2);
 101   x = x1 + x2 + x3;
 102   z = rand(sz,cls);
 103   valTest1 = (mu ./ (mu + x)); # calculate test 1 value
 104   valTest2 = (mu ./ (mu + x)); # calculate test 2 value
 105   posTest1 = find(z <= valTest1); # perform test 1, save positions where test 1 true
 106   posTest2 = find(z > valTest2); # perform test 2, save positions where test 2 true
 107   ## indposTest1 = transpose(posTest1) # debug: list positions
 108   ## indposTest2 = transpose(posTest2) # debug: list positions
 109   ## indTest1 = z <= valTest1 # debug: show test 1 truth table
 110   ## indTest2 = z > valTest2 # debug: show test 2 truth table
 111   rnd = NaN(sz); # Initialize return array
 112   rnd(posTest1) = x(posTest1); # populate return matrix with corresp. elements of x that satisfy test 1
 113   rnd(posTest2) = (mu(posTest2) .** 2 ./ x(posTest2)); # populate return matrix with corresp. elements of x that satisfy test 2
 114   k = ! isfinite (mu) | !(lambda >= 0) | !(lambda < Inf); # store position matrix indicating which parts of output are invalid based on
 115                                                           #   elements of the matrices: mu, lambda.
 116   rnd(k) = NaN; # mark invalid positions of output matrix with NaN
 117
 118 endfunction
 119
 120
 121 %!assert (size (invgaurnd (1,2)), [1, 1])
 122 %!assert (size (invgaurnd (ones (2,1), 2)), [2, 1])
 123 %!assert (size (invgaurnd (ones (2,2), 2)), [2, 2])
 124 %!assert (size (invgaurnd (1, 2*ones (2,1))), [2, 1])
 125 %!assert (size (invgaurnd (1, 2*ones (2,2))), [2, 2])
 126 %!assert (size (invgaurnd (1, 2, 3)), [3, 3])
 127 %!assert (size (invgaurnd (1, 2, [4 1])), [4, 1])
 128 %!assert (size (invgaurnd (1, 2, 4, 1)), [4, 1])
 129
 130 ## Test class of input preserved
 131 %!assert (class (invgaurnd (1, 2)), "double")
 132 %!assert (class (invgaurnd (single (1), 2)), "single")
 133 %!assert (class (invgaurnd (single ([1 1]), 2)), "single")
 134 %!assert (class (invgaurnd (1, single (2))), "single")
 135 %!assert (class (invgaurnd (1, single ([2 2]))), "single")
 136
 137 ## Test input validation
 138 %!error invgaurnd ()
 139 %!error invgaurnd (1)
 140 %!error invgaurnd (ones (3), ones (2))
 141 %!error invgaurnd (ones (2), ones (3))
 142 %!error invgaurnd (i, 2)
 143 %!error invgaurnd (2, i)
 144 %!error invgaurnd (1,2, -1)
 145 %!error invgaurnd (1,2, ones (2))
 146 %!error invgaurnd (1, 2, [2 -1 2])
 147 %!error invgaurnd (1,2, 1, ones (2))
 148 %!error invgaurnd (1,2, 1, -1)
 149 %!error invgaurnd (ones (2,2), 2, 3)
 150 %!error invgaurnd (ones (2,2), 2, [3, 2])
 151 %!error invgaurnd (ones (2,2), 2, 2, 3)