######################################################################## ## ## 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 = biinv3 (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); printf ("iter: 0, x: %.15g\n", x); tol = eps (class (x)); 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)); x(todo) += step; ind = (step != 0); todo = todo(ind); iter = 1; while (length (todo) > 0) printf ("iter: %d, # todo: %d, x: %.15g, max (df): %.15g\n", iter++, numel (todo), x, max (abs (step(:)))); step = -F (x(todo), a(todo), b(todo), y(todo)) ./ ... JF (x(todo), a(todo), b(todo), Bln(todo)); x(todo) += step; ind = abs (step) > tol; todo = todo(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")