function [p_upr] = bino_upr(n, x, aa)


%   upper limit of confidence interval for 100(1-aa)%

tol = 10^(-10);     p_upr = [];

for k = 1 : length(x)
    p10 = 1 - eps;      p20 = eps;      p = (p10 + p20)/2;
    lwr = binocdf(x(k,1), n, p);

    while aa < lwr | aa - lwr > tol
        if aa > lwr
            p10 = p;    p = (p10 + p20)/2;
        else
            p20 = p;    p = (p10 + p20)/2;
        end
        lwr = binocdf(x(k,1), n, p);
        if p < eps | p > 1-eps
            break;
        else
        end
    end
    p_upr = [p_upr; p];
end