########################################################################
##
## Copyright (C) 2017-2021 The Octave Project Developers
##
## See the file COPYRIGHT.md in the top-level directory of this
## distribution or .
##
## This file is part of Octave.
##
## Octave is free software: you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING. If not, see
## .
##
########################################################################
## -*- texinfo -*-
## @deftypefn {} {} betaincinv (@var{y}, @var{a}, @var{b})
## @deftypefnx {} {} betaincinv (@var{y}, @var{a}, @var{b}, "lower")
## @deftypefnx {} {} betaincinv (@var{y}, @var{a}, @var{b}, "upper")
## Compute the inverse of the normalized incomplete beta function.
##
## The normalized incomplete beta function is defined as
## @tex
## $$
## I_x (a, b) = {1 \over {B(a,b)}} \displaystyle{\int_0^x t^{a-1} (1-t)^{b-1} dt}
## $$
## @end tex
## @ifnottex
##
## @example
## @group
## x
## /
## |
## I_x (a, b) = | t^(a-1) (1-t)^(b-1) dt
## |
## /
## 0
## @end group
## @end example
##
## @end ifnottex
##
## If two inputs are scalar, then @code{betaincinv (@var{y}, @var{a}, @var{b})}
## is returned for each of the other inputs.
##
## If two or more inputs are not scalar, the sizes of them must agree, and
## @code{betaincinv} is applied element-by-element.
##
## The variable @var{y} must be in the interval [0,1], while @var{a} and
## @var{b} must be real and strictly positive.
##
## By default, @var{tail} is @qcode{"lower"} and the inverse of the incomplete
## beta function integrated from 0 to @var{x} is computed. If @var{tail} is
## @qcode{"upper"} then the complementary function integrated from @var{x} to 1
## is inverted.
##
## The function is computed by standard Newton's method, by solving
## @tex
## $$
## y - I_x (a, b) = 0
## $$
## @end tex
## @ifnottex
##
## @example
## @var{y} - betainc (@var{x}, @var{a}, @var{b}) = 0
## @end example
##
## @end ifnottex
##
## @seealso{betainc, beta, betaln}
## @end deftypefn
function x = biinv (y, a, b, tail = "lower")
if (nargin < 3)
print_usage ();
endif
[err, y, a, b] = common_size (y, a, b);
if (err > 0)
error ("betaincinv: Y, A, and B must be of common size or scalars");
endif
if (! (isfloat (y) && isfloat (a) && isfloat (b)
&& isreal (y) && isreal (a) && isreal (b)))
error ("betaincinv: Y, A, and B must be real, floating point values");
endif
## Remember original shape of data, but convert to column vector for calcs.
orig_sz = size (y);
y = y(:);
a = a(:);
b = b(:);
if (any ((y < 0) | (y > 1)))
error ("betaincinv: Y must be in the range [0, 1]");
endif
if (any (a <= 0))
error ("betaincinv: A must be strictly positive");
endif
if (any (b <= 0))
error ("betaincinv: B must be strictly positive");
endif
## If any of the arguments is single then the output should be as well.
if (isa (y, "single") || isa (a, "single") || isa (b, "single"))
y = single (y);
a = single (a);
b = single (b);
endif
## Initialize output array
x = zeros (size (y), class (y));
if (strcmpi (tail, "lower"))
ys = y;
elseif (strcmpi (tail, "upper"))
ys = 1 - y; #only for computation of initial points, no loss of accuracy
else
error ("betaincinv: invalid value for TAIL");
endif
## If (a-1)*(b-1)>0, F has a point of inflection at x=(a-1)/(a+b-2).
## In this case, it is convex on (0,x) and concave on (x,1) if a>1, otherwise
## it is the other way round. If (a-1)*(b-1)<=0, there is no point of
## inflection, and it is everywhere convex for a>1 and concave otherwise.
## We thus choose our starting x for the Newton iterations so that we stay
## within a region of constant sign of curvature and on the correct side of
## the eventual solution, guaranteeing convergence. Curvatures above are to
## be understood under the condition tail=="lower".
x_i = y_i = x * 0;
i_swap = x < 0; #false
## Have point of inflection
ind = find ((a - 1) .* (b - 1) > 0);
if any (ind)
x_i(ind) = (a(ind) - 1) ./ (a(ind) + b(ind) - 2);
y_i(ind) = betainc (x_i(ind), a(ind), b(ind));
end
## Converge outwards
ind2 = a(ind) > 1;
if any (ind2)
x(ind(ind2)) = x_i(ind(ind2));
end
## Converge inwards
## To the left of inflection point
tmpind = ind(find ((a(ind) <= 1) & (y_i(ind) >= ys(ind))));
if length (tmpind) > 0
x(tmpind) = (ys(tmpind) ./ y_i(tmpind)) .^ (1 ./ a(tmpind)) .* x_i(tmpind);
end
## To the right of inflection point
tmpind = ind(find ((a(ind) <= 1) & (y_i(ind) < ys(ind))));
if length (tmpind) > 0
x(tmpind) = 1 - ((1 - ys(tmpind)) ./ (1 - y_i(tmpind))) .^ ...
(1 ./ b(tmpind)) .* (1 - x_i(tmpind));
end
## Have no point of inflection
ind = find ((a - 1) .* (b - 1) <= 0);
## Negative curvature
tmpind = ind(find (a(ind) < 1));
if length (tmpind) > 0
x(tmpind) = (ys(tmpind) .* beta (a(tmpind), b(tmpind)) .* ...
a(tmpind)) .^ (1 ./ a(tmpind));
end
## Positive curvature
tmpind = ind(find (a(ind) >= 1));
if length (tmpind) > 0
x(tmpind) = 1 - ((1 - ys(tmpind)) .* beta (a(tmpind), b(tmpind)) .* ...
b(tmpind)) .^ (1 ./ b(tmpind));
end
if (strcmpi (tail, "lower"))
x(y == 0) = 0;
x(y == 1) = 1;
F = @(x, a, b, y) y - betainc (x, a, b);
JF = @(x, a, b, Bln) -exp ((a-1) .* log (x) + (b-1) .* log1p (-x) - Bln);
elseif (strcmpi (tail, "upper"))
x(y == 0) = 1;
x(y == 1) = 0;
F = @(x, a, b, y) y - betainc (x, a, b, "upper");
JF = @(x, a, b, Bln) exp ((a-1) .* log (x) + (b-1) .* log1p (-x) - Bln);
endif
x = newton_method (F, JF, x, a, b, y);
## Restore original shape
x = reshape (x, orig_sz);
endfunction
function x = newton_method (F, JF, x, a, b, y);
Bln = betaln (a, b);
## Special values have been already computed.
todo = find((y != 0) & (y != 1));
step = -F (x(todo), a(todo), b(todo), y(todo)) ./ ...
JF (x(todo), a(todo), b(todo), Bln(todo));
df = (x(todo) + step) - x(todo);
x(todo) += step;
ind = df != 0;
todo = todo(ind);
df = df(ind);
iter = 1;
while (length (todo) > 0)
printf ("iter: %d, # todo: %d, x: %.15g, max (df): %.15g\n", iter++, numel (todo), x, max (abs (df(:))));
step = -F (x(todo), a(todo), b(todo), y(todo)) ./ ...
JF (x(todo), a(todo), b(todo), Bln(todo));
df_new = (x(todo) + step) - x(todo);
x(todo) += step;
ind = df_new .* df > 0;
todo = todo(ind);
df = df_new(ind);
endwhile
endfunction
%!test
%! x = linspace (0.1, 0.9, 11);
%! a = [2, 3, 4];
%! [x,a,b] = ndgrid (x,a,a);
%! xx = betaincinv (betainc (x, a, b), a, b);
%! assert (xx, x, 3e-15);
%!test
%! x = linspace (0.1, 0.9, 11);
%! a = [2, 3, 4];
%! [x,a,b] = ndgrid (x,a,a);
%! xx = betaincinv (betainc (x, a, b, "upper"), a, b, "upper");
%! assert (xx, x, 3e-15);
%!test
%! x = linspace (0.1, 0.9, 11);
%! a = [0.1:0.1:1];
%! [x,a,b] = ndgrid (x,a,a);
%! xx = betaincinv (betainc (x, a, b), a, b);
%! assert (xx, x, 5e-15);
%!test
%! x = linspace (0.1, 0.9, 11);
%! a = [0.1:0.1:1];
%! [x,a,b] = ndgrid (x,a,a);
%! xx = betaincinv (betainc (x, a, b, "upper"), a, b, "upper");
%! assert (xx, x, 5e-15);
## Test the conservation of the input class
%!assert (class (betaincinv (0.5, 1, 1)), "double")
%!assert (class (betaincinv (single (0.5), 1, 1)), "single")
%!assert (class (betaincinv (0.5, single (1), 1)), "single")
%!assert (class (betaincinv (0.5, 1, single (1))), "single")
%!assert <*60528> (betaincinv (1e-6, 1, 3), 3.3333344444450657e-07, 5*eps)
%!assert <*60528> (betaincinv (1-1e-6, 3, 1), 0.999999666666556, 5*eps)
## Test input validation
%!error betaincinv ()
%!error betaincinv (1)
%!error betaincinv (1,2)
%!error
%! betaincinv (ones (2,2), ones (1,2), 1);
%!error betaincinv ('a', 1, 2)
%!error betaincinv (0, int8 (1), 1)
%!error betaincinv (0, 1, true)
%!error betaincinv (0.5i, 1, 2)
%!error betaincinv (0, 1i, 1)
%!error betaincinv (0, 1, 1i)
%!error betaincinv (-0.1,1,1)
%!error betaincinv (1.1,1,1)
%!error
%! y = ones (1, 1, 2);
%! y(1,1,2) = -1;
%! betaincinv (y,1,1);
%!error betaincinv (0.5,0,1)
%!error
%! a = ones (1, 1, 2);
%! a(1,1,2) = 0;
%! betaincinv (1,a,1);
%!error betaincinv (0.5,1,0)
%!error
%! b = ones (1, 1, 2);
%! b(1,1,2) = 0;
%! betaincinv (1,1,b);
%!error betaincinv (1,2,3, "foobar")