home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 41 / IOPROG_41.ISO / soft / c++ / NUMCPP11.ZIP / noncdist.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1999-05-09  |  2.3 KB  |  96 lines
  1. //===================================================================
  2. // noncdist.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 "domain.hpp"
  20. #include "mathx.h"
  21. #include "noncdist.h"
  22. #include "normdist.h"
  23. #include "numerror.h"
  24.  
  25. NUM_BEGIN
  26.  
  27. static const double R2DIVPI = sqrt(2.0 / NUMERICS_PI);
  28. static const double LOGRPI  = log(sqrt(NUMERICS_PI));
  29.  
  30. NUMERICS_EXPORT double ncstudtp(double x, double df, double delta)
  31. {
  32.     double tnc = 0.0;
  33.     double del = delta;
  34.     bool negdel = false;
  35.     
  36.     if(!isPositive(df)){
  37.         throw Exception("ncstudtp", "Invalid df");
  38.     } else {
  39.         if(x < 0.0){
  40.             negdel = true;
  41.             del *= -1;
  42.         }
  43.     
  44.         int en = 1;
  45.         double xx = x*x/(x*x+df);
  46.     
  47.         if(xx > 0.0){
  48.             double lambda = del*del;
  49.             double p = .5*exp(-.5*lambda);
  50.             double q = R2DIVPI * p * del;
  51.             double s = .5 - p, a = .5, b = .5*df;
  52.             double rxb = pow(1.0-xx, b);
  53.             double albeta = LOGRPI + gammln(b) - gammln(a+b);
  54.             double xodd = betai(a, b, xx);
  55.             double godd = 2.0*rxb*exp(a*log(x)-albeta);
  56.             double xeven = 1.0-rxb, geven = b*xx*rxb;
  57.             double errbd;
  58.             tnc = p * xodd + q * xeven;
  59.             do {
  60.                 a += 1.0;
  61.                 xodd -= godd;
  62.                 xeven -= geven;
  63.                 godd *= (xx*(a+b-1.0)/a);
  64.                 geven *= (xx*(a+b-.5)/(a+.5));
  65.                 p *= (lambda/(2.0*en));
  66.                 q *= (lambda/(2.0*en+1.0));
  67.                 s -= p;
  68.                 ++en;
  69.                 tnc += p*xodd + q*xeven;
  70.                 errbd = 2.0 * s * (xodd-godd);
  71.             } while(errbd > NUMERICS_MAX_ERROR && en <= NUMERICS_ITMAX);
  72.         }
  73.     
  74.         if(en > NUMERICS_ITMAX){
  75.             throw Exception("ncstudtp", "Iteration failed to converge");
  76.         } else {
  77.             tnc += (1.0 - normalp(del, 0.0, 1.0));
  78.             if(negdel){
  79.                 tnc = 1.0 - tnc;
  80.             }
  81.         }
  82.     }    
  83.     return tnc;
  84. }
  85.  
  86. NUM_END
  87.  
  88. //===================================================================
  89. // Revision History
  90. //
  91. // Version 1.0 - 08/28/1998 - New.
  92. // Version 1.1 - 04/10/1999 - Added Numerics namespace.
  93. //                            Added use of Exception class.
  94. //                            Added use of domain calls.
  95. //===================================================================
  96.