gibney.org
:
Technology
:
Javascript
:
Experiments
:
Four Fields Test
(Entry Nr. 323, by user 39 |
edit
)
<html> <script language="javascript"> 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)); } function hyper0(ai,n1_i,n_1i,ni) { if(!(n1_i|n_1i|ni)) { if(!(ai % 10 == 0)) { if(ai==sa+1) { sprob *= ((sn1_-sa)/(ai))*((sn_1-sa)/(ai+sn-sn1_-sn_1)); sa = ai; return sprob; } if(ai==sa-1) { sprob *= ((sa)/(sn1_-ai))*((sa+sn-sn1_-sn_1)/(sn_1-ai)); sa = ai; return sprob; } } sa = ai; } else { sa = ai; sn1_=n1_i; sn_1=n_1i; sn=ni; } sprob = hyper_323(sa,sn1_,sn_1,sn); return sprob; } 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