octave_packages/m/statistics/distributions/betainv.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} {} betainv (@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 Beta distribution with parameters @var{a}
  24 ## and @var{b}.
  25 ## @end deftypefn
  26
  27 ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
  28 ## Description: Quantile function of the Beta distribution
  29
  30 function inv = betainv (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 ("betainv: 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 ("betainv: 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) | !(a > 0) | !(b > 0) | isnan (x);
  54   inv(k) = NaN;
  55
  56   k = (x == 1) & (a > 0) & (b > 0);
  57   inv(k) = 1;
  58
  59   k = find ((x > 0) & (x < 1) & (a > 0) & (b > 0));
  60   if (any (k))
  61     if (!isscalar (a) || !isscalar (b))
  62       a = a(k);
  63       b = b(k);
  64       y = a ./ (a + b);
  65     else
  66       y = a / (a + b) * ones (size (k));
  67     endif
  68     x = x(k);
  69
  70     if (isa (y, "single"))
  71       myeps = eps ("single");
  72     else
  73       myeps = eps;
  74     endif
  75
  76     l = find (y < myeps);
  77     if (any (l))
  78       y(l) = sqrt (myeps) * ones (length (l), 1);
  79     endif
  80     l = find (y > 1 - myeps);
  81     if (any (l))
  82       y(l) = 1 - sqrt (myeps) * ones (length (l), 1);
  83     endif
  84
  85     y_old = y;
  86     for i = 1 : 10000
  87       h     = (betacdf (y_old, a, b) - x) ./ betapdf (y_old, a, b);
  88       y_new = y_old - h;
  89       ind   = find (y_new <= myeps);
  90       if (any (ind))
  91         y_new (ind) = y_old (ind) / 10;
  92       endif
  93       ind = find (y_new >= 1 - myeps);
  94       if (any (ind))
  95         y_new (ind) = 1 - (1 - y_old (ind)) / 10;
  96       endif
  97       h = y_old - y_new;
  98       if (max (abs (h)) < sqrt (myeps))
  99         break;
 100       endif
 101       y_old = y_new;
 102     endfor
 103
 104     inv(k) = y_new;
 105   endif
 106
 107 endfunction
 108
 109
 110 %!shared x
 111 %! x = [-1 0 0.75 1 2];
 112 %!assert(betainv (x, ones(1,5), 2*ones(1,5)), [NaN 0 0.5 1 NaN]);
 113 %!assert(betainv (x, 1, 2*ones(1,5)), [NaN 0 0.5 1 NaN]);
 114 %!assert(betainv (x, ones(1,5), 2), [NaN 0 0.5 1 NaN]);
 115 %!assert(betainv (x, [1 0 NaN 1 1], 2), [NaN NaN NaN 1 NaN]);
 116 %!assert(betainv (x, 1, 2*[1 0 NaN 1 1]), [NaN NaN NaN 1 NaN]);
 117 %!assert(betainv ([x(1:2) NaN x(4:5)], 1, 2), [NaN 0 NaN 1 NaN]);
 118
 119 %% Test class of input preserved
 120 %!assert(betainv ([x, NaN], 1, 2), [NaN 0 0.5 1 NaN NaN]);
 121 %!assert(betainv (single([x, NaN]), 1, 2), single([NaN 0 0.5 1 NaN NaN]));
 122 %!assert(betainv ([x, NaN], single(1), 2), single([NaN 0 0.5 1 NaN NaN]));
 123 %!assert(betainv ([x, NaN], 1, single(2)), single([NaN 0 0.5 1 NaN NaN]));
 124
 125 %% Test input validation
 126 %!error betainv ()
 127 %!error betainv (1)
 128 %!error betainv (1,2)
 129 %!error betainv (1,2,3,4)
 130 %!error betainv (ones(3),ones(2),ones(2))
 131 %!error betainv (ones(2),ones(3),ones(2))
 132 %!error betainv (ones(2),ones(2),ones(3))
 133 %!error betainv (i, 2, 2)
 134 %!error betainv (2, i, 2)
 135 %!error betainv (2, 2, i)
 136