home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 41 / IOPROG_41.ISO / soft / c++ / NUMCPP11.ZIP / discdist.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1999-05-09  |  7.7 KB  |  308 lines
  1. //===================================================================
  2. // discdist.cpp
  3. //
  4. // Version 1.1
  5. //
  6. // Written by:
  7. //   Brent Worden
  8. //   WordenWare
  9. //   email:  Brent@Worden.org
  10. //
  11. // Copyright (c) 1998-1999 WordenWare
  12. //
  13. // Created:  August 28, 1998
  14. // Revised:  April 10, 1999
  15. //===================================================================
  16.  
  17. #include <cmath>
  18.  
  19. #include "discdist.h"
  20. #include "domain.hpp"
  21. #include "mathx.h"
  22. #include "normdist.h"
  23. #include "numerror.h"
  24.  
  25. NUM_BEGIN
  26.  
  27. NUMERICS_EXPORT double binomialp(int x, int n, double p)
  28. {
  29.     if(!isZeroOne(p)){
  30.         throw Exception("binomialp", "Invalid p parameter");
  31.     } else if(n < 1){
  32.         throw Exception("binomialp", "Invalid n parameter");
  33.     } else if(isNegative(x)){
  34.         return 0.0;
  35.     } else if(x >= n){
  36.         return 1.0;
  37.     }
  38.     return 1.0 - betai(x + 1, n - x, p);
  39. }
  40.  
  41. NUMERICS_EXPORT void binomialv(double pp, int n, double p, int *b0, int *b1)
  42. {
  43.     int b;
  44.  
  45.     if(!isZeroOne(p)){
  46.         throw Exception("binomialv", "Invalid p parameter");
  47.     } else if(!isZeroOne(pp)){
  48.         throw Exception("binomialv", "Invalid pp parameter");
  49.     } else if(n < 1){
  50.         throw Exception("binomialv", "Invalid n parameter");
  51.     } else {
  52.         b = int(n*p);
  53.         while(binomialp(b, n, p) > pp) b--;
  54.         while(binomialp(b, n, p) <= pp) b++;
  55.         if(b == 0){
  56.             *b0 = b;
  57.         } else *b0 = b-1;
  58.         *b1 = b;
  59.     }
  60. }
  61.  
  62. NUMERICS_EXPORT double geomp(int x, double p)
  63. {
  64.     if(!isZeroOne(p)){
  65.         throw Exception("geomp", "Invalid p parameter");
  66.     } else if(isNonPositive(x)){
  67.         return 0.0;
  68.     }
  69.     return (1.0 - pow(1.0 - p, x));
  70. }
  71.  
  72. NUMERICS_EXPORT void geomv(double pp, double p, int *x0, int *x1)
  73. {
  74.     double x;
  75.     
  76.     if(!isZeroOne(p)){
  77.         throw Exception("geomv", "Invalid p parameter");
  78.     } else if(!isZeroOne(pp)){
  79.         throw Exception("geomv", "Invalid pp parameter");
  80.     } else {
  81.         x = log(1.0 - pp) / log(1.0 - p);
  82.         *x0 = (int)floor(x);
  83.         if(*x0 == 0) *x0 = 1;
  84.         *x1 = (int)ceil(x);
  85.     }
  86. }
  87.  
  88. NUMERICS_EXPORT double hyperp(int x, int n, int N, int M)
  89. {
  90.     int i1, i2, i3, i4, i5, mnk1, i, j, k, l, m, nn, k1, n1;
  91.     double ret, r1, r2, mn, p, pt, sig;
  92.     bool dir;
  93.     
  94.     k = N + 1;
  95.     l = x + 1;
  96.     m = M + 1;
  97.     nn = n + 1;
  98.     dir = true;
  99.     
  100.     ret = 0.0;
  101.     if(nn < 1 || m < nn || k < 1 || k > m) {
  102.         return ret;
  103.     }
  104.     
  105.     if(l < 1 || k - l > m - nn) {
  106.         return ret;
  107.     }
  108.     
  109.     ret = 1.0;
  110.     if (l > nn || l > k) {
  111.         return ret;
  112.     }
  113.     if (k == 1 || k == m || nn == 1 || nn == m) {
  114.         return ret;
  115.     }
  116.     if(x == ((N < n) ? N : n)) {
  117.         return ret;
  118.     }
  119.     
  120.     p = (double)n / (double)(M - n);
  121.     i1 = N; i2 = M - N;
  122.     r1 = p; r2 = 1.0 / (double)p;
  123.     if(((i1 < i2) ? i1 : i2) > ((r1 > r2) ? r1 : r2) * 16.0 && M > 1000) {
  124.         mn = (double)(N * n) / (double)M;
  125.         sig = sqrt(mn * ((double)(M - n) / (double)M) * ((double)(M - N) / (double)(M - 1.0)));
  126.         r1 = (double)(x + .5 - mn) / sig;
  127.         ret = normalp(r1, 0.0, 1.0);
  128.     } else {
  129.         i1 = k - 1, i2 = m - k;
  130.         i3 = nn - 1, i4 = m - nn;
  131.         if(((i1 < i2) ? i1 : i2) > ((i3 < i4) ? i4 : i3)){
  132.             i = k;
  133.             k = nn;
  134.             nn = i;
  135.         }
  136.         if(m - k < k - 1) {
  137.             dir = ! dir;
  138.             l = nn - l + 1;
  139.             k = m - k + 1;
  140.         }
  141.         if(M > 600){
  142.             i1 = M - N;
  143.             i2 = M - n;
  144.             i3 = n - x;
  145.             i4 = N - x;
  146.             i5 = M - n - N + x;
  147.             p = factln(n) - factln(M) + factln(i1) + factln(N) +
  148.                 factln(i2) - factln(x) - factln(i3) - factln(i4) - factln(i5);
  149.             ret = 0.0;
  150.             if (p >= -88.0) {
  151.                 ret = exp(p);
  152.             }
  153.         } else {
  154.             i1 = l - 1;
  155.             for(i = 1; i <= i1; ++i) {
  156.                 ret *= (double)((k - i) * (nn - i)) / (double)((l - i) * (m - i));
  157.             }
  158.             if (l != k) {
  159.                 j = m - nn + l;
  160.                 i1 = k - 1;
  161.                 for (i = l; i <= i1; ++i) {
  162.                     ret *= (double) (j - i) / (double) (m - i);
  163.                 }
  164.             }
  165.         }
  166.         if (ret == 0.0) {
  167.             if (M <= 600) {
  168.                 i1 = M - n;
  169.                 i2 = n - x;
  170.                 i3 = N - x;
  171.                 i4 = M - n - N + x;
  172.                 i5 = M - N;
  173.                 p = factln(n) - factln(M) + factln(N) + factln(i1)
  174.                     - factln(x) - factln(i2) - factln(i3) -
  175.                     factln(i4) + factln(i5);
  176.             }
  177.             p += log(1.0e35);
  178.             if (p < -88.0) {
  179.                 if (x > (n * N + n + N + 1) / (M + 2)){
  180.                     ret = 1.0;
  181.                 }
  182.                 return ret;
  183.             } else {
  184.                 p = exp(p);
  185.             }
  186.         } else {
  187.             p = ret * 1.0e35;
  188.         }
  189.         pt = 0.0;
  190.         n1 = nn - l;
  191.         k1 = k - l;
  192.         mnk1 = m - nn - k1 + 1;
  193.         if (l <= k1) {
  194.             i1 = l - 1;
  195.             for (i = 1; i <= i1; ++i) {
  196.                 p *= (double)((l - i) * (mnk1 - i)) / (double)((n1 + i) * (k1 + i));
  197.                 pt += p;
  198.             }
  199.             if (p == 0.0) {
  200.                 //*ifault = 3;
  201.             }
  202.         } else {
  203.             dir = ! dir;
  204.             i1 = k1 - 1;
  205.             for (j = 0; j <= i1; ++j) {
  206.                 p *= (double)((n1 - j) * (k1 - j)) / (double)((l + j) * (mnk1 + j));
  207.                 pt += p;
  208.             }
  209.             if (p == (float)0.) {
  210.                 //*ifault = 3;
  211.             }
  212.         }
  213.         if (dir) {
  214.             ret += pt / 1.0e35;
  215.         } else {
  216.             ret = 1.0 - pt / 1.0e35;
  217.         }
  218.     }
  219.  
  220.     return ret;
  221. }
  222.  
  223. NUMERICS_EXPORT void hyperv(double p, int n, int N, int M, int *x0, int *x1)
  224. {
  225.     int x;
  226.     
  227.     if(!isZeroOne(p)){
  228.         throw Exception("hyperv", "Invalid p parameter");
  229.     } else {
  230.         x = (int)((double)(n * N) / (double)M);
  231.         while(hyperp(x, n, N, M) > p) x--;
  232.         while(hyperp(x, n, N, M) <= p) x++;
  233.         if(x <= 0){
  234.             *x0 = x;
  235.         } else {
  236.             *x0 = x-1;
  237.         }
  238.         *x1 = x;
  239.     }
  240. }
  241.  
  242. NUMERICS_EXPORT double negbnlp(int x, int r, double p)
  243. {
  244.     if(!isZeroOne(p)){
  245.         throw Exception("negbnlp", "Invalid p parameter");
  246.     } else if(isNonPositive(r)){
  247.         throw Exception("negbnlp", "Invalid r parameter");
  248.     } else if(isNegative(x)){
  249.         return 0.0;
  250.     }
  251.     return 1.0 - betai(x + 1.0, r - 1.0, p);
  252. }
  253.  
  254. NUMERICS_EXPORT void negbnlv(double p, int r, double pp, int *x0, int *x1)
  255. {
  256.     int x;
  257.     
  258.     if(!isZeroOne(p)){
  259.         throw Exception("negbnlv", "Invalid p parameter");
  260.     } else if(!isZeroOne(pp)){
  261.         throw Exception("negbnlv", "Invalid pp parameter");
  262.     } else if(isNonPositive(r)){
  263.         throw Exception("negbnlv", "Invalid r parameter");
  264.     } else {
  265.         x = (int)(r*(1.0-p)/p);
  266.         while(negbnlp(x, r, pp) > p) x--;
  267.         while(negbnlp(x, r, pp) <= p) x++;
  268.         if(x == 0){
  269.             *x0 = x;
  270.         } else *x0 = x-1;
  271.         *x1 = x;
  272.     }
  273. }
  274.  
  275. NUMERICS_EXPORT double poissonp(int x, double l)
  276. {
  277.     return gammq(x + 1.0, l);
  278. }
  279.  
  280. NUMERICS_EXPORT void poissonv(double p, double l, int *x0, int *x1)
  281. {
  282.     int x;
  283.     
  284.     if(!isZeroOne(p)){
  285.         throw Exception("poissonv", "Invalid p parameter");
  286.     } else if(!isPositive(l)){
  287.         throw Exception("poissonv", "Invalid l parameter");
  288.     } else {
  289.         x = (int)l;
  290.         while(poissonp(x, l) > p) x--;
  291.         while(poissonp(x, l) <= p) x++;
  292.         if(x == 0){
  293.             *x0 = x;
  294.         } else *x0 = x-1;
  295.         *x1 = x;
  296.     }
  297. }
  298.  
  299. NUM_END
  300.  
  301. //===================================================================
  302. // Revision History
  303. //
  304. // Version 1.0 - 08/28/1998 - New.
  305. // Version 1.1 - 04/10/1999 - Added Numerics namespace.
  306. //                            Added use of domain calls.
  307. //===================================================================
  308.