Probability.cpp

Go to the documentation of this file.

// Probability.cpp: implementation of the Probability class.
//
 
#include "Probability.h"
#include <boost/serialization/export.hpp>
#include <boost/archive/binary_oarchive.hpp>
#include <boost/archive/binary_iarchive.hpp>
BOOST_CLASS_EXPORT_IMPLEMENT(Probability)
using namespace std;
 
// Construction/Destruction
 
Probability::Probability() {
  probability = -1.0;
  alpha_parameter = 2.0;
  beta_parameter = 100.0;
  failure_rate = 0.0;
}
 
Probability::~Probability() {
 
}
void Probability::setProbability(double x) { probability = x; }
double Probability::getProbability() { return probability; }
void Probability::copyProbability(Probability &_pro) {
  probability = _pro.getProbability();
  alpha_parameter = _pro.getAlpha_Parameter();
  beta_parameter = _pro.getBeta_Parameter();
  failure_rate = _pro.getFailureRate();
}
void Probability::saveProbabilityToFile(ofstream &fout) {
  fout << probability << ' '
       << alpha_parameter << ' '
       << beta_parameter << '\n';
}
void Probability::retrieveProbabilityFromFile(ifstream &fin) {
  fin >> probability
      >> alpha_parameter
      >> beta_parameter;
}
void Probability::setAlpha_Parameter(double x) { alpha_parameter = x; }
double Probability::getAlpha_Parameter() { return alpha_parameter; }
void Probability::setBeta_Parameter(double x) { beta_parameter = x; }
double Probability::getBeta_Parameter() { return beta_parameter; }
 
double Probability::r8_epsilon()
//
//  Purpose:
//
//    R8_EPSILON returns the R8 round off unit.
//
//  Discussion:
//
//    R8_EPSILON is a number R which is a power of 2 with the property that,
//    to the precision of the computer's arithmetic,
//      1 < 1 + R
//    but 
//      1 = ( 1 + R / 2 )
//
//  Modified:
//
//    06 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double R8_EPSILON, the floating point round-off unit.
//
{
  double r;
 
  r = 1.0;
 
  while (1.0 < (double) (1.0 + r)) {
    r = r / 2.0;
  }
 
  return (2.0 * r);
}
 
double Probability::r8_huge()
//
//  Purpose:
//
//    R8_HUGE returns a "huge" R8.
//
//  Discussion:
//
//    HUGE_VAL is the largest representable legal real number, and is usually
//    defined in math.h, or sometimes in stdlib.h.
//
//  Modified:
//
//    08 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double R8_HUGE, a "huge" real value.
//
{
  return HUGE_VAL;
}
 
double Probability::beta(double x, double y)
//
//  Purpose:
//
//    BETA returns the value of the Beta function.
//
//  Formula:
//
//    BETA(X,Y) = ( GAMMA(X) * GAMMA(Y) ) / GAMMA(X+Y)
//
//  Properties:
//
//    BETA(X,Y) = BETA(Y,X).
//    BETA(X,Y) = Integral ( 0 <= T <= 1 ) T**(X-1) (1-T)**(Y-1) dT.
//    BETA(X,Y) = GAMMA(X) * GAMMA(Y) / GAMMA(X+Y)
//
//  Restrictions:
//
//    Both X and Y must be greater than 0.
//
//  Modified:
//
//    09 May 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, Y, the two parameters that define the Beta function.
//    X and Y must be greater than 0.
//
//    Output, double BETA, the value of the Beta function.
//
{
  if (x <= 0.0 || y <= 0.0) {
    cout << "\n";
    cout << "BETA - Fatal error!\n";
    cout << "  Both X and Y must be greater than 0.\n";
    exit(1);
  }
 
  return (exp(gamma_log(x) + gamma_log(y) - gamma_log(x + y)));
}
double Probability::beta_inc(double a, double b, double x)
//
//  Purpose:
//
//    BETA_INC returns the value of the incomplete Beta function.
//
//  Discussion:
//
//    This calculation requires an iteration.  In some cases, the iteration
//    may not converge rapidly, or may become inaccurate.
//
//  Formula:
//
//    BETA_INC(A,B,X)
//
//      =   Integral ( 0 <= T <= X ) T**(A-1) (1-T)**(B-1) dT
//        / Integral ( 0 <= T <= 1 ) T**(A-1) (1-T)**(B-1) dT
//
//      =   Integral ( 0 <= T <= X ) T**(A-1) (1-T)**(B-1) dT
//        / BETA(A,B)
//
//  Modified:
//
//    10 October 2004
//
//  Author:
//
//    Majumder and Bhattacharjee.
//
//    C++ translation by John Burkardt.
//
//  Reference:
//
//    Majumder, Bhattacharjee,
//    Algorithm AS63,
//    Applied Statistics,
//    1973, volume 22, number 3.
//
//  Parameters:
//
//    Input, double A, B, the parameters of the function.
//    0.0 < A,
//    0.0 < B.
//
//    Input, double X, the argument of the function.
//    Normally, 0.0 <= X <= 1.0.
//
//    Output, double BETA_INC, the value of the function.
//
{
  double cx = 0.0;
  int i = 0;
  int it = 0;
  int it_max = 1000;
  bool indx = false;
  int ns = 0;
  double pp = 0;
  double psq = 0.0;
  double qq = 0.0;
  double rx = 0.0;
  double temp = 0.0;
  double term = 0.0;
  double tol = 1.0E-07;
  double value = 0.0;
  double xx = 0.0;
 
  if (a <= 0.0) {
    cout << "\n";
    cout << "BETA_INC - Fatal error!\n";
    cout << "  A <= 0.\n";
    exit(1);
  }
 
  if (b <= 0.0) {
    cout << "\n";
    cout << "BETA_INC - Fatal error!\n";
    cout << "  B <= 0.\n";
    exit(1);
  }
 
  if (x <= 0.0) {
    value = 0.0;
    return value;
  } else if (1.0 <= x) {
    value = 1.0;
    return value;
  }
//
//  Change tail if necessary and determine S.
//
  psq = a + b;
 
  if (a < (a + b) * x) {
    xx = 1.0 - x;
    cx = x;
    pp = b;
    qq = a;
    indx = true;
  } else {
    xx = x;
    cx = 1.0 - x;
    pp = a;
    qq = b;
    indx = false;
  }
 
  term = 1.0;
  i = 1;
  value = 1.0;
 
  //ns = ( int ) ( qq + cx * ( a + b ) );
  ns = static_cast<int>( qq + cx * (a + b));
//
//  Use Soper's reduction formulas.
//
  rx = xx / cx;
 
  //temp = qq - ( double ) i;
  temp = qq - static_cast<double> (i);
  if (ns == 0) {
    rx = xx;
  }
 
  it = 0;
 
  for (;;) {
    it = it + 1;
 
    if (it_max < it) {
      cout << "\n";
      cout << "BETA_INC - Fatal error!\n";
      cout << "  Maximum number of iterations exceeded!\n";
      cout << "  IT_MAX = " << it_max << "\n";
      exit(1);
    }
 
    //term = term * temp * rx / ( pp + ( double ) ( i ) );
    term = term * temp * rx / (pp + static_cast<double>( i ));
    value = value + term;
    temp = fabs(term);
 
    if (temp <= tol && temp <= tol * value) {
      break;
    }
 
    i = i + 1;
    ns = ns - 1;
 
    if (0 <= ns) {
      //temp = qq - ( double ) i;
      temp = qq - static_cast<double> (i);
      if (ns == 0) {
        rx = xx;
      }
    } else {
      temp = psq;
      psq = psq + 1.0;
    }
  }
//
//  Finish calculation.
//
  value = value * exp(pp * log(xx)
                          + (qq - 1.0) * log(cx)) / (beta(a, b) * pp);
 
  if (indx) {
    value = 1.0 - value;
  }
 
  return value;
}
 
double Probability::beta_pdf(double x, double a, double b)
//
//
//  Purpose:
//
//    BETA_PDF evaluates the Beta PDF.
//
//  Formula:
//
//    PDF(A,B;X) = X**(A-1) * (1-X)**(B-1) / BETA(A,B).
//
//  Discussion:
//
//    A = B = 1 yields the Uniform distribution on [0,1].
//    A = B = 1/2 yields the Arcsin distribution.
//        B = 1 yields the power function distribution.
//    A = B -> Infinity tends to the Normal distribution.
//
//  Modified:
//
//    10 October 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, the argument of the PDF.
//    0.0 <= X <= 1.0.
//
//    Input, double A, B, the parameters of the PDF.
//    0.0 < A,
//    0.0 < B.
//
//    Output, double BETA_PDF, the value of the PDF.
//
{
  double pdf;
 
  if (x < 0.0 || 1.0 < x) {
    pdf = 0.0;
  } else {
    pdf = pow(x, (a - 1.0)) * pow((1.0 - x), (b - 1.0)) / beta(a, b);
  }
 
  return pdf;
}
 
double Probability::beta_cdf(double x, double a, double b)
//
//
//  Purpose:
//
//    BETA_CDF evaluates the Beta CDF.
//
//  Modified:
//
//    10 October 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, double X, the argument of the CDF.
//
//    Input, double A, B, the parameters of the PDF.
//    0.0 < A,
//    0.0 < B.
//
//    Output, double BETA_CDF, the value of the CDF.
//
{
  double cdf;
 
  if (x <= 0.0) {
    cdf = 0.0;
  } else if (x <= 1.0) {
    cdf = beta_inc(a, b, x);
  } else {
    cdf = 1.0;
  }
 
  return cdf;
}
 
double Probability::beta_cdf_inv(double _cdf, double _a, double _b)
// find the inverse of the cumulative beta distribution
{
  double p_low, p_high, p_est, F_est;
  double err, err_max;
  int max_count = 1000; // set maximum number of iterations
  p_low = 0.0;
  p_high = 1.0;
  err = 100.0;
  err_max = 1.0e-8; // set maximum error
  //p_est = a/(a+b);
  while (fabs(err) > err_max && max_count != 0) {
    p_est = 0.5 * (p_low + p_high);
    F_est = beta_cdf(p_est, _a, _b);
    err = F_est - _cdf;
    if (err <= 0.0) {
      p_low = fmax(0.0, p_low + 0.5 * fabs(p_est - p_low));
    } else {
      p_high = fmin(1.0, p_high - 0.5 * fabs(p_high - p_est));
    }
    max_count -= 1;
  }
  return p_est;
}
 
double Probability::gamma_log(double x)
//
//  Purpose:
//
//    GAMMA_LOG calculates the natural logarithm of GAMMA ( X ) for positive X.
//
//  Discussion:
//
//    The program uses rational functions that theoretically approximate
//    LOG(GAMMA(X)) to at least 18 significant decimal digits.  The
//    approximation for 12 < X is from Hart et al, while approximations
//    for X < 12.0 are similar to those in Cody and Hillstrom, but are 
//    unpublished.
//
//    The accuracy achieved depends on the arithmetic system, the compiler,
//    intrinsic functions, and proper selection of the machine-dependent
//    constants.
//
//  Modified:
//
//    12 May 2003
//
//  Author:
//
//    William Cody and L. Stoltz
//    Argonne National Laboratory
//
//    C++ translation by John Burkardt.
//
//  Reference:
//
//    William Cody, Kenneth Hillstrom,
//    Chebyshev Approximations for the Natural Logarithm of the Gamma Function,
//    Mathematics of Computation,
//    Volume 21, 1967, pages 198-203.
//
//    Kenneth Hillstrom,
//    ANL/AMD Program ANLC366S, DGAMMA/DLGAMA,
//    May 1969.
//
//    Hart, Ward Cheney, Charles Lawson, Maehly, 
//    Charles Mesztenyi, John Rice, Thacher, Witzgall,
//    Computer Approximations,
//    Wiley and sons, New York, 1968.
//
//  Parameters:
//
//    Input, double X, the argument of the Gamma function.  X must be positive.
//
//    Output, double GAMMA_LOG, the logarithm of the Gamma function of X.
//    If X <= 0.0, or if overflow would occur, the program returns the
//    value XINF, the largest representable floating point number.
//
//  Machine-dependent constants:
//
//  BETA   - radix for the floating-point representation.
//
//  MAXEXP - the smallest positive power of BETA that overflows.
//
//  XBIG   - largest argument for which LN(GAMMA(X)) is representable
//           in the machine, i.e., the solution to the equation
//             LN(GAMMA(XBIG)) = BETA**MAXEXP.
//
//  FRTBIG - Rough estimate of the fourth root of XBIG
//
//
//  Approximate values for some important machines are:
//
//                            BETA      MAXEXP         XBIG
//
//  CRAY-1        (S.P.)        2        8191       9.62E+2461
//  Cyber 180/855
//    under NOS   (S.P.)        2        1070       1.72E+319
//  IEEE (IBM/XT,
//    SUN, etc.)  (S.P.)        2         128       4.08E+36
//  IEEE (IBM/XT,
//    SUN, etc.)  (D.P.)        2        1024       2.55D+305
//  IBM 3033      (D.P.)       16          63       4.29D+73
//  VAX D-Format  (D.P.)        2         127       2.05D+36
//  VAX G-Format  (D.P.)        2        1023       1.28D+305
//
//
//                           FRTBIG
//
//  CRAY-1        (S.P.)   3.13E+615
//  Cyber 180/855
//    under NOS   (S.P.)   6.44E+79
//  IEEE (IBM/XT,
//    SUN, etc.)  (S.P.)   1.42E+9
//  IEEE (IBM/XT,
//    SUN, etc.)  (D.P.)   2.25D+76
//  IBM 3033      (D.P.)   2.56D+18
//  VAX D-Format  (D.P.)   1.20D+9
//  VAX G-Format  (D.P.)   1.89D+76
//
{
  double c[7] = {
      -1.910444077728E-03,
      8.4171387781295E-04,
      -5.952379913043012E-04,
      7.93650793500350248E-04,
      -2.777777777777681622553E-03,
      8.333333333333333331554247E-02,
      5.7083835261E-03};
  double corr;
  double d1 = -5.772156649015328605195174E-01;
  double d2 = 4.227843350984671393993777E-01;
  double d4 = 1.791759469228055000094023;
  double frtbig = 1.42E+09;
  int i;
  double p1[8] = {
      4.945235359296727046734888,
      2.018112620856775083915565E+02,
      2.290838373831346393026739E+03,
      1.131967205903380828685045E+04,
      2.855724635671635335736389E+04,
      3.848496228443793359990269E+04,
      2.637748787624195437963534E+04,
      7.225813979700288197698961E+03};
  double p2[8] = {
      4.974607845568932035012064,
      5.424138599891070494101986E+02,
      1.550693864978364947665077E+04,
      1.847932904445632425417223E+05,
      1.088204769468828767498470E+06,
      3.338152967987029735917223E+06,
      5.106661678927352456275255E+06,
      3.074109054850539556250927E+06};
  double p4[8] = {
      1.474502166059939948905062E+04,
      2.426813369486704502836312E+06,
      1.214755574045093227939592E+08,
      2.663432449630976949898078E+09,
      2.940378956634553899906876E+010,
      1.702665737765398868392998E+011,
      4.926125793377430887588120E+011,
      5.606251856223951465078242E+011};
  double pnt68 = 0.6796875;
  double q1[8] = {
      6.748212550303777196073036E+01,
      1.113332393857199323513008E+03,
      7.738757056935398733233834E+03,
      2.763987074403340708898585E+04,
      5.499310206226157329794414E+04,
      6.161122180066002127833352E+04,
      3.635127591501940507276287E+04,
      8.785536302431013170870835E+03};
  double q2[8] = {
      1.830328399370592604055942E+02,
      7.765049321445005871323047E+03,
      1.331903827966074194402448E+05,
      1.136705821321969608938755E+06,
      5.267964117437946917577538E+06,
      1.346701454311101692290052E+07,
      1.782736530353274213975932E+07,
      9.533095591844353613395747E+06};
  double q4[8] = {
      2.690530175870899333379843E+03,
      6.393885654300092398984238E+05,
      4.135599930241388052042842E+07,
      1.120872109616147941376570E+09,
      1.488613728678813811542398E+010,
      1.016803586272438228077304E+011,
      3.417476345507377132798597E+011,
      4.463158187419713286462081E+011};
  double res;
  double sqrtpi = 0.9189385332046727417803297;
  double xbig = 4.08E+36;
  double xden;
  double xm1;
  double xm2;
  double xm4;
  double xnum;
  double xsq;
//
//  Return immediately if the argument is out of range.
//
  if (x <= 0.0 || xbig < x) {
    return r8_huge();
  }
 
  if (x <= r8_epsilon()) {
    res = -log(x);
  } else if (x <= 1.5) {
    if (x < pnt68) {
      corr = -log(x);
      xm1 = x;
    } else {
      corr = 0.0;
      xm1 = (x - 0.5) - 0.5;
    }
 
    if (x <= 0.5 || pnt68 <= x) {
      xden = 1.0;
      xnum = 0.0;
 
      for (i = 0; i < 8; i++) {
        xnum = xnum * xm1 + p1[i];
        xden = xden * xm1 + q1[i];
      }
 
      res = corr + (xm1 * (d1 + xm1 * (xnum / xden)));
    } else {
      xm2 = (x - 0.5) - 0.5;
      xden = 1.0;
      xnum = 0.0;
      for (i = 0; i < 8; i++) {
        xnum = xnum * xm2 + p2[i];
        xden = xden * xm2 + q2[i];
      }
 
      res = corr + xm2 * (d2 + xm2 * (xnum / xden));
 
    }
  } else if (x <= 4.0) {
    xm2 = x - 2.0;
    xden = 1.0;
    xnum = 0.0;
    for (i = 0; i < 8; i++) {
      xnum = xnum * xm2 + p2[i];
      xden = xden * xm2 + q2[i];
    }
 
    res = xm2 * (d2 + xm2 * (xnum / xden));
  } else if (x <= 12.0) {
    xm4 = x - 4.0;
    xden = -1.0;
    xnum = 0.0;
    for (i = 0; i < 8; i++) {
      xnum = xnum * xm4 + p4[i];
      xden = xden * xm4 + q4[i];
    }
 
    res = d4 + xm4 * (xnum / xden);
  } else {
    res = 0.0;
 
    if (x <= frtbig) {
 
      res = c[6];
      xsq = x * x;
 
      for (i = 0; i < 6; i++) {
        res = res / xsq + c[i];
      }
 
    }
 
    res = res / x;
    corr = log(x);
    res = res + sqrtpi - 0.5 * corr;
    res = res + x * (corr - 1.0);
 
  }
 
  return res;
}
 
double Probability::calculateProbability() {
  double cdf, p;
/*  
    cdf = static_cast<double>(rand())/ static_cast<double>(RAND_MAX);
    //cout << cdf << '\n';
    if (cdf < 1.0e-6) 
    {
        cdf = 1.0e-6;
    }
    if (cdf > 0.999999)
    {
        cdf = 0.999999;
    }
    p = beta_cdf_inv (cdf, alpha_parameter, beta_parameter);
*/
  p = alpha_parameter / (alpha_parameter + beta_parameter);
  return p;
}
double Probability::getFailureRate() const {
  return failure_rate;
}
void Probability::setFailureRate(double failure_rate) {
  Probability::failure_rate = failure_rate;
}
double Probability::calculateSurvivalProbability(double deltaT) {
  double p;
  p = exp(-failure_rate*deltaT);
  return p;
}