"Biostatistics Documentation": C:/Documents and Settings/Stephen/Desktop/Biostatistics Library and Program/Biostatistics/Distributions.cs Source File

00001 using System;
00002 
00003 namespace StephenAshley.Biostatistics
00004 {
00005    using NUMBER = Decimal;
00011    public static class Distributions
00012    {
00024       public static NUMBER ProbabilityF(NUMBER f, Int32 degreesOfFreedomBetweenGroups,
00025          Int32 degreesOfFreedomWithinGroups)
00026       {
00027          NUMBER f1, f2, x, ft, vp, xx, theta;
00028          NUMBER sth, cth, sts, cts, a, b, xi, gamma, cbrf;
00029          const Int32 maxn = 500;
00030 
00031          NUMBER F = f;
00032          if (F == 0.0m)
00033             return (1.0m);
00034 
00035          f1 = degreesOfFreedomBetweenGroups;
00036          f2 = degreesOfFreedomWithinGroups;
00037          ft = 0.0m;
00038 
00039          x = f2 / (f2 + f1 * F);
00040          vp = f1 + f2 - 2.0m;
00041 
00042          if ((!(((1 & degreesOfFreedomBetweenGroups) == 1))) &&
00043             (degreesOfFreedomBetweenGroups <= maxn))
00044          {
00045             xx = 1.0m - x;
00046             f1 = f1 - 2.0m;
00047             while (f1 >= 1.0m)
00048             {
00049                vp = vp - 2.0m;
00050                ft = xx * vp / f1 * (1.0m + ft);
00051                f1 = f1 - 2.0m;
00052             }
00053 
00054             ft = DecMath.Exp(0.5m * f2 * DecMath.Log(x)) * (1.0m + ft);
00055             // ft = exp(0.5*f2*log(x))*(1.0+ft);
00056          }
00057 
00058          else if ((!(((1 & degreesOfFreedomWithinGroups) == 1))) && (degreesOfFreedomWithinGroups <= maxn))
00059          {
00060             f2 = f2 - 2.0m;
00061             while (f2 >= 1.0m)
00062             {
00063                vp = vp - 2.0m;
00064                ft = x * vp / f2 * (1.0m + ft);
00065                f2 = f2 - 2.0m;
00066             }
00067 
00068             ft = 1.0m - DecMath.Exp(0.5m * f1 * DecMath.Log(1.0m - x)) *
00069                (1.0m + ft);
00070          }
00071 
00072          else if (degreesOfFreedomBetweenGroups + degreesOfFreedomWithinGroups <= maxn)
00073          {
00074             theta = DecMath.Atan(DecMath.Sqrt(f1 * F / f2));
00075             sth = DecMath.Sin(theta);
00076             cth = DecMath.Cos(theta);
00077             sts = sth * sth;
00078             cts = cth * cth;
00079             a = 0.0m;
00080             b = 0.0m;
00081             if (degreesOfFreedomWithinGroups > 1)
00082             {
00083                f2 = f2 - 2.0m;
00084                while (f2 >= 2.0m)
00085                {
00086                   a = cts * (f2 - 1.0m) / f2 * (1.0m + a);
00087                   f2 = f2 - 2.0m;
00088                }
00089 
00090                a = sth * cth * (1.0m + a);
00091             }
00092 
00093             a = theta + a;
00094             if (degreesOfFreedomBetweenGroups > 1)
00095             {
00096                f1 = f1 - 2.0m;
00097                while (f1 >= 2.0m)
00098                {
00099                   vp = vp - 2.0m;
00100                   b = sts * vp / f1 * (1.0m + b);
00101                   f1 = f1 - 2.0m;
00102                }
00103 
00104                gamma = 1.0m;
00105                f2 = 0.5m * degreesOfFreedomWithinGroups;
00106                xi = 1.0m;
00107                while (xi <= f2)
00108                {
00109                   gamma = xi * gamma / (xi - 0.5m);
00110                   xi = xi + 1.0m;
00111                }
00112 
00113                b = gamma * sth * DecMath.Exp(degreesOfFreedomWithinGroups * DecMath.Log(cth)) * (1.0m + b);
00114             }
00115 
00116             ft = 1.0m + (2m / DecMath.PI) * (b - a);     /* 0.6366...=2/pi */
00117          }
00118          else
00119          {
00120             f1 = 2.0m / (9.0m * f1);
00121             f2 = 2.0m / (9.0m * f2);
00122             cbrf = DecMath.Exp(0.33333333333333333333333333333m * DecMath.Log(F));
00123             ft = ALNorm(-((1.0m - f2) * cbrf + f1 - 1.0m) / DecMath.Sqrt(f2 * (cbrf * cbrf) + f1), false);
00124          }
00125 
00126          if (ft < 0.001m)
00127             return (0.000m);
00128          else
00129             return (ft);
00130       }
00131 
00132       private static Decimal ALNorm(Decimal z, Boolean upper)
00133       {
00134          /* Reference:
00135             I. D. Hill 
00136             Algorithm AS 66: "The Normal Integral"
00137             Applied Statistics
00138          */
00139          double X = 3.8052e-8;
00140          Decimal LTONE = 7.0m,
00141                 UTZERO = 18.66m,
00142                 CON = 1.28m,
00143                 a1 = 0.398942280444m,
00144                 a2 = 0.399903438504m,
00145                 a3 = 5.75885480458m,
00146                 a4 = 29.8213557808m,
00147                 a5 = 2.62433121679m,
00148                 a6 = 48.6959930692m,
00149                 a7 = 5.92885724438m,
00150                 b1 = 0.398942280385m,
00151                 b2 = (Decimal)X,
00152                 b3 = 1.00000615302m,
00153                 b4 = 3.98064794e-4m,
00154                 b5 = 1.986153813664m,
00155                 b6 = 0.151679116635m,
00156                 b7 = 5.29330324926m,
00157                 b8 = 4.8385912808m,
00158                 b9 = 15.1508972451m,
00159                 b10 = 0.742380924027m,
00160                 b11 = 30.789933034m,
00161                 b12 = 3.99019417011m;
00162          Decimal y, alnorm;
00163 
00164          if (z < 0)
00165          {
00166             upper = !upper;
00167             z = -z;
00168          }
00169          if (z <= LTONE || upper && z <= UTZERO)
00170          {
00171             y = 0.5m * z * z;
00172             if (z > CON)
00173             {
00174                alnorm = b1 * DecMath.Exp(-y) / (z - b2 + b3 / (z + b4 + b5 /
00175                   (z - b6 + b7 / (z + b8 - b9 / (z + b10 + b11 / (z + b12))))));
00176             }
00177             else
00178             {
00179                alnorm = 0.5m - z * (a1 - a2 * y / (y + a3 - a4 / (y + a5 + a6 /
00180                   (y + a7))));
00181             }
00182          }
00183          else
00184          {
00185             alnorm = 0;
00186          }
00187          if (!upper) alnorm = 1 - alnorm;
00188          return (alnorm);
00189       }
00190 
00191 
00199       public static NUMBER ProbabilityT(NUMBER t, Int32 degreesOfFreedom)
00200       {
00201          NUMBER decF, decA, decB, decS, decC, decFk, decPst;
00202          Int32 iM2, iOe, iKs, k;
00203 
00204          decF = degreesOfFreedom;
00205          decA = DecMath.Abs(t) / DecMath.Sqrt(decF);  /* absolute value added to make two tailed */
00206          decB = decF / (decF + (t * t));
00207          iM2 = degreesOfFreedom - 2;
00208          iOe = degreesOfFreedom - 2 * (Int32)(degreesOfFreedom / 2);
00209          decS = 1.0m;
00210          decC = 1.0m;
00211          iKs = 2 + iOe;
00212          decFk = iKs;
00213          if (iM2 >= 2)
00214          {
00215             k = iKs;
00216             while ((k == iKs) || (k <= iM2))
00217             {
00218                decC = decC * decB * (decFk - 1.0m) / decFk;
00219                decS = decS + decC;
00220                decFk = decFk + 2.0m;
00221                k = k + 2;
00222             }
00223          }
00224 
00225          if (iOe > 0)
00226          {
00227             if (degreesOfFreedom == 1)
00228                decS = 0;
00229             decPst = 2.0m * (1.0m - (0.5m +
00230                (decA * decB * decS + DecMath.Atan(decA)) * 0.3183098862m)); /* 2 * to get two tail */
00231          }
00232          else
00233             decPst = 2.0m * (1.0m - (0.5m + 0.5m * decA *
00234                DecMath.Sqrt(decB) * decS));   /* 2 * to get two tail */
00235 
00236          if (decPst < 0.001m)
00237             return (0.0m);
00238          else
00239             return decPst;
00240       }
00241       /* PROBST */
00242 
00243       /* function */
00244       /* Uses Algorithm 346, F-test probabilities, by J. Morris, CACM 12:184-185,
00245       1969 to find upper tail probability associated with a given value of F. */
00246 
00254       public static double ProbabilityChiSq(double chiSquared, Int32 degreesOfFreedom)
00255       {
00256          double dP = 0.0;
00257          double dQ = 0.0;
00258          double dDF = (double)degreesOfFreedom;
00259          int iStatus = 0;
00260          double dBounds = 0.0;
00261 
00262          cdfchi(ref dP, ref dQ, ref chiSquared, ref dDF, ref iStatus, ref dBounds);
00263 
00264          return dQ;
00265       }
00266       //This subroutrine was revised extensively 3-20-01
00274       public static NUMBER ProbabilityRho(NUMBER rho, Int16 points)
00275       {
00276          NUMBER b, xx, y, u;
00277          NUMBER prho, rtemp, IS;
00278          if (points < 7) return (-99.0m);
00279          if (DecMath.Abs(rho) < 0.001m) return (1.0m);
00280          if (DecMath.Abs(rho) > 0.999m) return (0.001m);
00281          NUMBER NPOINTS = points;
00282          IS = NPOINTS * ((NPOINTS * NPOINTS) - 1.0m) * (1.0m - DecMath.Abs(rho)) / 6.0m;
00283 
00284          /* Evaluation by Edgeworth Series */
00285          b = 1.0m / NPOINTS;
00286          xx = (6.0m * (IS - 1.0m) * b / (1.0m / (b * b) - 1.0m) - 1.0m) *
00287             DecMath.Sqrt(1.0m / b - 1.0m);
00288          y = xx * xx;
00289          u = xx * b * (0.2274m + b * (0.2531m + 0.1745m * b) + y * (-0.0758m + b
00290             * (0.1033m + 0.3932m * b) - y * b * (0.0879m + 0.0151m * b - y
00291             * (0.0072m - 0.0831m * b + y * b * (0.0131m - 0.00046m * y)))));
00292 
00293          rtemp = u / DecMath.Exp(y / 2.0m) + ALNorm(xx, true);
00294          prho = 2.0m * (1.0m - rtemp);
00295          if (prho < 0.001m) prho = 0.001m;
00296          if (prho > 1.0m) prho = 1.0m;
00297          return (prho);
00298       }
00299 
00305       public static NUMBER ProbabilityZ(NUMBER z)
00306       {
00307          NUMBER pz;
00308 
00309          pz = 2.0m * ALNorm(DecMath.Abs(z), true);
00310          // TODO Why?
00311          if (pz < 0.001m)
00312             return (0.0m);
00313          else
00314             return (pz);
00315       }
00316 
00317       private static void cdfchi(ref double p, ref double q, ref double x, ref double df,
00318        ref int status, ref double bound)
00319       {
00320          cumchi(ref x, ref df, ref p, ref q);
00321          return;
00322       } // end of cdfchi
00323 
00324       private static void cumchi(ref double x, ref double df, ref double cum, ref double ccum)
00325       {
00326          double a, xx;
00327          a = df * 0.5e0;
00328          xx = x * 0.5e0;
00329          cumgam(ref xx, ref a, ref cum, ref ccum);
00330          return;
00331       } // end of cumchi
00332 
00333       private static void cumgam(ref double x, ref double a, ref double cum, ref double ccum)
00334       {
00335          int K1 = 0;
00336          if (!(x <= 0.0e0)) goto S10;
00337          cum = 0.0e0;
00338          ccum = 1.0e0;
00339          return;
00340       S10:
00341          gratio(ref a, ref x, ref cum, ref ccum, ref K1);
00342          return;
00343       }
00344 
00345       private static void gratio(ref double a, ref double x, ref double ans, ref double qans,
00346          ref int ind)
00347       {
00348          double alog10 = 2.30258509299405e0;
00349          double d10 = -.185185185185185e-02;
00350          double d20 = .413359788359788e-02;
00351          double d30 = .649434156378601e-03;
00352          double d40 = -.861888290916712e-03;
00353          double d50 = -.336798553366358e-03;
00354          double d60 = .531307936463992e-03;
00355          double d70 = .344367606892378e-03;
00356          double rt2pin = .398942280401433e0;
00357          double rtpi = 1.77245385090552e0;
00358          double third = .333333333333333e0;
00359          double[] acc0 = {
00360                 5.0e-15,5.0e-7,5.0e-4
00361             };
00362          double[] big = {
00363                 20.0e0,14.0e0,10.0e0
00364             };
00365          double[] d0 = {
00366                 .833333333333333e-01,-.148148148148148e-01,.115740740740741e-02,
00367                 .352733686067019e-03,-.178755144032922e-03,.391926317852244e-04,
00368                 -.218544851067999e-05,-.185406221071516e-05,.829671134095309e-06,
00369                 -.176659527368261e-06,.670785354340150e-08,.102618097842403e-07,
00370                 -.438203601845335e-08
00371             };
00372          double[] d1 = {
00373                 -.347222222222222e-02,.264550264550265e-02,-.990226337448560e-03,
00374                 .205761316872428e-03,-.401877572016461e-06,-.180985503344900e-04,
00375                 .764916091608111e-05,-.161209008945634e-05,.464712780280743e-08,
00376                 .137863344691572e-06,-.575254560351770e-07,.119516285997781e-07
00377             };
00378          double[] d2 = {
00379                 -.268132716049383e-02,.771604938271605e-03,.200938786008230e-05,
00380                 -.107366532263652e-03,.529234488291201e-04,-.127606351886187e-04,
00381                 .342357873409614e-07,.137219573090629e-05,-.629899213838006e-06,
00382                 .142806142060642e-06
00383             };
00384          double[] d3 = {
00385                 .229472093621399e-03,-.469189494395256e-03,.267720632062839e-03,
00386                 -.756180167188398e-04,-.239650511386730e-06,.110826541153473e-04,
00387                 -.567495282699160e-05,.142309007324359e-05
00388             };
00389          double[] d4 = {
00390                 .784039221720067e-03,-.299072480303190e-03,-.146384525788434e-05,
00391                 .664149821546512e-04,-.396836504717943e-04,.113757269706784e-04
00392             };
00393          double[] d5 = {
00394                 -.697281375836586e-04,.277275324495939e-03,-.199325705161888e-03,
00395                 .679778047793721e-04
00396             };
00397          double[] d6 = {
00398                 -.592166437353694e-03,.270878209671804e-03
00399             };
00400          double[] e00 = {
00401                 .25e-3,.25e-1,.14e0
00402             };
00403          double[] x00 = {
00404                 31.0e0,17.0e0,9.7e0
00405             };
00406          int K1 = 1;
00407          int K2 = 0;
00408          double a2n, a2nm1, acc, am0, amn, an, an0, apn, b2n, b2nm1, c, c0, c1, c2, c3, c4,
00409             c5, c6, cma, e, e0, g, h, j, l, r, rta, rtx, s, sum, t, t1, tol, twoa, u, w, x0,
00410             y, z;
00411          int i, iop, m, max, n;
00412          double[] wk = new double[20];
00413          double T3;
00414          int T4, T5;
00415          double T6, T7;
00416 
00417          e = spmpar(ref K1);
00418          if (a < 0.0e0 || x < 0.0e0) goto S430;
00419          if (a == 0.0e0 && x == 0.0e0) goto S430;
00420          if (a * x == 0.0e0) goto S420;
00421          iop = ind + 1;
00422          if (iop != 1 && iop != 2) iop = 3;
00423 
00424          acc = Math.Max(acc0[iop - 1], e);
00425          acc = Math.Max(acc0[iop - 1], e);
00426          e0 = e00[iop - 1];
00427          x0 = x00[iop - 1];
00428          if (a >= 1.0e0) goto S10;
00429          if (a == 0.5e0) goto S390;
00430          if (x < 1.1e0) goto S160;
00431          t1 = a * Math.Log(x) - x;
00432          u = a * Math.Exp(t1);
00433          if (u == 0.0e0) goto S380;
00434          r = u * (1.0e0 + gam1(ref a));
00435          goto S250;
00436       S10:
00437          if (a >= big[iop - 1]) goto S30;
00438          if (a > x || x >= x0) goto S20;
00439          twoa = a + a;
00440          m = Convert.ToInt32(twoa);
00441          if (twoa != (double)m) goto S20;
00442          i = m / 2;
00443          if (a == (double)i) goto S210;
00444          goto S220;
00445       S20:
00446          t1 = a * Math.Log(x) - x;
00447          r = Math.Exp(t1) / Xgamm(ref a);
00448          goto S40;
00449       S30:
00450          l = x / a;
00451          if (l == 0.0e0) goto S370;
00452          s = 0.5e0 + (0.5e0 - l);
00453          z = rlog(ref l);
00454          if (z >= 700.0e0 / a) goto S410;
00455          y = a * z;
00456          rta = Math.Sqrt(a);
00457          if (Math.Abs(s) <= e0 / rta) goto S330;
00458          if (Math.Abs(s) <= 0.4e0) goto S270;
00459          t = Math.Pow(1.0e0 / a, 2.0);
00460          t1 = (((0.75e0 * t - 1.0e0) * t + 3.5e0) * t - 105.0e0) / (a * 1260.0e0);
00461          t1 -= y;
00462          r = rt2pin * rta * Math.Exp(t1);
00463       S40:
00464          if (r == 0.0e0) goto S420;
00465          if (x <= Math.Max(a, alog10)) goto S50;
00466          if (x < x0) goto S250;
00467          goto S100;
00468       S50:
00469          apn = a + 1.0e0;
00470          t = x / apn;
00471          wk[0] = t;
00472          for (n = 2; n <= 20; n++)
00473          {
00474             apn += 1.0e0;
00475             t *= (x / apn);
00476             if (t <= 1.0e-3) goto S70;
00477             wk[n - 1] = t;
00478          }
00479          n = 20;
00480       S70:
00481          sum = t;
00482          tol = 0.5e0 * acc;
00483       S80:
00484          apn += 1.0e0;
00485          t *= (x / apn);
00486          sum += t;
00487          if (t > tol) goto S80;
00488          max = n - 1;
00489          for (m = 1; m <= max; m++)
00490          {
00491             n -= 1;
00492             sum += wk[n - 1];
00493          }
00494          ans = r / a * (1.0e0 + sum);
00495          qans = 0.5e0 + (0.5e0 - ans);
00496          return;
00497       S100:
00498          amn = a - 1.0e0;
00499          t = amn / x;
00500          wk[0] = t;
00501          for (n = 2; n <= 20; n++)
00502          {
00503             amn -= 1.0e0;
00504             t *= (amn / x);
00505             if (Math.Abs(t) <= 1.0e-3) goto S120;
00506             wk[n - 1] = t;
00507          }
00508          n = 20;
00509       S120:
00510          sum = t;
00511       S130:
00512          if (Math.Abs(t) <= acc) goto S140;
00513          amn -= 1.0e0;
00514          t *= (amn / x);
00515          sum += t;
00516          goto S130;
00517       S140:
00518          max = n - 1;
00519          for (m = 1; m <= max; m++)
00520          {
00521             n -= 1;
00522             sum += wk[n - 1];
00523          }
00524          qans = r / x * (1.0e0 + sum);
00525          ans = 0.5e0 + (0.5e0 - qans);
00526          return;
00527       S160:
00528          an = 3.0e0;
00529          c = x;
00530          sum = x / (a + 3.0e0);
00531          tol = 3.0e0 * acc / (a + 1.0e0);
00532       S170:
00533          an += 1.0e0;
00534          c = -(c * (x / an));
00535          t = c / (a + an);
00536          sum += t;
00537          if (Math.Abs(t) > tol) goto S170;
00538          j = a * x * ((sum / 6.0e0 - 0.5e0 / (a + 2.0e0)) * x + 1.0e0 / (a + 1.0e0));
00539          z = a * Math.Log(x);
00540          h = gam1(ref a);
00541          g = 1.0e0 + h;
00542          if (x < 0.25e0) goto S180;
00543          if (a < x / 2.59e0) goto S200;
00544          goto S190;
00545       S180:
00546          if (z > -.13394e0) goto S200;
00547       S190:
00548          w = Math.Exp(z);
00549          ans = w * g * (0.5e0 + (0.5e0 - j));
00550          qans = 0.5e0 + (0.5e0 - ans);
00551          return;
00552       S200:
00553          l = rexp(ref z);
00554          w = 0.5e0 + (0.5e0 + l);
00555          qans = (w * j - l) * g - h;
00556          if (qans < 0.0e0) goto S380;
00557          ans = 0.5e0 + (0.5e0 - qans);
00558          return;
00559       S210:
00560          sum = Math.Exp(-x);
00561          t = sum;
00562          n = 1;
00563          c = 0.0e0;
00564          goto S230;
00565       S220:
00566          rtx = Math.Sqrt(x);
00567          sum = erfc1(ref K2, ref rtx);
00568          t = Math.Exp(-x) / (rtpi * rtx);
00569          n = 0;
00570          c = -0.5e0;
00571       S230:
00572          if (n == i) goto S240;
00573          n += 1;
00574          c += 1.0e0;
00575          t = x * t / c;
00576          sum += t;
00577          goto S230;
00578       S240:
00579          qans = sum;
00580          ans = 0.5e0 + (0.5e0 - qans);
00581          return;
00582       S250:
00583          tol = Math.Max(5.0e0 * e, acc);
00584          a2nm1 = a2n = 1.0e0;
00585          b2nm1 = x;
00586          b2n = x + (1.0e0 - a);
00587          c = 1.0e0;
00588       S260:
00589          a2nm1 = x * a2n + c * a2nm1;
00590          b2nm1 = x * b2n + c * b2nm1;
00591          am0 = a2nm1 / b2nm1;
00592          c += 1.0e0;
00593          cma = c - a;
00594          a2n = a2nm1 + cma * a2n;
00595          b2n = b2nm1 + cma * b2n;
00596          an0 = a2n / b2n;
00597          if (Math.Abs(an0 - am0) >= tol * an0) goto S260;
00598          qans = r * an0;
00599          ans = 0.5e0 + (0.5e0 - qans);
00600          return;
00601       S270:
00602          if (Math.Abs(s) <= 2.0e0 * e && a * e * e > 3.28e-3) goto S430;
00603          c = Math.Exp(-y);
00604          T3 = Math.Sqrt(y);
00605          w = 0.5e0 * erfc1(ref K1, ref T3);
00606          u = 1.0e0 / a;
00607          z = Math.Sqrt(z + z);
00608          if (l < 1.0e0) z = -z;
00609          T4 = iop - 2;
00610          if (T4 < 0) goto S280;
00611          else if (T4 == 0) goto S290;
00612          else goto S300;
00613       S280:
00614          if (Math.Abs(s) <= 1.0e-3) goto S340;
00615          c0 = ((((((((((((d0[12] * z + d0[11]) * z + d0[10]) * z + d0[9]) * z + d0[8]) *
00616             z + d0[7]) * z + d0[6]) * z + d0[5]) * z + d0[4]) * z + d0[3]) * z + d0[2]) *
00617             z + d0[1]) * z + d0[0]) * z - third;
00618          c1 = (((((((((((d1[11] * z + d1[10]) * z + d1[9]) * z + d1[8]) * z + d1[7]) * z
00619             + d1[6]) * z + d1[5]) * z + d1[4]) * z + d1[3]) * z + d1[2]) * z + d1[1]) * z
00620             + d1[0]) * z + d10;
00621          c2 = (((((((((d2[9] * z + d2[8]) * z + d2[7]) * z + d2[6]) * z + d2[5]) * z + d2[4])
00622             * z + d2[3]) * z + d2[2]) * z + d2[1]) * z + d2[0]) * z + d20;
00623          c3 = (((((((d3[7] * z + d3[6]) * z + d3[5]) * z + d3[4]) * z + d3[3]) * z + d3[2])
00624             * z + d3[1]) * z + d3[0]) * z + d30;
00625          c4 = (((((d4[5] * z + d4[4]) * z + d4[3]) * z + d4[2]) * z + d4[1]) * z + d4[0])
00626             * z + d40;
00627          c5 = (((d5[3] * z + d5[2]) * z + d5[1]) * z + d5[0]) * z + d50;
00628          c6 = (d6[1] * z + d6[0]) * z + d60;
00629          t = ((((((d70 * u + c6) * u + c5) * u + c4) * u + c3) * u + c2) * u + c1) * u + c0;
00630          goto S310;
00631       S290:
00632          c0 = (((((d0[5] * z + d0[4]) * z + d0[3]) * z + d0[2]) * z + d0[1]) * z + d0[0])
00633             * z - third;
00634          c1 = (((d1[3] * z + d1[2]) * z + d1[1]) * z + d1[0]) * z + d10;
00635          c2 = d2[0] * z + d20;
00636          t = (c2 * u + c1) * u + c0;
00637          goto S310;
00638       S300:
00639          t = ((d0[2] * z + d0[1]) * z + d0[0]) * z - third;
00640       S310:
00641          if (l < 1.0e0) goto S320;
00642          qans = c * (w + rt2pin * t / rta);
00643          ans = 0.5e0 + (0.5e0 - qans);
00644          return;
00645       S320:
00646          ans = c * (w - rt2pin * t / rta);
00647          qans = 0.5e0 + (0.5e0 - ans);
00648          return;
00649       S330:
00650          if (a * e * e > 3.28e-3) goto S430;
00651          c = 0.5e0 + (0.5e0 - y);
00652          w = (0.5e0 - Math.Sqrt(y) * (0.5e0 + (0.5e0 - y / 3.0e0)) / rtpi) / c;
00653          u = 1.0e0 / a;
00654          z = Math.Sqrt(z + z);
00655          if (l < 1.0e0) z = -z;
00656          T5 = iop - 2;
00657          if (T5 < 0) goto S340;
00658          else if (T5 == 0) goto S350;
00659          else goto S360;
00660       S340:
00661          c0 = ((((((d0[6] * z + d0[5]) * z + d0[4]) * z + d0[3]) * z + d0[2]) * z + d0[1])
00662             * z + d0[0]) * z - third;
00663          c1 = (((((d1[5] * z + d1[4]) * z + d1[3]) * z + d1[2]) * z + d1[1]) * z + d1[0])
00664             * z + d10;
00665          c2 = ((((d2[4] * z + d2[3]) * z + d2[2]) * z + d2[1]) * z + d2[0]) * z + d20;
00666          c3 = (((d3[3] * z + d3[2]) * z + d3[1]) * z + d3[0]) * z + d30;
00667          c4 = (d4[1] * z + d4[0]) * z + d40;
00668          c5 = (d5[1] * z + d5[0]) * z + d50;
00669          c6 = d6[0] * z + d60;
00670          t = ((((((d70 * u + c6) * u + c5) * u + c4) * u + c3) * u + c2) * u + c1) * u + c0;
00671          goto S310;
00672       S350:
00673          c0 = (d0[1] * z + d0[0]) * z - third;
00674          c1 = d1[0] * z + d10;
00675          t = (d20 * u + c1) * u + c0;
00676          goto S310;
00677       S360:
00678          t = d0[0] * z - third;
00679          goto S310;
00680       S370:
00681          /*
00682                               SPECIAL CASES
00683          */
00684          ans = 0.0e0;
00685          qans = 1.0e0;
00686          return;
00687       S380:
00688          ans = 1.0e0;
00689          qans = 0.0e0;
00690          return;
00691       S390:
00692          if (x >= 0.25e0) goto S400;
00693          T6 = Math.Sqrt(x);
00694          ans = erf1(ref T6);
00695          qans = 0.5e0 + (0.5e0 - ans);
00696          return;
00697       S400:
00698          T7 = Math.Sqrt(x);
00699          qans = erfc1(ref K2, ref T7);
00700          ans = 0.5e0 + (0.5e0 - qans);
00701          return;
00702       S410:
00703          if (Math.Abs(s) <= 2.0e0 * e) goto S430;
00704       S420:
00705          if (x <= a) goto S370;
00706          goto S380;
00707       S430:
00708          /*
00709                               ERROR RETURN
00710          */
00711          ans = 2.0e0;
00712          return;
00713       } // end of gratio
00714 
00715       private static double spmpar(ref int i)
00716       {
00717          int K1 = 4;
00718          int K2 = 8;
00719          int K3 = 9;
00720          int K4 = 10;
00721          double spmpar, b, binv, bm1, one, w, z;
00722          int emax, emin, ibeta, m;
00723 
00724          if (i > 1) goto S10;
00725          b = ipmpar(ref K1);
00726          m = ipmpar(ref K2);
00727          spmpar = Math.Pow(b, (double)(1 - m));
00728          return spmpar;
00729       S10:
00730          if (i > 2) goto S20;
00731          b = ipmpar(ref K1);
00732          emin = ipmpar(ref K3);
00733          one = 1.0;
00734          binv = one / b;
00735          w = Math.Pow(b, (double)(emin + 2));
00736          spmpar = w * binv * binv * binv;
00737          return spmpar;
00738       S20:
00739          ibeta = ipmpar(ref K1);
00740          m = ipmpar(ref K2);
00741          emax = ipmpar(ref K4);
00742          b = ibeta;
00743          bm1 = ibeta - 1;
00744          one = 1.0;
00745          z = Math.Pow(b, (double)(m - 1));
00746          w = ((z - one) * b + bm1) / (b * z);
00747          z = Math.Pow(b, (double)(emax - 2));
00748          spmpar = w * z * b * b;
00749          return spmpar;
00750       } // end of spmpar
00751 
00752       private static int ipmpar(ref int i)
00753       {
00754          int[] imach = new int[11];
00755          int ipmpar;
00756 
00757          imach[1] = 2;
00758          imach[2] = 31;
00759          imach[3] = 2147483647;
00760          imach[4] = 2;
00761          imach[5] = 24;
00762          imach[6] = -125;
00763          imach[7] = 128;
00764          imach[8] = 53;
00765          imach[9] = -1021;
00766          imach[10] = 1024;
00767 
00768 
00769          ipmpar = imach[i];
00770          return ipmpar;
00771       }
00772 
00773       private static double erfc1(ref int ind, ref double x)
00774       {
00775          double c = .564189583547756e0;
00776          double[] a = {
00777                 .771058495001320e-04,-.133733772997339e-02,.323076579225834e-01,
00778                 .479137145607681e-01,.128379167095513e+00
00779             };
00780          double[] b = {
00781                 .301048631703895e-02,.538971687740286e-01,.375795757275549e+00
00782             };
00783          double[] p = {
00784                 -1.36864857382717e-07,5.64195517478974e-01,7.21175825088309e+00,
00785                 4.31622272220567e+01,1.52989285046940e+02,3.39320816734344e+02,
00786                 4.51918953711873e+02,3.00459261020162e+02
00787             };
00788          double[] q = {
00789                 1.00000000000000e+00,1.27827273196294e+01,7.70001529352295e+01,
00790                 2.77585444743988e+02,6.38980264465631e+02,9.31354094850610e+02,
00791                 7.90950925327898e+02,3.00459260956983e+02
00792             };
00793          double[] r = {
00794                 2.10144126479064e+00,2.62370141675169e+01,2.13688200555087e+01,
00795                 4.65807828718470e+00,2.82094791773523e-01
00796             };
00797          double[] s = {
00798                 9.41537750555460e+01,1.87114811799590e+02,9.90191814623914e+01,
00799                 1.80124575948747e+01
00800             };
00801          int K1 = 1;
00802          double erfc1, ax, bot, e, t, top, w;
00803 
00804          ax = Math.Abs(x);
00805          if (ax > 0.5e0) goto S10;
00806          t = x * x;
00807          top = (((a[0] * t + a[1]) * t + a[2]) * t + a[3]) * t + a[4] + 1.0e0;
00808          bot = ((b[0] * t + b[1]) * t + b[2]) * t + 1.0e0;
00809          erfc1 = 0.5e0 + (0.5e0 - x * (top / bot));
00810          if (ind != 0) erfc1 = Math.Exp(t) * erfc1;
00811          return erfc1;
00812       S10:
00813          if (ax > 4.0e0) goto S20;
00814          top = ((((((p[0] * ax + p[1]) * ax + p[2]) * ax + p[3]) * ax + p[4]) * ax + p[5]) * ax + p[6]) * ax + p[
00815            7];
00816          bot = ((((((q[0] * ax + q[1]) * ax + q[2]) * ax + q[3]) * ax + q[4]) * ax + q[5]) * ax + q[6]) * ax + q[
00817            7];
00818          erfc1 = top / bot;
00819          goto S40;
00820       S20:
00821          if (x <= -5.6e0) goto S60;
00822          if (ind != 0) goto S30;
00823          if (x > 100.0e0) goto S70;
00824          if (x * x > -exparg(ref K1)) goto S70;
00825       S30:
00826          t = Math.Pow(1.0e0 / x, 2.0);
00827          top = (((r[0] * t + r[1]) * t + r[2]) * t + r[3]) * t + r[4];
00828          bot = (((s[0] * t + s[1]) * t + s[2]) * t + s[3]) * t + 1.0e0;
00829          erfc1 = (c - t * top / bot) / ax;
00830       S40:
00831          if (ind == 0) goto S50;
00832          if (x < 0.0e0) erfc1 = 2.0e0 * Math.Exp(x * x) - erfc1;
00833          return erfc1;
00834       S50:
00835          w = x * x;
00836          t = w;
00837          e = w - t;
00838          erfc1 = (0.5e0 + (0.5e0 - e)) * Math.Exp(-t) * erfc1;
00839          if (x < 0.0e0) erfc1 = 2.0e0 - erfc1;
00840          return erfc1;
00841       S60:
00842          erfc1 = 2.0e0;
00843          if (ind != 0) erfc1 = 2.0e0 * Math.Exp(x * x);
00844          return erfc1;
00845       S70:
00846          erfc1 = 0.0e0;
00847          return erfc1;
00848       }
00849 
00850       private static double gam1(ref double a)
00851       {
00852          double s1 = .273076135303957e+00;
00853          double s2 = .559398236957378e-01;
00854          double[] p = {
00855                 .577215664901533e+00,-.409078193005776e+00,-.230975380857675e+00,
00856                 .597275330452234e-01,.766968181649490e-02,-.514889771323592e-02,
00857                 .589597428611429e-03
00858             };
00859          double[] q = {
00860                 .100000000000000e+01,.427569613095214e+00,.158451672430138e+00,
00861                 .261132021441447e-01,.423244297896961e-02
00862             };
00863          double[] r = {
00864                 -.422784335098468e+00,-.771330383816272e+00,-.244757765222226e+00,
00865                 .118378989872749e+00,.930357293360349e-03,-.118290993445146e-01,
00866                 .223047661158249e-02,.266505979058923e-03,-.132674909766242e-03
00867             };
00868          double gam1, bot, d, t, top, w, T1;
00869          /*
00870               ..
00871               .. Executable Statements ..
00872          */
00873          t = a;
00874          d = a - 0.5e0;
00875          if (d > 0.0e0) t = d - 0.5e0;
00876          T1 = t;
00877          if (T1 < 0) goto S40;
00878          else if (T1 == 0) goto S10;
00879          else goto S20;
00880       S10:
00881          gam1 = 0.0e0;
00882          return gam1;
00883       S20:
00884          top = (((((p[6] * t + p[5]) * t + p[4]) * t + p[3]) * t + p[2]) * t + p[1])
00885             * t + p[0];
00886          bot = (((q[4] * t + q[3]) * t + q[2]) * t + q[1]) * t + 1.0e0;
00887          w = top / bot;
00888          if (d > 0.0e0) goto S30;
00889          gam1 = a * w;
00890          return gam1;
00891       S30:
00892          gam1 = t / a * (w - 0.5e0 - 0.5e0);
00893          return gam1;
00894       S40:
00895          top = (((((((r[8] * t + r[7]) * t + r[6]) * t + r[5]) * t + r[4]) * t + r[3])
00896             * t + r[2]) * t + r[1]) * t + r[0];
00897          bot = (s2 * t + s1) * t + 1.0e0;
00898          w = top / bot;
00899          if (d > 0.0e0) goto S50;
00900          gam1 = a * (w + 0.5e0 + 0.5e0);
00901          return gam1;
00902       S50:
00903          gam1 = t * w / a;
00904          return gam1;
00905       }
00906 
00907       private static double Xgamm(ref double a)
00908       {
00909          double d = .41893853320467274178e0;
00910          double pi = 3.1415926535898e0;
00911          double r1 = .820756370353826e-03;
00912          double r2 = -.595156336428591e-03;
00913          double r3 = .793650663183693e-03;
00914          double r4 = -.277777777770481e-02;
00915          double r5 = .833333333333333e-01;
00916          double[] p = {
00917              .539637273585445e-03,.261939260042690e-02,.204493667594920e-01,
00918              .730981088720487e-01,.279648642639792e+00,.553413866010467e+00,1.0e0
00919          };
00920          double[] q = {
00921              -.832979206704073e-03,.470059485860584e-02,.225211131035340e-01,
00922              -.170458969313360e+00,-.567902761974940e-01,.113062953091122e+01,1.0e0
00923          };
00924          int K2 = 3;
00925          int K3 = 0;
00926          double Xgamm, bot, g, lnx, t, top, w, x, z;
00927          int i, j, m, n, T1;
00928 
00929          Xgamm = 0.0e0;
00930          x = a;
00931          if (Math.Abs(a) >= 15.0e0) goto S110;
00932          t = 1.0e0;
00933          m = Convert.ToInt32(a) - 1;
00934          T1 = m;
00935          if (T1 < 0) goto S40;
00936          else if (T1 == 0) goto S30;
00937          else goto S10;
00938       S10:
00939          for (j = 1; j <= m; j++)
00940          {
00941             x -= 1.0e0;
00942             t = x * t;
00943          }
00944       S30:
00945          x -= 1.0e0;
00946          goto S80;
00947       S40:
00948          t = a;
00949          if (a > 0.0e0) goto S70;
00950          m = -m - 1;
00951          if (m == 0) goto S60;
00952          for (j = 1; j <= m; j++)
00953          {
00954             x += 1.0e0;
00955             t = x * t;
00956          }
00957       S60:
00958          x += (0.5e0 + 0.5e0);
00959          t = x * t;
00960          if (t == 0.0e0) return Xgamm;
00961       S70:
00962          if (Math.Abs(t) >= 1.0e-30) goto S80;
00963          if (Math.Abs(t) * spmpar(ref K2) <= 1.0001e0) return Xgamm;
00964          Xgamm = 1.0e0 / t;
00965          return Xgamm;
00966       S80:
00967          top = p[0];
00968          bot = q[0];
00969          for (i = 1; i < 7; i++)
00970          {
00971             top = p[i] + x * top;
00972             bot = q[i] + x * bot;
00973          }
00974          Xgamm = top / bot;
00975          if (a < 1.0e0) goto S100;
00976          Xgamm *= t;
00977          return Xgamm;
00978       S100:
00979          Xgamm /= t;
00980          return Xgamm;
00981       S110:
00982          if (Math.Abs(a) >= 1.0e3) return Xgamm;
00983          double s = 0.0;
00984          if (a <= 0.0e0)
00985          {
00986             x = -a;
00987             n = Convert.ToInt32(x);
00988             t = x - (double)n;
00989             if (t > 0.9e0) t = 1.0e0 - t;
00990             s = 0.0;
00991             s = Math.Sin(pi * t) / pi;
00992 
00993             if ((n % 2) == 0)
00994                s = -s;
00995             if (s == 0.0e0)
00996                return Xgamm;
00997          }
00998 
00999 
01000          t = 1.0e0 / (x * x);
01001          g = ((((r1 * t + r2) * t + r3) * t + r4) * t + r5) / x;
01002          lnx = Math.Log(x);
01003          z = x;
01004          g = d + g + (z - 0.5e0) * (lnx - 1.0e0);
01005          w = g;
01006          t = g - w;
01007          if (w > 0.99999e0 * exparg(ref K3))
01008             return Xgamm;
01009          Xgamm = Math.Exp(w) * (1.0e0 + t);
01010 
01011          if (a < 0.0e0)
01012          {
01013             Xgamm = 1.0e0 / (Xgamm * s) / x;
01014          }
01015          return Xgamm;
01016       }
01017 
01018       private static double rlog(ref double x)
01019       {
01020          double a = .566749439387324e-01;
01021          double b = .456512608815524e-01;
01022          double p0 = .333333333333333e+00;
01023          double p1 = -.224696413112536e+00;
01024          double p2 = .620886815375787e-02;
01025          double q1 = -.127408923933623e+01;
01026          double q2 = .354508718369557e+00;
01027          double rlog, r, t, u, w, w1;
01028 
01029          if (x < 0.61e0 || x > 1.57e0) goto S40;
01030          if (x < 0.82e0) goto S10;
01031          if (x > 1.18e0) goto S20;
01032          u = x - 0.5e0 - 0.5e0;
01033          w1 = 0.0e0;
01034          goto S30;
01035       S10:
01036          u = x - 0.7e0;
01037          u /= 0.7e0;
01038          w1 = a - u * 0.3e0;
01039          goto S30;
01040       S20:
01041          u = 0.75e0 * x - 1.0e0;
01042          w1 = b + u / 3.0e0;
01043       S30:
01044          r = u / (u + 2.0e0);
01045          t = r * r;
01046          w = ((p2 * t + p1) * t + p0) / ((q2 * t + q1) * t + 1.0e0);
01047          rlog = 2.0e0 * t * (1.0e0 / (1.0e0 - r) - r * w) + w1;
01048          return rlog;
01049       S40:
01050          r = x - 0.5e0 - 0.5e0;
01051          rlog = r - Math.Log(x);
01052          return rlog;
01053       }
01054 
01055       private static double rexp(ref double x)
01056       {
01057          double p1 = .914041914819518e-09;
01058          double p2 = .238082361044469e-01;
01059          double q1 = -.499999999085958e+00;
01060          double q2 = .107141568980644e+00;
01061          double q3 = -.119041179760821e-01;
01062          double q4 = .595130811860248e-03;
01063          double rexp, w;
01064 
01065          if (Math.Abs(x) > 0.15e0) goto S10;
01066          rexp = x * (((p2 * x + p1) * x + 1.0e0) / ((((q4 * x + q3) * x + q2)
01067             * x + q1) * x + 1.0e0));
01068          return rexp;
01069       S10:
01070          w = Math.Exp(x);
01071          if (x > 0.0e0) goto S20;
01072          rexp = w - 0.5e0 - 0.5e0;
01073          return rexp;
01074       S20:
01075          rexp = w * (0.5e0 + (0.5e0 - 1.0e0 / w));
01076          return rexp;
01077       }
01078 
01079       private static double exparg(ref int l)
01080       {
01081          int K1 = 4;
01082          int K2 = 9;
01083          int K3 = 10;
01084          double exparg, lnb;
01085          int b, m;
01086          /*
01087               ..
01088               .. Executable Statements ..
01089          */
01090          b = ipmpar(ref K1);
01091          if (b != 2) goto S10;
01092          lnb = .69314718055995e0;
01093          goto S40;
01094       S10:
01095          if (b != 8) goto S20;
01096          lnb = 2.0794415416798e0;
01097          goto S40;
01098       S20:
01099          if (b != 16) goto S30;
01100          lnb = 2.7725887222398e0;
01101          goto S40;
01102       S30:
01103          lnb = Math.Log((double)b);
01104       S40:
01105          if (l == 0) goto S50;
01106          m = ipmpar(ref K2) - 1;
01107          exparg = 0.99999e0 * ((double)m * lnb);
01108          return exparg;
01109       S50:
01110          m = ipmpar(ref K3);
01111          exparg = 0.99999e0 * ((double)m * lnb);
01112          return exparg;
01113       }
01114 
01115       private static double erf1(ref double x)
01116       {
01117          double c = .564189583547756e0;
01118          double[] a = {
01119              .771058495001320e-04,-.133733772997339e-02,.323076579225834e-01,
01120              .479137145607681e-01,.128379167095513e+00
01121          };
01122          double[] b = {
01123              .301048631703895e-02,.538971687740286e-01,.375795757275549e+00
01124          };
01125          double[] p = {
01126              -1.36864857382717e-07,5.64195517478974e-01,7.21175825088309e+00,
01127              4.31622272220567e+01,1.52989285046940e+02,3.39320816734344e+02,
01128              4.51918953711873e+02,3.00459261020162e+02
01129          };
01130          double[] q = {
01131              1.00000000000000e+00,1.27827273196294e+01,7.70001529352295e+01,
01132              2.77585444743988e+02,6.38980264465631e+02,9.31354094850610e+02,
01133              7.90950925327898e+02,3.00459260956983e+02
01134          };
01135          double[] r = {
01136              2.10144126479064e+00,2.62370141675169e+01,2.13688200555087e+01,
01137              4.65807828718470e+00,2.82094791773523e-01
01138          };
01139          double[] s = {
01140              9.41537750555460e+01,1.87114811799590e+02,9.90191814623914e+01,
01141              1.80124575948747e+01
01142          };
01143          double erf1, ax, bot, t, top, x2;
01144          /*
01145               ..
01146               .. Executable Statements ..
01147          */
01148          ax = Math.Abs(x);
01149          if (ax > 0.5e0) goto S10;
01150          t = x * x;
01151          top = (((a[0] * t + a[1]) * t + a[2]) * t + a[3]) * t + a[4] + 1.0e0;
01152          bot = ((b[0] * t + b[1]) * t + b[2]) * t + 1.0e0;
01153          erf1 = x * (top / bot);
01154          return erf1;
01155       S10:
01156          if (ax > 4.0e0) goto S20;
01157          top = ((((((p[0] * ax + p[1]) * ax + p[2]) * ax + p[3]) * ax + p[4]) *
01158             ax + p[5]) * ax + p[6]) * ax + p[7];
01159          bot = ((((((q[0] * ax + q[1]) * ax + q[2]) * ax + q[3]) * ax + q[4]) *
01160             ax + q[5]) * ax + q[6]) * ax + q[7];
01161          erf1 = 0.5e0 + (0.5e0 - Math.Exp(-(x * x)) * top / bot);
01162          if (x < 0.0e0) erf1 = -erf1;
01163          return erf1;
01164       S20:
01165          if (ax >= 5.8e0) goto S30;
01166          x2 = x * x;
01167          t = 1.0e0 / x2;
01168          top = (((r[0] * t + r[1]) * t + r[2]) * t + r[3]) * t + r[4];
01169          bot = (((s[0] * t + s[1]) * t + s[2]) * t + s[3]) * t + 1.0e0;
01170          erf1 = (c - top / (x2 * bot)) / ax;
01171          erf1 = 0.5e0 + (0.5e0 - Math.Exp(-x2) * erf1);
01172          if (x < 0.0e0) erf1 = -erf1;
01173          return erf1;
01174       S30:
01175          erf1 = fifdsign(1.0e0, x);
01176          return erf1;
01177       }
01178 
01179       private static double fifdsign(double mag, double sign)
01180       {
01181          if (mag < 0) mag = -mag;
01182          if (sign < 0) mag = -mag;
01183          return mag;
01184       }
01185 
01186       /* function */
01187       /* Computes two-tail probabilities for Spearman’s decRho using Algorithm
01188       AS 69 from Applied Statistics 24:377-379, 1975 by Best and Roberts.  This
01189       only works for n>=7, where approximate solution is reasonably accurate. */
01190 
01191    }
01192 }