

function compute(){
var outwin=document.theform.textout.value;

var k=0,n=0, cl=95;
k=parseInt(document.theform.successes.value);
n=parseInt(document.theform.trials.value);
cl=parseFloat(document.theform.conf.value);


if (isNaN(n))
        {
        alert("The number of trials is not a number.");
        return;
        }
if (n<=0)
        {
        alert("The number of trials must be a positive integer.");
        return;
        }

if (isNaN(k))
        {
        alert("The number of successes is not a number.");
        return;
        }
if ((k<0)||(k>n))
        {
        alert("The number of successes must be between zero "+
                "and the number of trials, inclusive.");
        return;
        }

if (isNaN(cl))
        {
        alert("The confidence level is not a number.");
        return;
        }
if ((cl<=0)||(cl>=100))
        {
        alert("The confidence level must be between "+
                "zero and 100, exclusive.");
        return;
        }


var out="";
out +="Here "+k+" successes were observed in "+
        n+" trials.\n"
out +="The MLE of the success probability is "+
        fix(k/n)+".\n";
out +="The lower endpoint of a one sided "+
        cl+"% confidence interval for the ";
out +="success probability is "+
        fix(bplower(n,k,(100-cl)/100))+".\n";
out +="The upper endpoint of a one sided "+
        cl+"% confidence interval for the ";
out +="success probability is "+
        fix(bpupper(n,k,(100-cl)/100))+".\n";
out +="The endpoints of a two sided "+
        cl+"% confidence interval for the ";
out +="success probability are "+
        fix(bplower(n,k,(100-cl)/200))+" and ";
out +=""+fix(bpupper(n,k,(100-cl)/200))+".\n\n";

document.theform.textout.value=out+outwin;
}



function fix(x)
{
return (Math.round(100000*x)/100000);
}


function loggamma(x) {
var a1=Math.log(2*Math.PI)/2;
var a2= 1/1680;
var a3= 1/1260;
var a4= 1/360;
var a5= 1/12;
var f,xx,z;
f=1;
if (x <=0) return(Number.NaN);
for(xx=x;xx < 7;xx++) f=f*xx;
f=-Math.log(f);
z=1/(xx*xx);
return(f+(xx-0.5)*Math.log(xx)-xx+a1+(((-a2*z+a3)*z-a4)*z+a5)/xx);
}


function dfbeta(x, p, q)
{

var one=1.0;
var zero=0.0;
var acu=0.0000001;
if ( (p<=zero) || (q<=zero) ) return(Number.NaN);
if (x<=zero) return(0.0);
if (x>=one) return(1.0);


var ai, beta, betain, psq, cx, ns, pp, qq, rx, temp, term, xx ;
var index;
betain=x;
psq=p+q;
cx=one-x;
if (p>= psq*x)
        {
        xx=x;
        pp=p;
        qq=q;
        index=0;
        }
else
        {
        xx=cx;
        cx=x;
        pp=q;
        qq=p;
        index=1;
        }
term=one;
ai=one;
betain=one;
ns=qq+cx*psq;


rx=xx/cx;
temp=qq-ai;
if (Math.floor(ns)==0)  rx=xx;
for(;;)
        {
        term=term*temp*rx/(pp+ai);
        betain=betain+term;
        temp=Math.abs(term);
        if  ( (temp <=acu) && (temp<=acu*betain) ) break;
        ai=ai+one;
        ns=ns-1;
        if (Math.floor(ns) >=0)
                {
                temp=qq-ai;
                if (Math.floor(ns) ==0) rx=xx;
                }
        else
                {
                temp=psq;
                psq=psq+one;
                }
        }

beta=loggamma(p)+loggamma(q)-loggamma(p+q);
betain=betain*Math.exp( pp*Math.log(xx) +(qq-one)*Math.log(cx) -beta)/pp;
if (index==1) betain=one-betain;
return(betain);
}


function ppbeta(x, p, q)
{

var     acu=1.0E-14, lower=0.0001, upper=0.9999, const1=2.30753;
var     const2=0.27061, const3=0.99229, const4=0.04481;
var beta=loggamma(p)+loggamma(q)-loggamma(p+q);
var     xinbta=x;
var  a, adj, g, gg, gx, h, pp, prev, qq, r, s, sq, t, tx, y, yprev, w ;
var index;
if ( (p<0) || (q<0) ) return(Number.NaN);
if (x<=0) return(0.0);
if (x>=1) return(1.0);


if (x <=0.5)
        {
        a=x;
        pp=p;
        qq=q;
        index=0;
        }
else
        {
        a=1.0-x;
        pp=q;
        qq=p;
        index=1;
        }


for(;;){ 
r=Math.sqrt(-Math.log(a*a));
y=r-(const1+const2*r)/(1+(const3+const4*r)*r);
if ( (pp>1) && (qq>1) )
        {
        r=(y*y-3)/6;
        s=1/(pp+pp-1);
        t=1/(qq+qq-1);
        h=2/(s+t);
        w=y*Math.sqrt(h+r)/h-(t-s)*(r+5/6-2/(3*h));
        xinbta=pp/(pp+qq*Math.exp(w+w));
        break;
        }
r=qq+qq;
t=1/(9*qq);
t=r*Math.pow(1-t+y*Math.sqrt(t),3);
if (t <= 0)
        {
        xinbta=1-Math.exp((Math.log((1-a)*qq)+beta)/qq);
        break;
        }

t=(4*pp+r-2)/t;
if (t<=1)
        {
        xinbta=Math.exp((Math.log(a*pp)+beta)/pp);
        break;
        }

xinbta=1-2/(t+1);
break;
}


r=1-pp;
t=1-qq;
yprev=0;
sq=1;
prev=1;
if (xinbta <lower) xinbta=lower;
if (xinbta>upper)  xinbta=upper;

bigfor: 
for(;;)
        {
        y=dfbeta(xinbta, pp, qq);
        y=(y-a)*Math.exp(beta +r*Math.log(xinbta)+t*Math.log(1-xinbta));
        if (y*yprev <0) prev=sq;
        g=1;

        for(;;){
        adj=g*y;
        sq=adj*adj;
        if (sq >=prev)
                {
                g=g/3;
                continue;
                }
        tx=xinbta-adj;
        if ( (tx <0) || (tx >1) )
                {
                g=g/3;
                continue;
                }
        if (prev <= acu) break bigfor;
        if (y*y <= acu) break bigfor;
        if ( (tx ==0) || (tx ==1) )
                {
                g=g/3;
                continue;
                }
        if (tx==xinbta) break bigfor;
        xinbta=tx;
        yprev=y;
        break;
        }
        }

if (index==1) return(1-xinbta);
return(xinbta);
}


function bplower(trials, suc, alpha)
{
if (suc==0) return 0.0;
return ppbeta(alpha, suc, (trials-suc+1));
}


function bpupper(trials, suc, alpha)
{
if(suc==trials) return 1.0;
return  ppbeta(1-alpha, (suc+1), (trials-suc));
}


