Four Fields Test

gibney.org : Technology : Javascript : Experiments : Four Fields Test

(Entry Nr. 323, by user 39 | edit)

<html>

function Csp(x)
{
var r=Math.sqrt(x)
//----------------------------------------------------------------------//
//      A JAVASCRIPT PROGRAM FOR THE NORMAL DISTRIBUTION  INTEGRAL      //
//   coded by William Knight - University of New Brunswick - Canada     //
//----------------------------------------------------------------------//

//   The results appear in these global variables:
//       f   the probability density
//       F   the cumulative function
//       G   the complementary cumulative,  G = 1-F

var  exp = Math.exp(-r*r/2) ;
     f =  exp/2.506628274631 ;       // f = probability density

// COEFFICIENTS USED IN APPROXIMATIONS
     // Coefficient sets 1 and 2 are taken from page 91 or
     //       Statistical Computing
     //       by William J. Kennedy, Jr. and James E. Gentle
     // Coefficient set 3 was derived by the author

var  P10 = 2.4266795523053075E2;     var  Q10 = 2.1505887586986120E2;
     var  P11 = 2.1979261618294152E1;     var  Q11 = 9.1164905404514901E1;
     var  P12 = 6.9963834886191355  ;     var  Q12 = 1.5082797630407787E1;
     var  P13 =-3.5609843701815385E-2;    var  Q13 = 1;

var  P20 = 3.004592610201616005E2;   var  Q20 = 3.004592609569832933E2;
     var  P21 = 4.519189537118729422E2;   var  Q21 = 7.909509253278980272E2;
     var  P22 = 3.393208167343436870E2;   var  Q22 = 9.313540948506096211E2;
     var  P23 = 1.529892850469404039E2;   var  Q23 = 6.389802644656311665E2;
     var  P24 = 4.316222722205673530E1;   var  Q24 = 2.775854447439876434E2;
     var  P25 = 7.211758250883093659E0;   var  Q25 = 7.700015293522947295E1;
     var  P26 = 5.641955174789739711E-1;  var  Q26 = 1.278272731962942351E1;
     var  P27 =-1.368648573827167067E-7;  var  Q27 = 1;

var  P30 = 0.00026731935609105996;   var  P31 =  5.96764655143274;
     var  P32 = 1.195388047407337;        var  P33 = -63.013348444126016;
     var  P34 = 194.3681790673782;        var  P35 = -212.55469380263868;
     // APPROXIMATIONS
     if (r < 0.58)
     {
        var  xn = r/1.4142135623730952 ;
        var  s = xn*xn ;
        var  P = P10+s*(P11+s*(P12+s*P13)) ;  // Numerator
        var  Q = Q10+s*(Q11+s*(Q12+s*Q13)) ;  // Denominator
        F = (1+xn*(P/Q))/2
        G = 1-F ;
     }
     else

if (r < 5.63) // (0.58 < r < 5.63)
     {
        var  xn = r/1.4142135623730952 ;
        var  P = P20+xn*(P21+xn*(P22+xn*(P23+xn*(P24+xn*(P25+xn*(P26+xn*P27)))))) ;
        var  Q = Q20+xn*(Q21+xn*(Q22+xn*(Q23+xn*(Q24+xn*(Q25+xn*(Q26+xn*Q27)))))) ;
        G = (P/Q)*exp/2;
        F = 1-G ;
     }
     else         // (5.63 < r)
     {
        var  xn = 1/r ;
        var  P = P30+xn*(P31+xn*(P32+xn*(P33+xn*(P34+xn*P35)))) ;
        var  R = x+1/(x+2/(x+3/(x+4/(x+5/(x+P))))) ;
        G = f/R ;
        F = 1-G ;

};

// SWAP SIDES FOR NEGATIVE x
     if (x<0) { H=F ; F=G ; G=H ;  }
return 2*G
 };

function BinomialP( p, n, k ) {
   if (n >= 1000)  return BinomialPF( p, n, k );
   else {
      // term-by-term
      if ((k > n) || (p >= 1))  return 1;
      else {
         var  q = 1 - p;
         var  n1p = (n+1) * p;

var  t = n * Math.log(q);  // k = 0
         var  r = Math.exp(t);
         var  j = 1;
         while (j <= k) {
            t += Math.log( 1 + (n1p - j) / (j * q) );
            r += Math.exp(t);
            j++;
            }

return  r;
         }
      }
   }
function BinomialPF( p, n, k ) {
   // by Normal approximation }
   // Peizer & Pratt 1968, JASA 63: 1416-1456
   var  inv2 = 1/2, inv3 = 1/3, inv6 = 1/6;

if (k >= n)  z = 1;
   else {
      var q = 1 - p;

var s = k + inv2;
      var t = n - k - inv2;

var d1 = s + inv6 - (n + inv3) * p;
      var d2 = q /(s+inv2)  -  p/(t+inv2)  +  (q-inv2)/(n+1);
      d2 = d1 + 0.02 * d2;

var num = 1  +  q * g(s/(n*p))  +  p * g(t/(n*q));
      var den = (n + inv6) * p * q;
      var z = num / den;
      z = d2 * Math.sqrt(z);
      z = NormalP( z );
      }

return z;
   }
function g( x ) {
   // Peizer & Pratt 1968, JASA 63: 1416-1456
   var  switchlev = 0.1, z;

if (x == 0)  z = 1;
   else
      if (x == 1)  z = 0;
      else {

var d = 1 - x;

if (Math.abs(d) > switchlev)
            z = (1 - (x * x) + (2 * x * Math.log(x))) / (d * d);
         else {

z = d / 3;  // first term
            var di = d;  // d**1

for (var i = 2; i <= 7; i++) {
               di *= d;  // d**i
               z += (2 * di) / ((i+1) * (i+2));
               }
            }
         }
   return z;
   }

function NormalP( x ) {
   // Abramowitz & Stegun 26.2.19
   var
      d1 = 0.0498673470,
      d2 = 0.0211410061,
      d3 = 0.0032776263,
      d4 = 0.0000380036,
      d5 = 0.0000488906,
      d6 = 0.0000053830,
      a = Math.abs(x),
      t;

t = 1.0 + a*(d1+a*(d2+a*(d3+a*(d4+a*(d5+a*d6)))));

// to 16th power
   t *= t;  t *= t;  t *= t;  t *= t;
   t = 1.0 / (t+t);  // the MINUS 16th

if (x >= 0)  t = 1-t;
   return t;
   }

function Fmt(x) {
var v
if(x>=0) { v=''+(x+0.000005) } else { v=''+(x-0.000005) }
return v.substring(0,v.indexOf('.')+6)
}

function lngamm(z)
// Reference: "Lanczos, C. 'A precision approximation
// of the gamma function', J. SIAM Numer. Anal., B, 1, 86-96, 1964."
// Translation of  Alan Miller's FORTRAN-implementation
// See httphiddenlib.stat.cmu.edu/apstat/245
{
  var x = 0;
  x += 0.1659470187408462e-06/(z+7);
  x += 0.9934937113930748e-05/(z+6);
  x -= 0.1385710331296526    /(z+5);
  x += 12.50734324009056     /(z+4);
  x -= 176.6150291498386     /(z+3);
  x += 771.3234287757674     /(z+2);
  x -= 1259.139216722289     /(z+1);
  x += 676.5203681218835     /(z);
  x += 0.9999999999995183;
  return(Math.log(x)-5.58106146679532777-z+(z-0.5)*Math.log(z+6.5));
}

function lnfact(n)
{
  if(n<=1) return(0);
  return(lngamm(n+1));
}

function lnbico(n,k)
{
  return(lnfact(n)-lnfact(k)-lnfact(n-k));
}

function hyper_323(a,n1_,n_1,n)
{
  return(Math.exp(lnbico(n1_,a)+lnbico(n-n1_,n_1-a)-lnbico(n,n_1)));
}

var sa,sn1_,sn_1,sn,sprob;
function hyper(a)
{
  return(hyper0(a,0,0,0));
}

var sleft,sright,sless,slarg;
function exact(a,n1_,n_1,n)
{
  var p,i,j,prob;
  var max=n1_;
  if(n_1<max) max=n_1;
  var min = n1_+n_1-n;
  if(min<0) min=0;
  if(min==max)
  {
    sless = 1;
    sright= 1;
    sleft = 1;
    slarg = 1;
    return 1;
  }
  prob=hyper0(a,n1_,n_1,n);
  sleft=0;
  p=hyper(min);
  for(i=min+1; p<0.99999999999*prob; i++)
  {
    sleft += p;
    p=hyper(i);
  }
  i--;
  if(p<1.00000000001*prob) sleft += p;
  else i--;
  sright=0;
  p=hyper(max);
  for(j=max-1; p<0.99999999999*prob; j--)
  {
    sright += p;
    p=hyper(j);
  }
  j++;
  if(p<1.00000000001*prob) sright += p;
  else j++;
  if(Math.abs(i-a)<Math.abs(j-a))
  {
    sless = sleft;
    slarg = 1 - sleft + prob;
  }
  else
  {
    sless = 1 - sright + prob;
    slarg = sright;
  }
  return prob;
}
var newline="\n";
if (navigator.appVersion.length < 1) newline = "\r\n";
if (navigator.appVersion.lastIndexOf('Win') != -1) newline="\r\n";

var left,right,twotail;
var aold=-1;
var bold=-1;
var cold=-1;
var dold=-1;

n=5000;

x=0;

for (j=0; j<500; j++){
 x1=x2=0;
 for (i=0; i<n; i++)
  if(Math.random()>=.5) x1++;
 for (i=0; i<n; i++)
  if(Math.random()>=.5) x2++;
 
 a=x1; b= n-x1;
 c=x2; d= n-x2;

var n1_ = a+b;
  var n_1 = a+c;
  var n_   = a +b +c +d;
  var prob=exact(a,n1_,n_1,n_);

left = sless;
 right = slarg;
 twotail = sleft+sright; if(twotail > 1) twotail=1;

document.write('...'+a+' '+' '+b+' '+c+' '+d+' ');
  document.write(twotail+'  ');

}

</script>

</html>

Create a new entry at this position