Changeset fb55369 in tol

Timestamp:

Aug 3, 2012, 9:04:16 AM (14 years ago)

Author:

Víctor de Buen <vdebuen@…>

Branches:

1941-docker-image-to-test, devel, master

Children:

3f203ae

Parents:

902a6ae

Message:

Refs #1590

git-svn-id: ​https://www.tol-project.org/svn/tolp/trunk@5005 be4846ae-f051-483a-831f-b64e7d08527f

Location:

tol/bmath/bprdist

Files:

• prd.cpp (modified) (14 diffs)
• tol_bprdist.h (modified) (5 diffs)

tol/bmath/bprdist/prd.cpp

@@ -2415 +2415 @@
 
 //--------------------------------------------------------------------
-BBinomialDist::BBinomialDist(BInt N, BDat p)
+BBinomialDist::BBinomialDist(BDat N, BDat p)
 
 /*! Creates an object aimed to work whith a binomial distribution. 
@@ -2430 +2430 @@
     logp_ = Log(p_);
     logq_ = Log(q_);
-    logGammaN_ = LogFact(N_);
+    logGammaN_ = LogGamma(1+N_);
     if((N_<=0)||(p_<0)||(p_>1))
     {
@@ -2451 +2451 @@
     BInt k = (BInt)Floor(x).Value();
     if((k<0)||(k>N_)) { return(0); }
-    return(Exp(logGammaN_-LogFact(N_-k)-LogFact(k)+k*logp_+(N_-k)*logq_));
+    return(Exp(logGammaN_-LogGamma(1+N_-k)-LogFact(k)+k*logp_+(N_-k)*logq_));
 }
 
@@ -2471 +2471 @@
         
         s = k;
-        xn = N_;
+        xn = N_.Value();
         pr = p_.Value();
         ompr = 1 - pr;
@@ -2494 +2494 @@
     p = x.Value();
     q = 1-p;
-    xn = (double) N_;
+    xn = N_.Value();
     pr = p_.Value();
     ompr = 1-pr;
@@ -2542 +2542 @@
 
 //--------------------------------------------------------------------
-BNegBinomialDist::BNegBinomialDist(BInt N, BDat P)
+BNegBinomialDist::BNegBinomialDist(BDat N, BDat P)
 
 /*! Creates an object aimed to work whith a negative binomial distribution. 
- *  A negative binomial distribution counts the number of failures
- *  before a fixed number of successes whith a fixed success probability.
+ *  A negative binomial distribution counts the number of successes
+ *  before a fixed number of failures whith a fixed success probability.
  * \param N Number of succeses.
  * \param P Success probability.
@@ -2576 +2576 @@
 
 /*! Returns NegBinomial(x,N_,P_) 
- *  = prob of getting x failures before the N_-th success
+ *  = prob of getting x succeses before the N_-th failures
  *
  *    ALGORITHM: NB(x,N_,P_) = P_ * DensBinomial(N_-1,k+N_-1,P_)
@@ -2583 +2583 @@
 //--------------------------------------------------------------------
 {
-    if(wrongParameter_ || x.IsUnknown()) { return(BDat::Unknown()); }
-    BInt k = (BInt)Floor(x).Value();
-    if((k<0)) { return(0); }
-    if((k==0)&&(N_==1)){return P_;}
-    else{
-        BBinomialDist b(k+N_-1,P_);
-        return(b.Dens(N_-1)*P_);    
-    } 
+  if(wrongParameter_ || x.IsUnknown()) { return(BDat::Unknown()); }
+  BInt k = (BInt)Floor(x).Value();
+  return(gsl_ran_negative_binomial_pdf(k,P_.Value(),N_.Value()));
 }
 
@@ -2601 +2596 @@
 //--------------------------------------------------------------------
 {
+  if(wrongParameter_ || x.IsUnknown()) { return(BDat::Unknown()); }
+  BInt k = (BInt)Floor(x).Value();
+//return(1-gsl_sf_beta_inc (k+1, N_.Value(), 1-P_.Value()));
+//return(gsl_cdf_negative_binomial_P(k,P_.Value(),N_.Value()));
+/* */
     double p, q, s, xn, pr, ompr, bound;
     int which, status; 
     s = x.Value();
-    xn = (double) N_;
+    xn = N_.Value();
     pr = P_.Value();
     ompr = 1-pr;
@@ -2613 +2613 @@
     else
         return(BDat::Unknown());
+/* */
 }
 
@@ -2623 +2624 @@
 {
   if(wrongParameter_ || x.IsUnknown()) { return(BDat::Unknown()); }
+//return(BDiscreteDist::Inverse(x,tolerance));
+
     double p, q, s, xn, pr, ompr, bound;
     int which, status; 
     p = x.Value();
     q = 1-p;
-    xn = (double) N_;
+    xn = N_.Value();
     pr = P_.Value();
     ompr = 1-pr;
@@ -2641 +2644 @@
         
         cumnbn(&x0,&xn,&pr,&ompr,&fx0,&q);
-        /* may be, there is an incremental way */
         cumnbn(&x1,&xn,&pr,&ompr,&fx1,&q);
         int K1 = 1;
         return fabs(fx1-p)<=spmpar(&K1) ? x1 : x0;
     }
+
 }
 
@@ -2778 +2781 @@
 //--------------------------------------------------------------------
 {
+  if(wrongParameter_ || p.IsUnknown()) { return(BDat::Unknown()); }
   return(BDiscreteDist::Inverse(p,tolerance));
 /*
@@ -2865 +2869 @@
 
 
+
+
+//--------------------------------------------------------------------
+BGenCountDist::BGenCountDist(BDat a, BDat v)
+
+/*! Creates an object aimed to deal with a Generic Counting distribution
+ *  depending just in average and variance parameters and that could be 
+ *    - Binomial:                 average>variance
+ *    - Poisson:                  average=variance
+ *    - Negative Binomial:        average<variance
+ */
+//--------------------------------------------------------------------
+:  BDiscreteDist(), a_(a), v_(v), 
+   cd_(NULL), m_(BDat::Unknown()), p_(BDat::Unknown())
+{
+  subType_ = Sign(a_-v_).Value();
+  m_ = a_*a_/Abs(a_-v_);
+  p_ = 1.0 - min(v_/a_,a_/v_);
+  if(subType_==1)
+  {
+    cd_ = new BBinomialDist(m_.Value(),p_);
+  }
+  else if(subType_==0)
+  {
+    cd_ = new BPoissonDist(a_);
+  }
+  if(subType_==-1)
+  {
+    cd_ = new BNegBinomialDist(m_.Value(),1-p_);
+  }
+}
+
+//--------------------------------------------------------------------
+ BGenCountDist::~BGenCountDist()
+//--------------------------------------------------------------------
+{
+  if(cd_) { delete cd_; }
+}
+
+//--------------------------------------------------------------------
+BDat BGenCountDist::Dens(BDat x)
+
+/*! Returns PDF
+ */
+//--------------------------------------------------------------------
+{
+  BDat d;
+  if(cd_) { d = cd_->Dens(x); }
+  return(d);
+}
+
+
+//--------------------------------------------------------------------
+BDat BGenCountDist::Dist(BDat x)
+
+/*! Returns CDF
+ */
+//--------------------------------------------------------------------
+{
+  BDat d;
+  if(cd_) { d = cd_->Dist(x); }
+  return(d);
+}
+
+//--------------------------------------------------------------------
+BDat BGenCountDist::Inverse(BDat p, BDat tolerance)
+
+/*! Returns the maximum x matching Poisson(x,fi) <= p
+ */
+//--------------------------------------------------------------------
+{
+  BDat i;
+  if(cd_) { i = cd_->Inverse(p,tolerance); }
+  return(i);
+}
+
 //--------------------------------------------------------------------
 BHypergeometricDist::BHypergeometricDist(BInt N, BInt P, BInt n)

tol/bmath/bprdist/tol_bprdist.h

@@ -873 +873 @@
 {
 private:
-  BInt N_;
+  BDat N_;
   BDat p_;
   BDat q_;
@@ -882 +882 @@
 public:
   // contructors and destructors:
-  BBinomialDist(BInt N, BDat p);
+  BBinomialDist(BDat N, BDat p);
 
     // Implementation: prd.cpp
@@ -901 +901 @@
 {
 private:
-  BInt N_;
+  BDat N_;
   BDat p_;
   BDat q_;
@@ -912 +912 @@
 public:
   // contructors and destructors:
-  BNegBinomialDist(BInt N, BDat P);
+  BNegBinomialDist(BDat N, BDat P);
 
   // Implementation: prd.cpp
@@ -959 +959 @@
 };
 
+//--------------------------------------------------------------------
+class TOL_API BGenCountDist : public BDiscreteDist
+
+/*! Generic counting probability distribution definition.
+ */
+//--------------------------------------------------------------------
+{
+private:
+  BDat a_; //average
+  BDat v_; //variance
+  int subType_;
+  BDiscreteDist* cd_;
+  BDat m_;
+  BDat p_;
+
+public:
+  // contructors and destructors:
+  BGenCountDist(BDat a, BDat v);
+ ~BGenCountDist();
+
+  // Implementation: prd.cpp
+  BDat Average  () { return(a_); }
+  BDat Varianze () { return(v_); }
+  BDat Dist     (BDat x);
+  BDat Dens     (BDat x);
+  BDat Inverse  (BDat prob, BDat tolerance = 0.1);
+};
 
 //--------------------------------------------------------------------