home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 41 / IOPROG_41.ISO / soft / c++ / NUMCPP11.ZIP / normdist.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1999-05-09  |  6.1 KB  |  213 lines
  1. //===================================================================
  2. // normdist.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 "contdist.h"
  20. #include "domain.hpp"
  21. #include "mathx.h"
  22. #include "normdist.h"
  23. #include "numerror.h"
  24. #include "utility.hpp"
  25.  
  26. NUM_BEGIN
  27.  
  28. NUMERICS_EXPORT double chisqp(double x2, double v)
  29. {
  30.     if(!isPositive(v)){
  31.         throw Exception("chisqp", "Invalid v parameter");
  32.     } else if(x2 <= 0.0){
  33.         return 0.0;
  34.     }
  35.     return gammp(.5 * v, .5 * x2);
  36. }
  37.  
  38. NUMERICS_EXPORT double chisqv(double p, double v)
  39. {
  40.     double ch = 0.0;
  41.  
  42.     if(!isZeroOne(p)){
  43.         throw Exception("chisqv", "Invalid p parameter");
  44.     } else if(!isPositive(v)){
  45.         throw Exception("chisqv", "Invalid v parameter");
  46.     } else {
  47.         double g = gammln(.5*v), aa = .6931471806, e = 5e-7, r1, a, b, c, q, t, x;
  48.         double p1, p2, s1, s2, s3, s4, s5, s6, xx;
  49.         long i;
  50.         xx = .5 * v;
  51.         c = xx - 1.0;
  52.         if(v < -1.24 * log(p)){
  53.             ch = pow(p * xx * exp(g + xx * aa), 1.0 / xx);
  54.             if (ch < e) return ch;
  55.         } else {
  56.             if(v < .32){
  57.                 ch = .4;
  58.                 a = log(1.0 - p);
  59.                 do {
  60.                     q = ch;
  61.                     p1 = 1.0 + ch * (4.67 + ch);
  62.                     p2 = ch * (6.73 + ch * (6.66 + ch));
  63.                     t = -.5 + (4.67 + 2.0 * ch) / p1 - (6.73 + ch * (13.32 + 3.0 * ch)) / p2;
  64.                     ch -= (1.0 - exp(a + g + .5 * ch + c * aa) * p2 / p1) / t;
  65.                 } while(!areEqual(q / ch, 1.0));
  66.             } else {
  67.                 x = normalv(p, 0.0, 1.0);
  68.                 p1 = .222222 / v;
  69.                 r1 = x * sqrt(p1) + 1.0 - p1;
  70.                 ch = v * r1 * r1 * r1;
  71.                 if(ch > 2.2 * v + 6.0) ch = -2.0 * (log(1.0 - p) - c * log(.5 * ch) + g);
  72.             }
  73.         }
  74.         i = 0;
  75.         do {
  76.             q = ch;
  77.             p1 = .5 * ch;
  78.             p2 = p - gammp(xx, p1);
  79.             t = p2 * exp(xx * aa + g + p1 - c * log(ch));
  80.             b = t / ch;
  81.             a = .5 * t - b * c;
  82.             s1 = (210.0 + a * (140.0 + a * (105.0 + a * (84.0 + a * (70.0 + 60.0 * a))))) / 420.0;
  83.             s2 = (420.0 + a * (735.0 + a * (966.0 + a * (1141.0 + 1278.0 * a)))) / 2520.0;
  84.             s3 = (210.0 + a * (462.0 + a * (707.0 + 932.0 * a))) / 2520.0;
  85.             s4 = (252.0 + a * (672.0 + 1182.0 * a) + c * (294.0 + a * (889.0 + 1740.0 * a))) / 5040.0;
  86.             s5 = (84.0 + 2264.0 * a + c * (1175.0 + 606.0 * a)) / 2520.0;
  87.             s6 = (120. + c * (346.0 + 127.0 * c)) / 5040.0;
  88.             ch += t * (1.0 + .5 * t * s1 - b * c * (s1 - b * (s2 - b * (s3 - b * (s4 - b * (s5 - b * s6))))));
  89.         } while(!areEqual(q / ch, 1.0) && i++ < 20);
  90.     }
  91.     return ch;
  92. }
  93.  
  94. NUMERICS_EXPORT double fratiop(double f, double v1, double v2)
  95. {
  96.     if(!isPositive(v1)){
  97.         throw Exception("fratiop", "Invalid v1 parameter");
  98.     } else if(!isPositive(v2)){
  99.         throw Exception("fratiop", "Invalid v2 parameter");
  100.     } else if(isNonPositive(f)){
  101.         return 0.0;
  102.     }
  103.     return 1.0 - betai(v2 * .5, v1 * .5, v2 / (v2 + v1 * f));
  104. }
  105.  
  106. NUMERICS_EXPORT double fratiov(double p, double v1, double v2)
  107. {
  108.     double f0, f1, f2;
  109.     long i = 0;
  110.     
  111.     if(!isPositive(v1)){
  112.         throw Exception("fratiov", "Invalid v1 parameter");
  113.     } else if(!isPositive(v2)){
  114.         throw Exception("fratiov", "Invalid v2 parameter");
  115.     } else if(!isZeroOne(p)){
  116.         throw Exception("fratiov", "Invalid p parameter");
  117.     } else {
  118.         f0 = f1 = 5.0;
  119.         while(fratiop(f0,v1,v2) > p) f0 -= .5;
  120.         while(fratiop(f1,v1,v2) < p) f1 += .5;
  121.         if(f0 < 0.0) f0 = 0.0;
  122.         do {
  123.             f2 = f0 + (f1-f0)/2.0;
  124.             if(fratiop(f2,v1,v2) > p) f1 = f2;
  125.             else f0 = f2;
  126.         } while(!areEqual(f0, f1) && ++i < 100);
  127.     }
  128.     return f2;
  129. }
  130.  
  131. NUMERICS_EXPORT double normalp(double x, double m, double v)
  132. {
  133.     if(!isPositive(v)){
  134.         throw Exception("normalp", "Invalid v parameter");
  135.     } else if(areEqual(x, m)){
  136.         return .5;
  137.     }
  138.     return .5*(erff((x - m) / (sqrt(2.0 * v))) + 1.0);
  139. }
  140.  
  141. NUMERICS_EXPORT double normalv(double p, double m, double v)
  142. {
  143.     double q, r, z;
  144.     
  145.     q = p - .5;
  146.     if(fabs(q) <= .42){
  147.         r = q * q;
  148.         z = q * (((-25.4410604963 * r + 41.39119773534) * r - 18.61500062529) * r + 2.50662823884)
  149.             / ((((3.13082909833 * r -21.06224101826) * r + 23.08336743743) * r - 8.47351093090) * r + 1.0);
  150.         return z * sqrt(v) + m;
  151.     }
  152.     r = p;
  153.     if(q > 0.0) r = 1.0 - p;
  154.     if(r > 0.0){
  155.         r = sqrt(-log(r));
  156.         z = (((2.32121276858 * r + 4.85014127135) * r - 2.29796479134) * r - 2.78718931138)
  157.             / ((1.63706781897 * r + 3.54388924762) * r + 1.0);
  158.         if(q < 0.0) return -z * sqrt(v) + m;
  159.         return z * sqrt(v) + m;
  160.     }
  161.     return 0.0;
  162. }
  163.  
  164. NUMERICS_EXPORT double studtp(double t, double v)
  165. {
  166.     if(!isPositive(v)){
  167.         throw Exception("studtp", "Invalid v parameter");
  168.     } else if(areEqual(t, 0.0)){
  169.         return .5;
  170.     } else if(!isPositive(t)){
  171.         return 0.5 * betai(.5 * v, .5, v / (v + t * t));
  172.     }
  173.     return 1.0 - 0.5 * betai(.5 * v, .5, v / (v + t * t));
  174. }
  175.  
  176. NUMERICS_EXPORT double studtv(double p, double v)
  177. {
  178.     double t0, t1, t2, xx, sum, x;
  179.     int i;
  180.     
  181.     x = fabs(normalv(p, 0.0, 1.0));
  182.     xx = x * x;
  183.     sum = x + ((xx + 1.0) * x) / (4.0 * v);
  184.     sum += (((5.0 * xx + 16.0) * xx + 3.0) * x) / (96.0 * v * v);
  185.     sum += ((((3.0 * xx + 19.0) * xx + 17.0) * xx - 15.0) * x) / (384.0 * v * v * v);
  186.     sum += (((((79.0 * xx + 776.0) * xx + 1482.0) * xx - 1920.0) * xx - 945.0) * x) / (92160.0 * v * v * v * v);
  187.     if (p >= .5) t1 = sum;
  188.     else t1 = -sum;
  189.     t0 = t1;
  190.     while(studtp(t0,v) > p) t0 -= .5;
  191.     while(studtp(t1,v) < p) t1 += .5;
  192.     i = 0;
  193.     do {
  194.         t2 = t0 + (t1-t0)/2.0;
  195.         if(studtp(t2,v) > p) t1 = t2;
  196.         else t0 = t2;
  197.     } while(!areEqual(t0, t1) && ++i < 100);
  198.     return t2;
  199. }
  200.  
  201. NUM_END
  202.  
  203. //===================================================================
  204. // Revision History
  205. //
  206. // Version 1.0 - 08/28/1998 - New.
  207. // Version 1.1 - 04/10/1999 - Added Numerics namespace.
  208. //                            Use areEqual to test for floating-point
  209. //                            equality.
  210. //                            Added use of Exception class.
  211. //                            Added use of domain calls.
  212. //===================================================================
  213.