// **** The aurthor is grateful to John C. Pezzullo
// for publishing these routines and making them available

// original conent at http://statpages.org/scicalc.html


var Pi=Math.PI; var PiD2=Pi/2; var PiD4=Pi/4; var Pi2=2*Pi
var e=exp(1); var e10 = exp(0.1)
var Deg=180/Pi

function fact(n) {
    if(n==0 | n==1) { return 1 }
    if(n<0) { return fact(n+1)/(n+1) }
    if(n>1) { return n*fact(n-1) }
    if(n<0.5) { r = n } else { r = 1 - n }
    var r = 1 / (1 + r*( 0.577215664819072 + r*(-0.655878067489187 + r*(-0.042002698827786 + r*(0.166538990722800 + r*(-0.042197630554869 + r*(-0.009634403818022 + r*(0.007285315490429 + r*(-0.001331461501875 ) ) ) ) ) ) ) ) )
    if( n > 0.5 ) { r = n*(1-n)*Pi / (r*Math.sin(Pi*n)) }
    return r
    }
    
function gamma(n) { return fact(n-1) }
function abs(x) { return Math.abs(x) }
function sqrt(x) { return Math.sqrt(x) }
function exp(x) { return Math.exp(x) }
function ln(x) { return Math.log(x) }
function log10(x) { return Math.log(x)/Math.log(10) }
function log2(x) { return Math.log(x)/Math.log(2) }
function power(x,y) { return Math.pow(x,y) }
function sin(x) { return Math.sin(x) }
function cos(x) { return Math.cos(x) }
function tan(x) { return Math.sin(x)/Math.cos(x) }
function cot(x) { return Math.cos(x)/Math.sin(x) }
function sec(x) { return 1/Math.cos(x) }
function csc(x) { return 1/Math.sin(x) }
function atan(x) { return Math.atan(x) }
function asin(x) { return Math.asin(x) }
function acos(x) { return Math.acos(x) }
function acot(x) { return atan(1/x) }
function asec(x) { return acos(1/x) }
function acsc(x) { return asin(1/x) }

function sinh(x) { return (Math.exp(x)-Math.exp(-x))/2 }
function cosh(x) { return (Math.exp(x)+Math.exp(-x))/2 }
function tanh(x) { return sinh(x)/cosh(x) }
function coth(x) { return 1/tanh(x) }
function sech(x) { return 1/cosh(x) }
function csch(x) { return 1/sinh(x) }
function asinh(x) { return Math.log(x+Math.sqrt(x*x+1)) }
function acosh(x) { return Math.log(x+Math.sqrt(x*x-1)) }
function atanh(x) { return 0.5*Math.log((1+x)/(1-x)) }
function acoth(x) { return 0.5*Math.log((x+1)/(x-1)) }
function asech(x) { return Math.log(1/x+Math.sqrt(1/(x*x)+1)) }
function acsch(x) { return Math.log(1/x+Math.sqrt(1/(x*x)-1)) }

function chisq(x,n) {
    if(x>1000 | n>1000) { var q=norm((Math.pow(x/n,1/3)+2/(9*n)-1)/sqrt(2/(9*n)))/2; if (x>n) {return q}{return 1-q} }
    var p=Math.exp(-0.5*x); if((n%2)==1) { p=p*Math.sqrt(2*x/Pi) }
    var k=n; while(k>=2) { p=p*x/k; k=k-2 }
    var t=p; var a=n; while(t>1e-15*p) { a=a+2; t=t*x/a; p=p+t }
    return 1-p
    }
function norm(z) { var q=z*z
    if(abs(z)>7) {return (1-1/q+3/(q*q))*exp(-q/2)/(abs(z)*sqrt(PiD2))} else {return chisq(q,1) }
    }

function normd (x,m,s) {
return (1/(s*Math.sqrt(Pi2))*Math.exp(-(x-m)*(x-m)/(2*s*s)))
} 

function normstd(x) {
return (normd(x,0,1))
}

function normalpdf(x,m,s) {
if (s == undefined) {
		if (m == undefined) return (normd(x,0,1));
		else return(normd(x,m,1));
		} 
else return normd (x,m,s);
}


function gauss(z) { return ( (z<0) ? ( (z<-10) ? 0 : chisq(z*z,1)/2 ) : ( (z>10) ? 1 : 1-chisq(z*z,1)/2 ) ) }

function erf(z) { return ( (z<0) ? (2*gauss(sqrt(2)*z)-1) : (1-2*gauss(-sqrt(2)*z)) ) }

function studt(t,n) {
    t=Math.abs(t); var w=t/Math.sqrt(n); var th=Math.atan(w)
    if(n==1) { return 1-th/PiD2 }
    var sth=Math.sin(th); var cth=Math.cos(th)
    if((n%2)==1)
        { return 1-(th+sth*cth*StatCom(cth*cth,2,n-3,-1))/PiD2 }
        else
        { return 1-sth*StatCom(cth*cth,1,n-3,-1) }
    }
function fishf(f,n1,n2) {
    var x=n2/(n1*f+n2)
    if((n1%2)==0) { return StatCom(1-x,n2,n1+n2-4,n2-2)*Math.pow(x,n2/2) }
    if((n2%2)==0){ return 1-StatCom(x,n1,n1+n2-4,n1-2)*Math.pow(1-x,n1/2) }
    var th=Math.atan(Math.sqrt(n1*f/n2)); var a=th/PiD2; var sth=Math.sin(th); var cth=Math.cos(th)
    if(n2>1) { a=a+sth*cth*StatCom(cth*cth,2,n2-3,-1)/PiD2 }
    if(n1==1) { return 1-a }
    var c=4*StatCom(sth*sth,n2+1,n1+n2-4,n2-2)*sth*Math.pow(cth,n2)/Pi
    if(n2==1) { return 1-a+c/2 }
    var k=2; while(k<=(n2-1)/2) {c=c*k/(k-.5); k=k+1 }
    return 1-a+c
    }
function StatCom(q,i,j,b) {
    var zz=1; var z=zz; var k=i; while(k<=j) { zz=zz*q*k/(k-b); z=z+zz; k=k+2 }
    return z
    }
function anorm(p) { var v=0.5; var dv=0.5; var z=0
    while(dv>1e-15) { z=1/v-1; dv=dv/2; if(norm(z)>p) { v=v-dv } else { v=v+dv } }
    return z
    }

function agauss(p) { if(p>0.5) { return( sqrt(achisq(2*(1-p),1)) ) } else { return( -sqrt(achisq(2*p,1)) ) } }

function aerf(p) { return agauss(p/2+0.5)/sqrt(2) }

function achisq(p,n) { var v=0.5; var dv=0.5; var x=0
    while(dv>1e-15) { x=1/v-1; dv=dv/2; if(chisq(x,n)>p) { v=v-dv } else { v=v+dv } }
    return x
    }
function astudt(p,n) { var v=0.5; var dv=0.5; var t=0
    while(dv>1e-15) { t=1/v-1; dv=dv/2; if(studt(t,n)>p) { v=v-dv } else { v=v+dv } }
    return t
    }
function afishf(p,n1,n2) { var v=0.5; var dv=0.5; var f=0
    while(dv>1e-15) { f=1/v-1; dv=dv/2; if(fishf(f,n1,n2)>p) { v=v-dv } else { v=v+dv } }
    return f
    }