function [ci3] = bino_wilson(n, x, aa)

%   Wilson's 2-sided shrinkage interval of binomial

p_h = x./n;      za = norminv(1-aa/2,0,1); 
tm1 = p_h + za^2 / (2*n);
tm2 = za*sqrt((p_h.*(1-p_h)+za^2/(4*n))/n);
tm3 = 1 + za^2 / n;

p1 = (tm1 - tm2)/tm3;
p2 = (tm1 + tm2)/tm3;

if sum(p1 < 0) >=1
    for k = 1 : length(x)
        if p1(k,1) < 0
            p1(k,1) = 0;
        else
        end
    end
else if sum(p2 > 1) >= 1
        for k = 1 : length(x)
            if p1(k,1) > 1
                p1(k,1) = 1;
            else
            end
        end
    else
    end
end

ci3 = [p1 p2];