home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 41 / IOPROG_41.ISO / soft / c++ / NUMCPP11.ZIP / contdist.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1999-05-17  |  9.8 KB  |  419 lines
  1. //===================================================================
  2. // contdist.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 "numerror.h"
  23. #include "utility.hpp"
  24.  
  25. NUM_BEGIN
  26.  
  27. NUMERICS_EXPORT double betap(double x, double a, double b)
  28. {
  29.     double ret = -1.0;
  30.     
  31.     if(!isPositive(a)){
  32.         throw Exception("betap", "Invalid a parameter");
  33.     } else if(!isPositive(b)){
  34.         throw Exception("betap", "Invalid b parameter");
  35.     } else if(!isPositive(x)){
  36.         ret = 0.0;
  37.     } else if(x >= 1.0){
  38.         ret = 1.0;
  39.     } else {
  40.         ret = betai(a, b, x);
  41.     }
  42.     
  43.     return ret;
  44. }
  45.  
  46. NUMERICS_EXPORT double betav(double p, double a, double b)
  47. {
  48.     double b0, b1, b2;
  49.     int i;
  50.     double ret = -1.0;
  51.     
  52.     if(!isZeroOne(p)){
  53.         throw Exception("betav", "Invalid p parameter");
  54.     } else if(!isPositive(a)){
  55.         throw Exception("betav", "Invalid a parameter");
  56.     } else if(!isPositive(b)){
  57.         throw Exception("betav", "Invalid b parameter");
  58.     } else {
  59.         b0 = b1 = .5;
  60.         while(betap(b0, a, b) > p){
  61.             b0 -= .1;
  62.         }
  63.         while(betap(b1, a, b) < p){
  64.             b1 += .1;
  65.         }
  66.         if(b0 < 0.0){
  67.             b0 = 0.0;
  68.         }
  69.         if(b1 > 1.0){
  70.             b1 = 1.0;
  71.         }
  72.         
  73.         i = 0;
  74.         do {
  75.             b2 = b0 + (b1 - b0) / 2.0;
  76.             if(betap(b2, a, b) > p){
  77.                 b1 = b2;
  78.             } else {
  79.                 b0 = b2;
  80.             }
  81.         } while(!areEqual(b0, b1) && ++i < NUMERICS_ITMAX);
  82.         
  83.         if(i >= NUMERICS_ITMAX){
  84.             throw Exception("betav", "Iteration failed to converge");
  85.         } else {
  86.             ret = b2;
  87.         }
  88.     }
  89.     
  90.     return ret;
  91. }
  92.  
  93. NUMERICS_EXPORT double cauchyp(double x, double m, double s)
  94. {
  95.     double ret = -1.0;
  96.     
  97.     if(!isPositive(s)){
  98.         throw Exception("cauchyp", "Invalid s parameter");
  99.     } else {
  100.         ret = atan((x - m) / s) / NUMERICS_PI + .5;
  101.     }
  102.     
  103.     return ret;
  104. }
  105.  
  106. NUMERICS_EXPORT double cauchyv(double p, double m, double s)
  107. {
  108.     double ret = 0.0;
  109.     
  110.     if(!isPositive(s)){
  111.         throw Exception("cauchyv", "Invalid s parameter");
  112.     } else if(!isZeroOne(p)){
  113.         throw Exception("cauchyv", "Invalid p parameter");
  114.     } else {
  115.         ret = m + s * tan(NUMERICS_PI * (p - .5));
  116.     }
  117.     
  118.     return ret;
  119. }
  120.  
  121. NUMERICS_EXPORT double dblexpp(double x, double m, double s)
  122. {
  123.     double ret = 0.0;
  124.     double sn = sign(x - m);
  125.     
  126.     if(!isPositive(s)){
  127.         throw Exception("dblexpp", "Invalid s parameter");
  128.     } else {
  129.         ret = (1.0 + sn * (1.0 - exp(-sn * (x - m) / s))) * .5;
  130.     }
  131.     
  132.     return ret;
  133. }
  134.  
  135. NUMERICS_EXPORT double dblexpv(double p, double m, double s)
  136. {
  137.     double sn = sign(2.0 * p - 1.0);
  138.     double ret = 0.0;
  139.     
  140.     if(!isPositive(s)){
  141.         throw Exception("dblexpv", "Invalid s parameter");
  142.     } else if(!isZeroOne(p)){
  143.         throw Exception("dblexpv", "Invalid p parameter");
  144.     } else {
  145.         ret = m - s * sn * log(1.0 - sn * (2.0 * p - 1.0));
  146.     }
  147.     
  148.     return ret;
  149. }
  150.  
  151. NUMERICS_EXPORT double expp(double x, double b)
  152. {
  153.     double ret = -1.0;
  154.     
  155.     if(!isPositive(b)){
  156.         throw Exception("expp", "Invalid b parameter");
  157.     } else if(!isPositive(x)){
  158.         ret = 0.0;
  159.     } else {
  160.         ret = 1.0 - exp(-x / b);
  161.     }
  162.     
  163.     return ret;
  164. }
  165.  
  166. NUMERICS_EXPORT double expv(double p, double b)
  167. {
  168.     double ret = -1.0;
  169.  
  170.     if(!isZeroOne(p)){
  171.         throw Exception("expv", "Invalid p parameter");
  172.     } else if(!isPositive(b)){
  173.         throw Exception("expv", "Invalid b parameter");
  174.     } else {
  175.         ret = -log(1 - p) * b;
  176.     }
  177.     
  178.     return ret;
  179. }
  180.  
  181. NUMERICS_EXPORT double gammap(double x, double a, double b)
  182. {
  183.     double ret = -1.0;
  184.     
  185.     if(!isPositive(a)){
  186.         throw Exception("gammp", "Invalid a parameter");
  187.     } else if(!isPositive(b)){
  188.         throw Exception("gammp", "Invalid b parameter");
  189.     } else if(!isPositive(x)){
  190.         ret = 0.0;
  191.     } else {
  192.         ret = 1.0 - gammq(a, x / b);
  193.     }
  194.     
  195.     return ret;
  196. }
  197.  
  198. NUMERICS_EXPORT double gammav(double p, double a, double b)
  199. {
  200.     int i;
  201.     double g0, g1, g2;
  202.     double ret = -1.0;
  203.     
  204.     if(!isPositive(a)){
  205.         throw Exception("gammav", "Invalid a parameter");
  206.     } else if(!isPositive(b)){
  207.         throw Exception("gammav", "Invalid b parameter");
  208.     } else if(!isZeroOne(p)){
  209.         throw Exception("gammav", "Invalid p parameter");
  210.     } else {
  211.         g0 = g1 = a*b;
  212.         while(gammap(g0, a, b) > p){
  213.             g0 -= .5;
  214.         }
  215.         while(gammap(g1, a, b) < p){
  216.             g1 += .5;
  217.         }
  218.         if(g0 < 0.0){
  219.             g0 = 0.0;
  220.         }
  221.         i = 0;
  222.         do {
  223.             g2 = g0 + (g1-g0)/2.0;
  224.             if(gammap(g2,a,b) > p){
  225.                 g1 = g2;
  226.             } else {
  227.                 g0 = g2;
  228.             }
  229.         } while(!areEqual(g0, g1) && ++i < 100);
  230.         if(i >= NUMERICS_ITMAX){
  231.             throw Exception("gammav", "Iteration failed to converge");
  232.         } else {
  233.             ret = g2;
  234.         }
  235.     }
  236.     
  237.     return ret;
  238. }
  239.  
  240. NUMERICS_EXPORT double logisticp(double x, double m, double b)
  241. {
  242.     double ret = -1.0;
  243.     
  244.     if(!isPositive(b)){
  245.         throw Exception("logisticp", "Invalid b parameter");
  246.     } else {
  247.         ret = 1.0 / (1.0 + exp((m - x) / b));
  248.     }
  249.     
  250.     return ret;
  251. }
  252.  
  253. NUMERICS_EXPORT double logisticv(double p, double m, double b)
  254. {
  255.     double ret = -1.0;
  256.     
  257.     if(!isPositive(b)){
  258.         throw Exception("logisticp", "Invalid b parameter");
  259.     } else if(!isZeroOne(p)){
  260.         throw Exception("logisticv", "Invalid p parameter");
  261.     } else {
  262.         ret = m - b * log(1.0 / p - 1.0);
  263.     }
  264.     
  265.     return ret;
  266. }
  267.  
  268. NUMERICS_EXPORT double lognormp(double x, double m, double v)
  269. {
  270.     double ret = -1.0;
  271.     
  272.     if(!isPositive(v)){
  273.         throw Exception("lognormp", "Invalid v parameter");
  274.     } else if(!isPositive(x)){
  275.         ret = 0.0;
  276.     } else {
  277.         ret = .5 * (erff((log(x) - m) / (sqrt(2.0 * v))) + 1.0);
  278.     }
  279.     
  280.     return ret;
  281. }
  282.  
  283. NUMERICS_EXPORT double lognormv(double p, double m, double v)
  284. {
  285.     double l0, l1, l2;
  286.     int i;
  287.     double ret = -1.0;
  288.     
  289.     if(!isPositive(v)){
  290.         throw Exception("lognormv", "Invalid v parameter");
  291.     } else if(!isZeroOne(p)){
  292.         throw Exception("lognormv", "Invalid p parameter");
  293.     } else {
  294.         l0 = l1 = exp(m+v*v*.5);
  295.         while(lognormp(l0, m, v) > p){
  296.             l0 -= .5;
  297.         }
  298.         while(lognormp(l1, m, v) < p){
  299.             l1 += .5;
  300.         }
  301.         if(l0 < 0.0){
  302.             l0 = 0.0;
  303.         }
  304.         i = 0;
  305.         do {
  306.             l2 = l0 + (l1-l0)/2.0;
  307.             if(lognormp(l2,m,v) > p){
  308.                 l1 = l2;
  309.             } else {
  310.                 l0 = l2;
  311.             }
  312.         } while(!areEqual(l0, l1) && ++i < 100);
  313.         if(i >= NUMERICS_ITMAX){
  314.             throw Exception("lognormv", "Iteration failed to converge");
  315.         } else {
  316.             ret = l2;
  317.         }
  318.     }
  319.     
  320.     return l2;
  321. }
  322.  
  323. NUMERICS_EXPORT double rayleighp(double x, double s)
  324. {
  325.     return 1.0 - exp(-x * x/ (2.0 * s * s));
  326. }
  327.  
  328. NUMERICS_EXPORT double rayleighv(double p, double s)
  329. {
  330.     double ret = -1.0;
  331.     
  332.     if(!isZeroOne(p)){
  333.         throw Exception("rayleighv", "Invalid p parameter");
  334.     } else {
  335.         ret = s * sqrt(log( 1.0 / ((1.0 - p) * (1.0 - p))));
  336.     }
  337.     
  338.     return ret;
  339. }
  340.  
  341. NUMERICS_EXPORT double uniformp(double u, double a, double b)
  342. {
  343.     double ret = -1.0;
  344.     
  345.     if(a > b){
  346.         throw Exception("uniformp", "Invalid a and b parameters");
  347.     } else if(u <= a){
  348.         ret = 0.0;
  349.     } else if(u > b){
  350.         ret = 1.0;
  351.     } else {
  352.         ret = (u - a) / (b - a);
  353.     }
  354.     
  355.     return ret;
  356. }
  357.  
  358. NUMERICS_EXPORT double uniformv(double p, double a, double b)
  359. {
  360.     double ret = -1.0;
  361.     
  362.     if(!isZeroOne(p)){
  363.         throw Exception("uniformv", "Invalid p parameter");
  364.     } else if(a > b){
  365.         throw Exception("uniformv", "Invalid a and b parameters");
  366.     } else {
  367.         ret = p * (b - a) + a;
  368.     }
  369.     
  370.     return ret;
  371. }
  372.  
  373. NUMERICS_EXPORT double weibullp(double x, double g, double b)
  374. {
  375.     double ret = -1.0;
  376.     
  377.     if(!isPositive(g)){
  378.         throw Exception("weibullp", "Invalid g parameter");
  379.     } else if(!isPositive(b)){
  380.         throw Exception("weibullp", "Invalid b parameter");
  381.     } else if(!isPositive(x)){
  382.         ret = 0.0;
  383.     } else {
  384.         ret = 1.0 - exp(-pow(x / b, g));
  385.     }
  386.     
  387.     return ret;
  388. }
  389.  
  390. NUMERICS_EXPORT double weibullv(double p, double g, double b)
  391. {
  392.     double ret = -1.0;
  393.     
  394.     if(!isPositive(g)){
  395.         throw Exception("weibullv", "Invalid g parameter");
  396.     } else if(!isPositive(b)){
  397.         throw Exception("weibullv", "Invalid b parameter");
  398.     } else if(!isZeroOne(p)){
  399.         throw Exception("weibullv", "Invalid p parameter");
  400.     } else {
  401.         ret = b * pow(-log(1.0 - p), (1.0 / g));
  402.     }
  403.     
  404.     return ret;
  405. }
  406.  
  407. NUM_END
  408.  
  409. //===================================================================
  410. // Revision History
  411. //
  412. // Version 1.0 - 08/28/1998 - New.
  413. // Version 1.1 - 04/10/1999 - Use areEqual to test for floating-point
  414. //                            equality. 
  415. //                            Added Numerics namespace.
  416. //                            Added use of Exception class.
  417. //                            Added use of domain calls.
  418. //===================================================================
  419.