//The Gamma, Beta and Erlang distribution online.
//copyright: Daan G Uitenbroek PhD

function ClearOutp()   
{  
var writeln="\n";  
if (navigator.appVersion.lastIndexOf('Win') != -1)
writeln="\r\n";   
document.form1.outp.value ="cleared"+writeln;  
}  
 
function array(len)  
{  
  this.length=len;  
   for (var i=1;i<=len;i++)  
      this[i]=0;  
  return this  
} 
  
function odd(number)  //determines if a number is odd  
{  
var isodd=0*1;  
if (Math.round(number/2+0.00001)==number/2)  isodd=1*1;  
return isodd  
}  
 
function intg(number)  //determines if a number is an integer 
{ 
var isint=0*1 
if (Math.round(number)==number) isint=1*1 
return isint 
} 

function factorial(start,startloop,fvalue) //factorial without recursion
{
factorialvalue=start
for (var index=startloop;index<=fvalue;index++)
 factorialvalue*=index

return factorialvalue 
}
 
function factorial2(start,andex,fvalue)  //returns the factorial value of a number 
{ 
factorialvalue=start 
if (andex<=fvalue)
factorial(factorialvalue*andex,andex*1+1*1,fvalue)  
 
return factorialvalue 
} // recursive procedure (doesn't work well in Netscape)

function logfactorial(start,startloop,fvalue) //logfactorial without recursion
{
factorialvalue=start
for (var index=startloop;index<=fvalue;index++)
 factorialvalue=factorialvalue*1+Math.log(index)*1

return factorialvalue 
}
 
function logfactorial2(start,andex,fvalue)  //returns the log factorial value of a number 
{ 
factorialvalue=start 
if (andex<=fvalue)  
  
logfactorial(factorialvalue*1+Math.log(andex)*1,andex*1+1*1,fvalue)  
 
return factorialvalue 
} 
 
function truncate(real,places)  // truncates 'real' numbers n 'places' behind the decimal dot 
{ 
var multi=Math.pow(10,places); 
return (Math.round(real*multi)/multi) 
} 
 
function remain(num)  //returns the remainder of a real number 
{ 
var remainder=num-Math.round(num-0.5) 
 
return remainder 
} 

function erlangcumulative(b,c,x)
{
var sum=0
for (var index=0;index<c;index++)
  sum+=Math.exp(Math.log(x/b)*index-logfactorial(0,1,index))
sum*=Math.exp(-1*x/b)
return sum=1-sum
}

function gammadensity(b,c,x)
{
var dens=Math.pow((x/b),(c-1))
dens*=Math.exp(-1*x/b)
dens/=(b*gamma(c))
return dens
}

function cumpoisson(p,x)
{
var p0=Math.exp(-x)
var tot=1-p0;
for (var index=1;index<p;index++) 
{
     p0*=x/index; 
     tot-=p0; 
}
return tot
}

function gammam(x,p,q,perc) 
{ 
writeln="\n";  
if (navigator.appVersion.lastIndexOf('Win') != -1)
writeln="\r\n";  
if (x>0)
document.form1.outp.value +=writeln+writeln+"Gamma (X) function value: "+gamma(x)

if ((p>0)&(x>0)) 
 document.form1.outp.value +=writeln+writeln+"Incomplete Gamma (c,X) value: "+cumpoisson(p,x)

if ((p>0)&(q>0))
{
document.form1.outp.value +=writeln+writeln+"Gamma Distribution mean "
     +(p*q)+writeln
     +"Gamma Distribution variance= "+(p*q*q)
document.form1.outp.value +=writeln+"Gamma Distribution mode "+(q*(p-1))
}

if ((p>0)&(q>0)&(x!=0)) 
{if (x<=0.0)  
    document.form1.outp.value+=writeln+"negative gamma-X coefficient is not allowed" 

document.form1.outp.value +=writeln+writeln+"Gamma Density: "+gammadensity(q,p,x)
document.form1.outp.value +=writeln+"Cummulative Probability (Erlang): "+erlangcumulative(q,p,x)
}

if ((p>0)&(q>0)&(perc!=0))
 {var helper= inverterlang(q,p,perc)
 document.form1.outp.value +=writeln+writeln+"X: "+helper+"; G(X): "+gamma(helper) 
}
}
 
function gamma(x)  //algorithm 221 combined with exact factorial 
{ 
writeln="\n";  
if (navigator.appVersion.lastIndexOf('Win') != -1)
writeln="\r\n";  

if ((x>3) & (intg(x)==1)) var gamma=factorial(1,1,x-1) 
else 
{ 
a=new array(11);  
a[0]=1.0 
a[1]=.4227843370 
a[2]=.4118402518 
a[3]=.0815782188 
a[4]=.0742379076 
a[5]=-.0002109075 
a[6]=.0109736958 
a[7]=-.0024667480 
a[8]=.0015397681 
a[9]=-.0003442342 
a[10]=.0000677106 
 
var t=(x<=1) ? x : (x<=2) ? x-1 : x-2 
var p=a[10] 
for (var index=9;index>-1;index-=1) var p=t*p+a[index]*1 
var gamma=(x<=1) ? p/(x*(x*1+1*1)) : (x<=2) ? p/x : p 
}
return(gamma)
} 

function betamean(p,q) 
{
return (p/(p*1+q*1))
}

function betavariance(p,q) 
{
return (p*q/((p*1+q*1)*(p*1+q*1)*(p*1+q*1+1*1)))
}

function betamode(p,q) 
{
return ((p-1)/(p*1+q*1-2))
}
 
function beta(x,p,q,perc) 
{ 
writeln="\n";  
if (navigator.appVersion.lastIndexOf('Win') != -1)
writeln="\r\n";  

if ((p>0)&(q>0))
{
var beta1=gamma(p)*gamma(q)/gamma(p*1+q*1)
document.form1.outp.value +=writeln+writeln+"Beta function value eq: "+beta1 
if ((p*1+q*1)>=50) 
{var
beta2=logfactorial(0,1,p-1)*1+logfactorial(0,1,q-1)*1-logfactorial(0,1,p*1+q*1-1) 
document.form1.outp.value +=writeln+"log-Beta eq: " 
                          +beta2+"-> beta: "+Math.exp(beta2)} 
document.form1.outp.value +=writeln+writeln+"Beta Distribution mean "
     +betamean(p,q)+writeln
     +"Beta Distribution variance= "+betavariance(p,q)
if ((p>1)&(q>1))
 document.form1.outp.value +=writeln+"Beta Distribution mode "+betamode(p,q)
}

if ((p>0)&(q>0)&(x!=0)) Ibeta(x,p,q)
if ((p>0)&(q>0)&(perc!=0)) IRbeta(perc,p,q)
} 
 
function Ibeta(x,p,q)
{ 
writeln="\n";  
if (navigator.appVersion.lastIndexOf('Win') != -1)
writeln="\r\n";  
 
if (x<=0.0)  
     document.form1.outp.value +=writeln+"negative beta-E coefficient is not allowed" 
else if (x>=1.0) 
     document.form1.outp.value +=writeln+"beta-E coefficient larger than one is not allowed" 

var beta1=(gamma(p)*gamma(q)/gamma(p*1+q*1))
var density=Math.pow(x,(p-1))*Math.pow((1-x),(q-1))/beta1

document.form1.outp.value +=writeln+writeln+"Density (if calculable): "+density
document.form1.outp.value +=writeln+"Cummulative Probability (I-Beta): "+bino2(x,p,(p*1+q*1-1)) 
} 
 
function IRbeta(prob,p,q) 
{ 
writeln="\n";  
if (navigator.appVersion.lastIndexOf('Win') != -1)
writeln="\r\n";  
 
if (prob<=0.0)  
     document.form1.outp.value +=writeln+"ERROR negative p value is not allowed" 
else if (prob>=1.0) 
     document.form1.outp.value +=writeln+"ERROR p-value larger than one is not allowed" 
 
document.form1.outp.value +=writeln+writeln+"Beta-coefficient (E-value): "+invertbino(prob,p,(p*1+q*1-1)) 
} 
 
function bino2(m,k,nk)  
{  
var p0=Math.exp(Math.log(1-m)*nk)    //(1/m)^nk  
var tot=1-p0;  
for (var index=0;index<(k-1);index++)   
     {p0*=((nk-index)/(index*1+1*1))*(m/(1-m));  
     tot-=p0   
     }  
return tot;  
} 
 
function invertbino(p,x,n)  
{ 
var minu=0.0*1  
var maxu=1.0*1  
 
var mean=x/n;    
for (var index=1;maxu-minu > 0.0000001;index+=1)  
  {  
   if (bino2(mean,x,n)>p)   
                 maxu=mean  
   else  
                  minu=mean  
   mean=(maxu*1 + minu*1)*0.5  
   }  
return mean 
}  

function inverterlang(b,c,p)  
{ 
var minu=0.0*1  
var maxu=10000.0*1  
 
var mean=(c*b*b);    
for (var index=1;maxu-minu > 0.0000001;index+=1)  
  {  
   if (erlangcumulative(b,c,mean)>p)   
                 maxu=mean  
   else  
                  minu=mean  
   mean=(maxu*1 + minu*1)*0.5  
   }  
return mean 
}  

function gammaone(x) 
{ 
writeln="\n";  
if (navigator.appVersion.lastIndexOf('Win') != -1)
writeln="\r\n";  
 
var gammaval=gamma(x) 
 
document.form1.outp.value+=writeln+"gamma= "+gammaval 
}