octave_packages/m/statistics/distributions/gaminv.m

   1 ## Copyright (C) 2012 Rik Wehbring
   2 ## Copyright (C) 1995-2012 Kurt Hornik
   3 ##
   4 ## This file is part of Octave.
   5 ##
   6 ## Octave is free software; you can redistribute it and/or modify it
   7 ## under the terms of the GNU General Public License as published by
   8 ## the Free Software Foundation; either version 3 of the License, or (at
   9 ## your option) any later version.
  10 ##
  11 ## Octave is distributed in the hope that it will be useful, but
  12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 ## General Public License for more details.
  15 ##
  16 ## You should have received a copy of the GNU General Public License
  17 ## along with Octave; see the file COPYING.  If not, see
  18 ## <http://www.gnu.org/licenses/>.
  19
  20 ## -*- texinfo -*-
  21 ## @deftypefn {Function File} {} gaminv (@var{x}, @var{a}, @var{b})
  22 ## For each element of @var{x}, compute the quantile (the inverse of
  23 ## the CDF) at @var{x} of the Gamma distribution with shape parameter
  24 ## @var{a} and scale @var{b}.
  25 ## @end deftypefn
  26
  27 ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
  28 ## Description: Quantile function of the Gamma distribution
  29
  30 function inv = gaminv (x, a, b)
  31
  32   if (nargin != 3)
  33     print_usage ();
  34   endif
  35
  36   if (!isscalar (a) || !isscalar (b))
  37     [retval, x, a, b] = common_size (x, a, b);
  38     if (retval > 0)
  39       error ("gaminv: X, A, and B must be of common size or scalars");
  40     endif
  41   endif
  42
  43   if (iscomplex (x) || iscomplex (a) || iscomplex (b))
  44     error ("gaminv: X, A, and B must not be complex");
  45   endif
  46
  47   if (isa (x, "single") || isa (a, "single") || isa (b, "single"))
  48     inv = zeros (size (x), "single");
  49   else
  50     inv = zeros (size (x));
  51   endif
  52
  53   k = ((x < 0) | (x > 1) | isnan (x)
  54        | !(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf));
  55   inv(k) = NaN;
  56
  57   k = (x == 1) & (a > 0) & (a < Inf) & (b > 0) & (b < Inf);
  58   inv(k) = Inf;
  59
  60   k = find ((x > 0) & (x < 1) & (a > 0) & (a < Inf) & (b > 0) & (b < Inf));
  61   if (any (k))
  62     if (!isscalar (a) || !isscalar (b))
  63       a = a(k);
  64       b = b(k);
  65       y = a .* b;
  66     else
  67       y = a * b * ones (size (k));
  68     endif
  69     x = x(k);
  70
  71     if (isa (x, "single"))
  72       myeps = eps ("single");
  73     else
  74       myeps = eps;
  75     endif
  76
  77     l = find (x < myeps);
  78     if (any (l))
  79       y(l) = sqrt (myeps) * ones (length (l), 1);
  80     endif
  81
  82     y_old = y;
  83     for i = 1 : 100
  84       h     = (gamcdf (y_old, a, b) - x) ./ gampdf (y_old, a, b);
  85       y_new = y_old - h;
  86       ind   = find (y_new <= myeps);
  87       if (any (ind))
  88         y_new (ind) = y_old (ind) / 10;
  89         h = y_old - y_new;
  90       endif
  91       if (max (abs (h)) < sqrt (myeps))
  92         break;
  93       endif
  94       y_old = y_new;
  95     endfor
  96
  97     inv(k) = y_new;
  98   endif
  99
 100 endfunction
 101
 102
 103 %!shared x
 104 %! x = [-1 0 0.63212055882855778 1 2];
 105 %!assert(gaminv (x, ones(1,5), ones(1,5)), [NaN 0 1 Inf NaN], eps);
 106 %!assert(gaminv (x, 1, ones(1,5)), [NaN 0 1 Inf NaN], eps);
 107 %!assert(gaminv (x, ones(1,5), 1), [NaN 0 1 Inf NaN], eps);
 108 %!assert(gaminv (x, [1 -Inf NaN Inf 1], 1), [NaN NaN NaN NaN NaN]);
 109 %!assert(gaminv (x, 1, [1 -Inf NaN Inf 1]), [NaN NaN NaN NaN NaN]);
 110 %!assert(gaminv ([x(1:2) NaN x(4:5)], 1, 1), [NaN 0 NaN Inf NaN]);
 111
 112 %% Test class of input preserved
 113 %!assert(gaminv ([x, NaN], 1, 1), [NaN 0 1 Inf NaN NaN], eps);
 114 %!assert(gaminv (single([x, NaN]), 1, 1), single([NaN 0 1 Inf NaN NaN]), eps("single"));
 115 %!assert(gaminv ([x, NaN], single(1), 1), single([NaN 0 1 Inf NaN NaN]), eps("single"));
 116 %!assert(gaminv ([x, NaN], 1, single(1)), single([NaN 0 1 Inf NaN NaN]), eps("single"));
 117
 118 %% Test input validation
 119 %!error gaminv ()
 120 %!error gaminv (1)
 121 %!error gaminv (1,2)
 122 %!error gaminv (1,2,3,4)
 123 %!error gaminv (ones(3),ones(2),ones(2))
 124 %!error gaminv (ones(2),ones(3),ones(2))
 125 %!error gaminv (ones(2),ones(2),ones(3))
 126 %!error gaminv (i, 2, 2)
 127 %!error gaminv (2, i, 2)
 128 %!error gaminv (2, 2, i)
 129