home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Programmierung / SOURCE.mdf / programm / msdos / c / yamp / dist.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1993-01-11  |  15.4 KB  |  582 lines
  1.  
  2. /*
  3. YAMP - Yet Another Matrix Program
  4. Version: 1.6
  5. Author: Mark Von Tress, Ph.D.
  6. Date: 01/11/93  
  7.  
  8. Copyright(c) Mark Von Tress 1993
  9. Portions of this code are (c) 1991 by Allen I. Holub and are used by
  10. permission
  11.  
  12. DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
  13. WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
  14. TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
  15. ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
  16. FROM USE OF THIS PROGRAM.
  17.  
  18. */
  19.  
  20.  
  21. #include "math.h"
  22. #include "stdio.h"
  23. #include "stdlib.h"
  24. #include "dist.h"
  25.  
  26. void nrerror(char error_text[])
  27. {
  28.     
  29.     printf("Numerical Recipes run-time error...\n");
  30.     printf("%s\n", error_text);
  31.     printf("...now exiting to system...\n");
  32.     exit(1);
  33. }
  34. double safelog( double x )
  35. {
  36.     if(!x) return -HUGE_VAL;
  37.     return log(x);
  38. }
  39.  
  40. double gammln(double xx)
  41. {
  42.     double x, tmp, ser;
  43.     static double cof[6] = { 76.18009173, - 86.50532033, 24.01409822,
  44.         - 1.231739516, 0.120858003e-2, - 0.536382e-5 };
  45.     int j;
  46.     
  47.     x = xx - 1.0;
  48.     tmp = x + 5.5;
  49.     tmp -= (x + 0.5) *safelog(tmp);
  50.     ser = 1.0;
  51.     for (j = 0; j <= 5; j++) {
  52.         x += 1.0;
  53.         ser += cof[j] / x;
  54.     }
  55.     return - tmp + safelog(2.50662827465*ser);
  56. }
  57.  
  58. #define ITMAX 100
  59. #define EPS 3.0e-9
  60.  
  61. double betacf(double a, double b, double x)
  62. {
  63.     double qap, qam, qab, em, tem, d;
  64.     double bz, bm = 1.0, bp, bpp;
  65.     double az = 1.0, am = 1.0, ap, app, aold;
  66.     int m;
  67.     
  68.     qab = a + b;
  69.     qap = a + 1.0;
  70.     qam = a - 1.0;
  71.     bz = 1.0- qab*x / qap;
  72.     for (m = 1; m <= ITMAX; m++) {
  73.         em = (double) m;
  74.         tem = em + em;
  75.         d = em*(b - em) *x / ((qam + tem) *(a + tem));
  76.         ap = az + d*am;
  77.         bp = bz + d*bm;
  78.         d = -(a + em) *(qab + em) *x / ((qap + tem) *(a + tem));
  79.         app = ap + d*az;
  80.         bpp = bp + d*bz;
  81.         aold = az;
  82.         am = ap / bpp;
  83.         bm = bp / bpp;
  84.         az = app / bpp;
  85.         bz = 1.0;
  86.         if (fabs(az - aold) < (EPS*fabs(az))) return az;
  87.     }
  88.     nrerror("BETACF: a or b too big, or ITMAX too small");
  89. }
  90.  
  91. #undef ITMAX
  92. #undef EPS
  93.  
  94. double betai(double a, double b, double x)
  95. {
  96.     double bt;
  97.     
  98.     if (x < 0.0 || x > 1.0) nrerror("BETAI: Bad x");
  99.     if (x == 0.0 || x == 1.0) return exp(x);
  100.     else
  101.         bt = (gammln(a + b) - gammln(a) - gammln(b)
  102.               + a*safelog(x) + b*safelog(1.0-x) );
  103.         
  104.     if (x < (a + 1.0) / (a + b + 2.0))
  105.         return (exp(bt + safelog(betacf(a, b, x)) - safelog(a)));
  106.     else
  107.         return (1.0-exp(bt + safelog(betacf(b, a, 1.0-x)) - safelog(b)));
  108. }
  109.  
  110.  
  111. double dmin(double x, double y)
  112. {
  113.     return ((x < y) ? x : y);
  114. }
  115.  
  116. double dmax(double x, double y)
  117. {
  118.     return ((x > y) ? x : y);
  119. }
  120.  
  121. double addlogc(double a, double b)
  122.     /* function to add two numbers a1 and b1 */
  123.         /* where a=log(a1) and b=log(b1) */
  124.         {
  125.             double del, sum;
  126.             del = dmax(a, b);
  127.             sum = safelog(exp(a - del) + exp(b - del)) + del;
  128.             return (sum);
  129. }
  130.  
  131. #define DMINEXP - 307
  132. #define DMAXEXP 308
  133. double ncbeta(double x, double a, double b, double lam)
  134. {
  135.     /*   noncentral beta distribution function
  136.     Vic Norton, Appl. Statist (1983), 32, #1, pp 84-85
  137.     this algorithm has been modified to work in logrithms
  138.     */
  139.     double s, r, p, n = 0, minn, elam;
  140.     double eps = 0.0000001, nmax = 500, leps;
  141.     
  142.     if ((a <= 0) || (b <= 0) || (x < 0) || (lam < 0))
  143.         nrerror("NCBETA: non-positive parameter");
  144.  
  145.     s = safelog(betai(a, b, x));
  146.     if (lam <= exp(DMINEXP)) return exp(s);     /* lam is too close to zero */
  147.     
  148.     r = safelog(lam);
  149.     if (lam > DMAXEXP) elam = exp(- DMAXEXP); else elam = exp(- lam);
  150.     leps = safelog(eps);
  151.     do {
  152.         n += 1;
  153.         p = safelog(betai((a + n), b, x));
  154.         s = addlogc(s, (r + p));
  155.         r += safelog(lam / (n + 1.0));
  156.         minn = (1 < (elam + lam / (n + 2))) ? 0 : safelog(lam / (n + 2)) - lam;
  157.     } while (((r + p + minn) > leps) && (n < nmax));
  158.     return (exp(- lam + s));
  159. }     /* ncbeta */
  160. #undef DMINEXP
  161. #undef DMAXEXP
  162.  
  163. double ncf(double f, double n1, double n2, double lam)
  164. {
  165.     double x, a, b, nc, lam2;
  166.     
  167.     if ((n1 <= 0) || (n2 <= 0) || (f < 0) || (lam < 0.0))
  168.         nrerror("NCF: non-positive parameter");
  169.     if (f == 0.0) return 0.0;
  170.     a = 0.5*n1;
  171.     b = 0.5*n2;
  172.     x = a*f / (a*f + b);
  173.     lam2 = 0.5*lam;
  174.     nc = ncbeta(x, a, b, lam2);
  175.     return (nc);
  176. }     /* non central f distribution using continuous degrees of freedom */
  177.  
  178.  
  179. double probf(double f, double n1, double n2)
  180. {
  181.     double a, b, x, lam2, nc;
  182.     if ((n1 <= 0.0) || (n2 <= 0.0) || (f < 0.0))
  183.         nrerror("PROBF: non-positive parameter");
  184.     if (f == 0.0) return (0.0);
  185.     a = 0.5*n1;
  186.     b = 0.5*n2;
  187.     x = a*f / (a*f + b);
  188.     
  189.     nc = ncbeta(x, a, b, 0.0);
  190.     return (nc);
  191. }     /* central f distribution using continuous degrees of freedom */
  192.  
  193.  
  194. /****************************** begin inverse ncentral f */
  195. /*          uses zbrak, and rtbis     */
  196.  
  197.  
  198. /***********************  normal distribution functions */
  199.  
  200. double probnormi(double u)     /* abramowitz and stegun 26.2.23 */
  201.     /* u in [0,1] -> returns normal quantiles */
  202.     {
  203.         double c0 = 2.515517, c1 = 0.802853, c2 = 0.010328;
  204.         double d1 = 1.432788, d2 = 0.189269, d3 = 0.001308;
  205.         double t, xp, num, den, ut;
  206.         double crit = 5;
  207.         if (u < 0.000001) return (- crit);
  208.         if (u > 0.999999) return (crit);
  209.         ut = ((u > 0.5) ? (1.0-u) : (u));
  210.         t = sqrt(- safelog(ut*ut));
  211.         num = (c2*t + c1) *t + c0;
  212.         den = ((d3*t + d2) *t + d1) *t + 1.0;
  213.         xp = t - num / den;
  214.         return (((u <= 0.5) ? (- xp) : xp));
  215. }
  216.  
  217. double erfc(double x)
  218. {
  219.     float t, z, ans;
  220.  
  221.     z = fabs(x);
  222.     t = 1.0 / (1.0+0.5*z);
  223.     ans = t*exp(- z*z - 1.26551223+ t*(1.00002368+ t*(0.37409196
  224.                  + t*(  0.09678418 +
  225.                  + t*(- 0.18628806 + t*(0.27886807+ t*(- 1.13520398
  226.                  + t*(  1.48851587
  227.                  + t*(- 0.82215223 + t*0.17087277)))))))));
  228.     ans = (x >= 0.0) ? ans : 2.0-ans;
  229.     return (1.0 - ans);
  230. }
  231. double probnorm(double x)
  232. {
  233.     return ((1.0+erfc((x) / 1414213562373095E-15)) * 0.5);
  234. }
  235.  
  236. /**************************** normal distribution functions */
  237. double cinv0(double p, double v)
  238.     /* approximate inverse chisquare */
  239.     {
  240.         double xp, chip;
  241.         
  242.         xp = probnormi(p);
  243.         chip = 1.0 - 2.0 / (9.0*v) + xp*sqrt(2.0 / (9.0*v));
  244.         chip = v*chip*chip*chip;
  245.         if (chip < 0) chip = 2.0*p*exp(2.0*gammln((0.5*v + 1.0)) / v);
  246.         return chip;
  247. }
  248.  
  249.  
  250. #define JMAX 75
  251.  
  252. double rtbis(double x1, double x2, double xacc, double v1,
  253.     double v2, double lam, double p,
  254.     double (*func) (double x1, double v1, double v2, double lam))
  255. {
  256.     int j;
  257.     double dx, f, fmid, xmid, rtb;
  258.     
  259.     f = (*func) (x1, v1, v2, lam) - p;
  260.     fmid = (*func) (x2, v1, v2, lam) - p;
  261.     if (f*fmid >= 0.0)
  262.         nrerror("RTBIS: Root must be bracketed for bisection");
  263.     rtb = f < 0.0 ? (dx = x2 - x1, x1) : (dx = x1 - x2, x2);
  264.     for (j = 1; j <= JMAX; j++) {
  265.         fmid = (*func) ((xmid = rtb + (dx *= 0.5)), v1, v2, lam) - p;
  266.         if (fmid <= 0.0) rtb = xmid;
  267.         if (fabs(dx) < xacc || fmid == 0.0) return rtb;
  268.     }
  269.     nrerror("RTBIS: Too many bisections");
  270. }
  271.  
  272. #undef JMAX
  273.  
  274.  
  275. #define FACTOR 1.6
  276. #define NTRY 75
  277.  
  278. int zbrac(double *x1, double *x2, double p, double v1,
  279.     double v2, double lam,
  280.     double (*func) (double x1, double v1, double v2, double lam))
  281. {
  282.     int j;
  283.     double f1, f2;
  284.     
  285.     if (*x1 == *x2) nrerror("ZBRAC: Bad initial range");
  286.     f1 = (*func) (*x1, v1, v2, lam) - p;
  287.     f2 = (*func) (*x2, v1, v2, lam) - p;
  288.     for (j = 1; j <= NTRY; j++) {
  289.         if (f1*f2 < 0.0) return 1;
  290.         if (fabs(f1) < fabs(f2))
  291.             f1 = (*func) ((*x1 += FACTOR*(*x1 - *x2)), v1, v2, lam) - p;
  292.         else
  293.             f2 = (*func) ((*x2 += FACTOR*(*x2 - *x1)), v1, v2, lam) - p;
  294.     }
  295.     return 0;
  296. }
  297.  
  298. #undef FACTOR
  299. #undef NTRY
  300.  
  301. double ncfinv(double p, double v1, double v2, double lam)
  302. {
  303.     double brak1, brak2, ncfi;
  304.     double xacc = 0.00001;
  305.     
  306.     if ((p < 0.0) || (p >= 1.0) || (v1 <= 0.0) || (v2 <= 0.0) || (lam < 0.0)){ 
  307.         printf("p, df1, df2, lam: %f %f %f %f\n", p, v1, v2, lam);
  308.         nrerror("NCFINV: invalid parameter");
  309.     }
  310.     
  311.     if (p == 0.0) return 0.0;
  312.  
  313.     brak1 = 0.000001;
  314.     brak2 = 2.0*cinv0(p, v1) / v1;
  315.     
  316.     if (!zbrac(&brak1, &brak2, p, v1, v2, lam, ncf))
  317.         nrerror("NCFINV: braketts not found");
  318.     
  319.     ncfi = rtbis(brak1, brak2, xacc, v1, v2, lam, p, ncf);
  320.     
  321.     return (ncfi);
  322. }
  323. double probfi(double p, double n1, double n2)
  324. {
  325.     if ((p < 0.0) || (p >= 1.0) || (n1 <= 0.0) || (n2 <= 0.0)) {
  326.         printf("p, df1, df2, lam: %f %f %f\n", p, n1, n2);
  327.         nrerror("PROBFI: invalid parameter");
  328.     }
  329.     return (ncfinv(p, n1, n2, 0.0));
  330. }
  331.  
  332. /******************************************     end ncfinv */
  333.  
  334. /************************************ non-central t distribution */
  335.  
  336.  
  337. double nct(double t, double df, double nc)
  338. {
  339.     double u = -0.2257913526, prob, f, g, q, r, s, a;
  340.     double tcnt = 0.0, p = 1.0, v = 0.0, w = 0.0, h, y, one = 1.0E7;
  341.     if ((df <= 0)) nrerror("NCT: non-positive parameter in nct");
  342.     if ( nc == 0 ) return probt(t, df);
  343.     
  344.     prob = (1.0+erfc((- nc) / 1414213562373095E-15)) * 0.5;
  345.     y = t*t;
  346.     f = (t < 0.0) ? - 1.0 : 1.0;
  347.     g = (nc < 0.0) ? - 1.0 : 1.0;
  348.     q = fabs(nc);
  349.     r = safelog(0.5) - q*q*0.5;
  350.     q = safelog(q);
  351.     s = f;
  352.     for (;;) {
  353.         tcnt++;
  354.         a = probf((y / tcnt), tcnt, df);
  355.         if (a <= 0.0) break;
  356.         h = exp(r + w + safelog(a));
  357.         prob += s*h;
  358.         if (prob >= (one * h)) break;
  359.         r += q;
  360.         p = -p;
  361.         if (p < 0.0) {
  362.             u -= safelog(tcnt);
  363.             w = u;
  364.             s = g;
  365.         } else {
  366.             v -= safelog(tcnt);
  367.             w = v;
  368.             s = f;
  369.         }
  370.     }     /* end for */
  371.     return prob;
  372. }     /* non central t distribution using continuous degrees of freedom */
  373.  
  374.  
  375. double probt(double t, double n1)
  376. {
  377.     double x;
  378.     if ((n1 <= 0)) nrerror("PROBT: non-positive parameter");
  379.     x = (t < 0) ? 0.5*(1.0- probf(t*t, 1.0, n1)) 
  380.                 : 0.5*(1.0+ probf(t*t, 1.0, n1));
  381.     return (x);
  382. }
  383. double probti(double p, double df)
  384. {
  385.     double one = 1.0, tx, q4;
  386.     if (df <= 0.0) nrerror("PROBTI: negative degrees of freedom");
  387.     if ((p < 0) || p > 1.0) nrerror("PROBTI: invalid probability p");
  388.     tx = 2.0*p - 1.0;
  389.     if (tx == 0.0) return tx;
  390.     q4 = probfi(fabs(tx), one, df);
  391.     tx = (tx < 0.0) ? - sqrt(q4) : sqrt(q4);
  392.     return tx;
  393. }
  394.  
  395. double ncti(double p, double df, double nc)
  396.     /* inverse non central t */
  397.     {
  398.         double tncx, q1, s, d1, r;
  399.         int cnt = 0;
  400.         if (df <= 0) nrerror("NCTI: negative degrees of freedom");
  401.         if ((p < 0) || p > 1.0) nrerror("NCTI: invalid probability p");
  402.         
  403.         tncx = probti(p, df) + nc;
  404.         q1 = nc / df;
  405.         if (df == 0.0) return tncx;
  406.         if (nc < 0.0) q1 = -q1;
  407.         L1 :
  408.         r = nct(tncx, df, nc) - p;
  409.         L2 :
  410.         cnt++;
  411.         tncx -=(r < 0.0) ? - q1 : q1;
  412.         s = nct(tncx, df, nc) - p;
  413.         if ((r*s == 0.0) || (cnt > 100)) return tncx;
  414.         if (r*s < 0.0) goto L3;
  415.         r = s;
  416.         goto L2;
  417.         L3 :
  418.         d1 = q1 / ((double) fabs((1.0 - (r / s))));
  419.         tncx += (r < 0.0) ? - d1 : d1;
  420.         q1 /= 64.0;
  421.         if ((d1 <= 1.0E-6) || (cnt > 100)) return tncx;
  422.         goto L1;
  423. }
  424. /*********************************** end non-central t */
  425. /*********************************** gamma, chi-square and nc */
  426.  
  427. #define ITMAX 100
  428. #define EPS 3.0e-7
  429.  
  430. void gcf(double *gammcf, double a, double x, double *gln)
  431. {
  432.     int n;
  433.     double gold = 0.0, g, fac = 1.0, b1 = 1.0;
  434.     double b0 = 0.0, anf, ana, an, a1, a0 = 1.0;
  435.     
  436.     *gln = gammln(a);
  437.     a1 = x;
  438.     for (n = 1; n <= ITMAX; n++) {
  439.         an = (double) n;
  440.         ana = an - a;
  441.         a0 = (a1 + a0*ana) *fac;
  442.         b0 = (b1 + b0*ana) *fac;
  443.         anf = an*fac;
  444.         a1 = x*a0 + anf*a1;
  445.         b1 = x*b0 + anf*b1;
  446.         if (a1) {
  447.             fac = 1.0 / a1;
  448.             g = b1*fac;
  449.             if (fabs((g - gold) / g) < EPS) {
  450.                 *gammcf = exp(- x + a*safelog(x) - (*gln)) *g;
  451.                 return;
  452.             }
  453.             gold = g;
  454.         }
  455.     }
  456.     nrerror("GCF: a too large, ITMAX too small");
  457. }
  458.  
  459. #undef ITMAX
  460. #undef EPS
  461.  
  462.  
  463. #define ITMAX 100
  464. #define EPS 3.0e-7
  465.  
  466. void gser(double *gamser, double a, double x, double *gln)
  467. {
  468.     int n;
  469.     double sum, del, ap;
  470.     
  471.     *gln = gammln(a);
  472.     if (x <= 0.0) {
  473.         if (x < 0.0) nrerror("GSER: x less than 0");
  474.         *gamser = 0.0;
  475.         return;
  476.     } else {
  477.         ap = a;
  478.         del = sum = 1.0 / a;
  479.         for (n = 1; n <= ITMAX; n++) {
  480.             ap += 1.0;
  481.             del *= x / ap;
  482.             sum += del;
  483.             if (fabs(del) < fabs(sum) *EPS) {
  484.                 *gamser = sum*exp(- x + a*safelog(x) - (*gln));
  485.                 return;
  486.             }
  487.         }
  488.         nrerror("GSER: a too large, ITMAX too small");
  489.         return;
  490.     }
  491. }
  492.  
  493. #undef ITMAX
  494. #undef EPS
  495.  
  496. double gammp(double x, double a, double b)
  497.     /* change from numerical recipes */
  498.     {
  499.         double gamser, gammcf, gln;
  500.         
  501.         x /= b;
  502.         if (x < 0.0 || a <= 0.0 || b <= 0.0)
  503.             nrerror("GAMMP: Invalid arguments");
  504.         if (x < (a + 1.0)) {
  505.             gser(&gamser, a, x, &gln);
  506.             return gamser;
  507.         } else {
  508.             gcf(&gammcf, a, x, &gln);
  509.             return 1.0- gammcf;
  510.         }
  511. }
  512. double ncgamma(double x, double a, double b, double lam)
  513.     /* algorithm AS 170 */
  514.     {
  515.         double eps = 1.0e-8, sum = 0, term = 0.0, c = 1.0, t = 0.0;
  516.  
  517.         if (x < 0.0 || a <= 0.0 || b <= 0.0 || lam < 0.0)
  518.             nrerror("NCGAMMA: Invalid arguments");
  519.         lam /= b;
  520.         sum = gammp(x, a, b);
  521.         do {
  522.             t++;
  523.             c *= lam / t;
  524.             a++;
  525.             term = c * gammp(x, a, b);
  526.             sum += term;
  527.         } while (term >= eps);
  528.         return (exp(safelog(sum) - lam));
  529. }
  530. double probchi(double x, double df)
  531. {
  532.     if (x < 0.0 || df <= 0.0)
  533.         nrerror("PROBCHI: Invalid arguments");
  534.     return (gammp(x, (0.5*df), 2.0));
  535. }
  536. double ncchi(double x, double df, double lam)
  537. {
  538.     if (x < 0.0 || df <= 0.0 || lam <= 0.0)
  539.         nrerror("NCCHI: Invalid arguments");
  540.     return (ncgamma(x, (0.5*df), 2.0, lam));
  541. }
  542.  
  543. double ncgammai(double p, double a, double b, double lam)
  544. {
  545.     double brak1, brak2, ncgamm;
  546.     double xacc = 0.000001;
  547.     
  548.     if ((p < 0.0) || (p >= 1.0) || (a <= 0.0) || (b <= 0.0) || (lam < 0.0)){            
  549.         printf("p, a, b, lam: %f %f %f %f\n", p, a, b, lam);
  550.         nrerror("NCGAMMAI: invalid parameter");
  551.     }
  552.     
  553.     if (p == 0.0) return 0.0;
  554.     
  555.     brak1 = 0.000001;
  556.     brak2 = 2.0*cinv0(p, a) / a;
  557.     
  558.     if (!zbrac(&brak1, &brak2, p, a, b, lam, ncgamma))
  559.         nrerror("ZBRAK: braketts not found");
  560.     
  561.     ncgamm = rtbis(brak1, brak2, xacc, a, b, lam, p, ncgamma);
  562.     
  563.     return (ncgamm);
  564. }
  565. double ncchii(double p, double a, double lam)
  566. {
  567.     
  568.     if ((p < 0.0) || (p >= 1.0) || (a <= 0.0) || (lam < 0.0)) {
  569.         printf("p, a, lam: %f %f %f\n", p, a, lam);
  570.         nrerror("NCCHII: invalid parameter");
  571.     }
  572.     return (ncgammai(p, (0.5*a), 2.0, lam));
  573. }
  574. double probchii(double p, double df)
  575. {
  576.     if ((p < 0.0) || (p >= 1.0) || (df <= 0.0)) {
  577.         printf("p, df: %f %f\n", p, df);
  578.         nrerror("PROBCHII: invalid parameter");
  579.     }
  580.     return (ncgammai(p, (0.5*df), 2.0, 0.0));
  581. }
  582.