custom_betainv.m

toolbox of BVQX, This is the access between BV and matlab. It will help you to analysis data from BV

function binv = custom_betainv(p, a, b)
% custom_betainv  - inverse function of the beta distribution
%
% FORMAT:       binv = custom_betainv(p, a, b)
%
% Input fields:
%
%       p           values for which the inverse beta CDF is computed
%       a, b        beta function parameters
%
% Output fields:
%
%       binv        inverse beta function for (p, a, b)

% Version:  v0.6f
% Build:    7070510
% Date:     Jul-05 2007, 10:07 AM CEST
% Author:   KH <[e-mail]> (OCTAVE)

% check arguments
if nargin < 3 || ...
   ~isa(p, 'double') || ...
   ~isa(a, 'double') || ...
   ~isa(b, 'double') || ...
    isempty(a) || ...
    isempty(b)
    error( ...
        'BVQXtools:BadArgument', ...
        'Missing or invalid argument given.' ...
    );
end
if isempty(p)
    binv = [];
    return;
end
osize = size(p);
if numel(p) == 1
    if numel(a) > 1
        osize = size(a);
    elseif numel(b) > 1
        osize = size(b);
    end
end
p = p(:);
a = a(:);
b = b(:);

% build output
binv = zeros(osize);

% special cases
binv(p < 0 | p > 1 | ~(a > 0) | ~(b > 0) | isnan(p)) = NaN;
binv(p == 1 & a > 0 & b > 0) = 1;

% find normal case
k = find(p > 0 & p < 1 & a > 0 & b > 0);
if ~isempty(k)
    if numel(a) > 1
        a = a(k);
    end
    if numel(b) > 1
        b = b(k);
    end
    if numel(a) ~= 1 || ...
        numel(b) ~= 1
    	y = a ./ (a + b);
    else
        y = a ./ (a + b) * ones(numel(k), 1);
    end
    p = p(k);
    
    % prepare loop
    y(y < eps) = sqrt(eps);
    y(y > (1 - eps)) = 1 - sqrt(eps);

    y_old = y;
    for lc = 1:10000
        h = (custom_betacdf(y_old, a, b) - p) ./ custom_betapdf(y_old, a, b);
        y_new = y_old - h;
        ind = find(y_new <= eps);
        if ~isempty(ind)
            y_new(ind) = y_old(ind) ./ 10;
        end
        ind = find(y_new >= (1 - eps));
        if ~isempty(ind)
            y_new(ind) = 1 - (1 - y_old(ind)) ./ 10;
        end
        h = y_old - y_new;
        if (max(abs(h)) < sqrt(eps))
            break;
        end
        y_old = y_new;
    end
    binv(k) = y_new;
end

% reshape
binv = reshape(binv, osize);