//+------------------------------------------------------------------+
//|                                           Distribution_class.mqh |
//|                                                           denkir |
//|                                                 denkir@gmail.com |
//+------------------------------------------------------------------+
#property copyright "denkir"
#property link      "denkir@gmail.com"

// definitions
#define ITMAX 100         //maximum allowed number of iterations
#define EPS DBL_EPSILON   //1.0+DBL_EPSILON != 1.0
#define FPMIN DBL_MIN/EPS //numeric_limits<double>::min()/EPS;
//+------------------------------------------------------------------+
//|                Error Function class definition                   |
//+------------------------------------------------------------------+
class CErf
  {
public:
   int               ncof;    //coefficient array size
   double            cof[28]; //array of Chebyshev coefficients
   //+------------------------------------------------------------------+
   //| CErf class constructor                                           |
   //+------------------------------------------------------------------+
   void CErf()
     {
      int Ncof=28;
      double Cof[28]=//Chebyshev coefficients
        {
         -1.3026537197817094,6.4196979235649026e-1,
         1.9476473204185836e-2,-9.561514786808631e-3,-9.46595344482036e-4,
         3.66839497852761e-4,4.2523324806907e-5,-2.0278578112534e-5,
         -1.624290004647e-6,1.303655835580e-6,1.5626441722e-8,-8.5238095915e-8,
         6.529054439e-9,5.059343495e-9,-9.91364156e-10,-2.27365122e-10,
         9.6467911e-11, 2.394038e-12,-6.886027e-12,8.94487e-13, 3.13092e-13,
         -1.12708e-13,3.81e-16,7.106e-15,-1.523e-15,-9.4e-17,1.21e-16,-2.8e-17
        };
      setCErf(Ncof,Cof);
     };
   //+------------------------------------------------------------------+
   //| Set-method for ncof                                              |
   //+------------------------------------------------------------------+
   void setCErf(int Ncof,double &Cof[])
     {
      ncof=Ncof;
      ArrayCopy(cof,Cof);
     };
   //+------------------------------------------------------------------+
   //| CErf class destructor                                            |
   //+------------------------------------------------------------------+
   void ~CErf(){};
   //+------------------------------------------------------------------+
   //| Error function                                                   |
   //+------------------------------------------------------------------+
   double erf(double x)
     {
      if(x>=0.0) return 1.0-erfccheb(x);
      else return erfccheb(-x)-1.0;
     }
   //+------------------------------------------------------------------+
   //| Complementary error function                                     |
   //+------------------------------------------------------------------+
   double erfc(double x)
     {
      if(x>=0.0) return erfccheb(x);
      else return 2.0-erfccheb(-x);
     }
   //+------------------------------------------------------------------+
   //| Chebyshev approximations of the error function                   |
   //+------------------------------------------------------------------+
   double erfccheb(double z)
     {
      int j;
      double t,ty,tmp,d=0.0,dd=0.0;
      if(z<0.) Alert("erfccheb requires nonnegative argument!");
      t=2.0/(2.0+z);
      ty=4.0*t-2.0;
      for(j=ncof-1;j>0;j--)
        {
         tmp=d;
         d=ty*d-dd+cof[j];
         dd=tmp;
        }
      return t*exp(-z*z+0.5*(cof[0]+ty*d)-dd);
     }
   //+------------------------------------------------------------------+
   //| Inverse complementary error function                             |
   //+------------------------------------------------------------------+
   double inverfc(double p)
     {
      double x,err,t,pp;
      if(p >= 2.0) return -100.0;
      if(p <= 0.0) return 100.0;
      pp=(p<1.0)? p : 2.0-p;
      t = sqrt(-2.*log(pp/2.0));
      x = -0.70711*((2.30753+t*0.27061)/(1.0+t*(0.99229+t*0.04481)) - t);
      for(int j=0;j<2;j++)
        {
         err=erfc(x)-pp;
         x+=err/(M_2_SQRTPI*exp(-pow(x,2))-x*err);
        }
      return(p<1.0? x : -x);
     }
   //+------------------------------------------------------------------+
   //| Inverse error function                                           |
   //+------------------------------------------------------------------+
   double inverf(double p)
     {return inverfc(1.0-p);}
  };
//+------------------------------------------------------------------+
double erfcc(const double x)
/* 
 complementary error function erfc(x) with
 relative error of 1.2 * 10^(-7)               
*/
  {
   double t,z=fabs(x),ans;
   t=2./(2.0+z);
   ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
             t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
             t*(-0.82215223+t*0.17087277)))))))));
   return(x>=0.0 ? ans : 2.0-ans);
  }
//+------------------------------------------------------------------+
//|                Normal Distribution class definition              |
//+------------------------------------------------------------------+
class CNormaldist : CErf // Erf class inheritance
  {
public:
   double            mu,//mean parameter (?)
   sig;                 //variance parameter (?)
   //+------------------------------------------------------------------+
   //| CNormaldist class constructor                                    |
   //+------------------------------------------------------------------+
   void  CNormaldist()
     {
      mu=0.0;sig=1.0; //parameters ? and ? by default
      if(sig<=0.) Alert("bad sig in Normal Distribution!");
     }
   //+------------------------------------------------------------------+
   //| Probability density function                                     |
   //| Probability density function (pdf)                               |
   //+------------------------------------------------------------------+
   double pdf(double x)
     {
      return(0.398942280401432678/sig)*exp(-0.5*pow((x-mu)/sig,2));
     }
   //+------------------------------------------------------------------+
   //| Cumulative distribution function                                 |
   //| Cumulative distribution function (cdf)                           |
   //+------------------------------------------------------------------+
   double cdf(double x)
     {
      return 0.5*erfc(-0.707106781186547524*(x-mu)/sig);
     }
   //+------------------------------------------------------------------+
   //| Inverse cumulative distribution function (quantile function)     |
   //| Inverse cumulative distribution function (invcdf)                |
   //+------------------------------------------------------------------+
   double invcdf(double p)
     {
      if(!(p>0. && p<1.))
         Alert("bad p in Normal Distribution!");
      return -1.41421356237309505*sig*inverfc(2.*p)+mu;
     }
   //+------------------------------------------------------------------+
   //| Reliability (survival) function                                  |
   //| Survival function (sf)                                           |
   //+------------------------------------------------------------------+
   double sf(double x)
     {
      return 1-cdf(x);
     }
  };
//+------------------------------------------------------------------+
//|                Lognormal Distribution class definition           |
//+------------------------------------------------------------------+
class CLognormaldist : CErf // Erf class inheritance
  {
public:
   double            mu,//location parameter (?)
   sig;                 //scale parameter (?)
   //+------------------------------------------------------------------+
   //| CLognormaldist class constructor                                 |
   //+------------------------------------------------------------------+
   void  CLognormaldist()
     {
      mu=0.0;sig=1.0; //parameters ? and ? by default
      if(sig<=0.) Alert("bad sig in Lognormal Distribution!");
     }
   //+------------------------------------------------------------------+
   //| Probability density function                                     |
   //| Probability density function (pdf)                               |
   //+------------------------------------------------------------------+
   double pdf(double x)
     {
      if(x<0.) Alert("bad x in Lognormal Distribution!");
      if(x==0.) return 0.;
      return(0.398942280401432678/(sig*x))*exp(-0.5*pow((log(x)-mu)/sig,2));
     }
   //+------------------------------------------------------------------+
   //| Cumulative distribution function                                 |
   //| Cumulative distribution function (cdf)                           |
   //+------------------------------------------------------------------+
   double cdf(double x)
     {
      if(x<0.) Alert("bad x in Lognormal Distribution!");
      if(x==0.) return 0.;
      return 0.5*erfc(-0.707106781186547524*(log(x)-mu)/sig);
     }
   //+------------------------------------------------------------------+
   //| Inverse cumulative distribution function (quantile function)     |
   //| Inverse cumulative distribution function (invcdf)                |
   //+------------------------------------------------------------------+
   double invcdf(double p)
     {
      if(!(p>0. && p<1.))
         Alert("bad p in Lognormal Distribution!");
      return exp(-1.41421356237309505*sig*inverfc(2.*p)+mu);
     }
   //+------------------------------------------------------------------+
   //| Reliability (survival) function                                  |
   //| Survival function (sf)                                           |
   //+------------------------------------------------------------------+
   double sf(double x)
     {
      return 1-cdf(x);
     }
  };
//+------------------------------------------------------------------+
//|                Cauchy Distribution class definition              |
//+------------------------------------------------------------------+
class CCauchydist //
  {
public:
   double            mu,//location parameter (?)
   sig;                 //scale parameter (?)
   //+------------------------------------------------------------------+
   //| CCauchydist class constructor                                    |
   //+------------------------------------------------------------------+
   void  CCauchydist()
     {
      mu=0.0;sig=1.0; //parameters ? and ? by default
      if(sig<=0.) Alert("bad sig in Cauchy Distribution!");
     }
   //+------------------------------------------------------------------+
   //| Probability density function                                     |
   //| Probability density function (pdf)                               |
   //+------------------------------------------------------------------+
   double pdf(double x)
     {
      return 0.318309886183790671/(sig*(1.+pow((x-mu)/sig,2)));
     }
   //+------------------------------------------------------------------+
   //| Cumulative distribution function                                 |
   //| Cumulative distribution function (cdf)                           |
   //+------------------------------------------------------------------+
   double cdf(double x)
     {
      return 0.5+0.318309886183790671*atan2(x-mu,sig); //todo      
     }
   //+------------------------------------------------------------------+
   //| Inverse cumulative distribution function (quantile function)     |
   //| Inverse cumulative distribution function (invcdf)                |
   //+------------------------------------------------------------------+
   double invcdf(double p)
     {
      if(!(p>0. && p<1.))
         Alert("bad p in Cauchy Distribution!");
      return mu+sig*tan(M_PI*(p-0.5));
     }
   //+------------------------------------------------------------------+
   //| Reliability (survival) function                                  |
   //| Survival function (sf)                                           |
   //+------------------------------------------------------------------+
   double sf(double x)
     {
      return 1-cdf(x);
     }
  };
//+------------------------------------------------------------------+
double atan2(double y,double x)
/*
 Returns the principal value of the arc tangent of y/x,
 expressed in radians. To compute the value, the function 
 uses the sign of both arguments to determine the quadrant.
 y - double value representing an y-coordinate.
 x - double value representing an x-coordinate. 
*/
  {
   double a;
   if(fabs(x)>fabs(y))
      a=atan(y/x);
   else
     {
      a=atan(x/y); // pi/4 <= a <= pi/4
      if(a<0.)
         a=-1.*M_PI_2-a; //a is negative, so we're adding
      else
         a=M_PI_2-a;
     }
   if(x<0.)
     {
      if(y<0.)
         a=a-M_PI;
      else
         a=a+M_PI;
     }
   return a;
  }
//+------------------------------------------------------------------+
//|        Hyperbolic Secant Distribution class definition           |
//+------------------------------------------------------------------+
class CHypersecdist //
  {
public:
   double            mu,// location parameter (?)
   sig;                 //scale parameter (?)
   //+------------------------------------------------------------------+
   //| CHypersecdist class constructor                                  |
   //+------------------------------------------------------------------+
   void  CHypersecdist()
     {
      mu=0.0;sig=1.0; //parameters ? and ? by default
      if(sig<=0.) Alert("bad sig in Hyperbolic Secant Distribution!");
     }
   //+------------------------------------------------------------------+
   //| Probability density function                                     |
   //| Probability density function (pdf)                               |
   //+------------------------------------------------------------------+
   double pdf(double x)
     {
      return sech((M_PI*(x-mu))/(2.*sig))/(2.*sig);
     }
   //+------------------------------------------------------------------+
   //| Cumulative distribution function                                 |
   //| Cumulative distribution function (cdf)                           |
   //+------------------------------------------------------------------+
   double cdf(double x)
     {
      return 2/M_PI*atan(exp((M_PI*(x-mu)/(2*sig))));
     }
   //+------------------------------------------------------------------+
   //| Inverse cumulative distribution function (quantile function)     |
   //| Inverse cumulative distribution function (invcdf)                |
   //+------------------------------------------------------------------+
   double invcdf(double p)
     {
      if(!(p>0. && p<1.))
         Alert("bad p in Hyperbolic Secant Distribution!");
      return(mu+(2.0*sig/M_PI*log(tan(M_PI/2.0*p))));
     }
   //+------------------------------------------------------------------+
   //| Reliability (survival) function                                  |
   //| Survival function (sf)                                           |
   //+------------------------------------------------------------------+
   double sf(double x)
     {
      return 1-cdf(x);
     }
  };
//+------------------------------------------------------------------+
//|              Hyperbolic Secant Function                          |
//+------------------------------------------------------------------+
double sech(double x)
  {
   double k=pow(M_E,x)+pow(M_E,-1.*x);
   return 2./k;
  }
//+------------------------------------------------------------------+
//|                 Gauleg18 class definition                        |
//+------------------------------------------------------------------+
class CGauleg18
//CGauleg18 class provides coefficients for Gauss-Legendre quadrature
  {
protected:
   int               ngau;
   double            y[18];
   double            w[18];
public:
   //+------------------------------------------------------------------+
   //| CGauleg18 class constructor                                      |
   //+------------------------------------------------------------------+
   void CGauleg18()
     {
      int Ngau=18;
      double Y[18]=
        {
         0.0021695375159141994,
         0.011413521097787704,0.027972308950302116,0.051727015600492421,
         0.082502225484340941, 0.12007019910960293,0.16415283300752470,
         0.21442376986779355, 0.27051082840644336, 0.33199876341447887,
         0.39843234186401943, 0.46931971407375483, 0.54413605556657973,
         0.62232745288031077, 0.70331500465597174, 0.78649910768313447,
         0.87126389619061517,0.95698180152629142
        };
      double W[18]=
        {
         0.0055657196642445571,
         0.012915947284065419,0.020181515297735382,0.027298621498568734,
         0.034213810770299537,0.040875750923643261,0.047235083490265582,
         0.053244713977759692,0.058860144245324798,0.064039797355015485,
         0.068745323835736408,0.072941885005653087,0.076598410645870640,
         0.079687828912071670,0.082187266704339706,0.084078218979661945,
         0.085346685739338721,0.085983275670394821
        };
      setCGauleg18(Ngau,Y,W);
     };
   void setCGauleg18(int Ngau,double &Y[],double &W[]) //set-method for CGauleg18
     {
      ngau=Ngau;
      ArrayCopy(y,Y);
      ArrayCopy(w,W);
     };
  };
//+------------------------------------------------------------------+
//|            Incomplete Beta Function class definition             |
//+------------------------------------------------------------------+
class CBeta : public CGauleg18
  {
private:
   int               Switch;                     //when to use the quadrature method
   double            Eps,Fpmin;
public:
   //+------------------------------------------------------------------+
   //| CBeta class constructor                                          |
   //+------------------------------------------------------------------+
   void CBeta()
     {
      int swi=3000;
      setCBeta(swi,EPS,FPMIN);
     };
   //+------------------------------------------------------------------+
   //| CBeta class set-method                                           |
   //+------------------------------------------------------------------+
   void setCBeta(int swi,double eps,double fpmin)
     {
      Switch=swi;
      Eps=eps;
      Fpmin=fpmin;
     };
   double            betai(const double a,const double b,const double x); //incomplete beta function Ix(a,b)
   double            betacf(const double a,const double b,const double x);//continued fraction for incomplete beta function
   double            betaiapprox(double a,double b,double x); //Incomplete beta by quadrature
   double            invbetai(double p,double a,double b);    //Inverse of incomplete beta function
  };
//+------------------------------------------------------------------+
//|           Incomplete beta function Ix(a,b)                       |
//+------------------------------------------------------------------+
double CBeta::betai(const double a,const double b,const double x)
/*
   Returns incomplete beta function Ix(a,b) for 
   positive a and b, and x between 0 and 1.   
*/
  {
   double bt;
   if(a<=0.0 || b<=0.0) Alert("Bad a or b in routine betai!");
   if(x< 0.0 || x>1.0) Alert("Bad x in routine betai!");
   if(x == 0.0 || x == 1.0) return x;
   if(a>Switch && b>Switch) return betaiapprox(a,b,x);
   bt=exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));
   if(x<(a+1.0)/(a+b+2.0)) return bt*betacf(a,b,x)/a;
   else return 1.0-bt*betacf(b,a,1.0-x)/b;
  }
//+------------------------------------------------------------------+
//|      Continued fraction for incomplete beta function             |
//+------------------------------------------------------------------+
double CBeta::betacf(const double a,const double b,const double x)
/*
   Evaluates continued fraction for incomplete beta function 
   by modified Lentz’s method. User should not call directly.
*/
  {
   int m,m2;
   double aa,c,d,del,h,qab,qam,qap;
   qab=a+b;
   qap=a+1.0;
   qam=a-1.0;
   c=1.0;
   d=1.0-qab*x/qap;
   if(fabs(d)<FPMIN) d=FPMIN;
   d=1.0/d;
   h=d;
   for(m=1;m<10000;m++)
     {
      m2=2*m;
      aa=m*(b-m)*x/((qam+m2)*(a+m2));
      d=1.0+aa*d;
      if(fabs(d)<FPMIN) d=FPMIN;
      c=1.0+aa/c;
      if(fabs(c)<FPMIN) c=FPMIN;
      d=1.0/d;
      h *= d*c;
      aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
      d=1.0+aa*d;
      if(fabs(d)<FPMIN) d=FPMIN;
      c=1.0+aa/c;
      if(fabs(c)<FPMIN) c=FPMIN;
      d=1.0/d;
      del=d*c;
      h *= del;
      if(fabs(del-1.0)<=EPS) break;
     }
   return h;
  }
//+------------------------------------------------------------------+
//|              Incomplete beta by quadrature                       |
//+------------------------------------------------------------------+
double CBeta::betaiapprox(double a,double b,double x)
/*
   Incomplete beta by quadrature. Returns Ix(a,b).
   User should not call directly.
*/
  {
   int j;
   double xu,t,sum,ans;
   double a1=a-1.0,b1=b-1.0,mu=a/(a+b);
   double lnmu=log(mu),lnmuc=log(1.-mu);
   t=sqrt(a*b/(pow(a+b,2)*(a+b+1.0)));
   if(x>a/(a+b))
     {
      if(x>=1.0) return 1.0;
      xu=fmin(1.,fmax(mu+10.*t,x+5.0*t));
        } else {
      if(x<=0.0) return 0.0;
      xu=fmax(0.,fmin(mu-10.*t,x-5.0*t));
     }
   sum=0;
   for(j=0;j<18;j++)
     {
      t=x+(xu-x)*y[j];
      sum+=w[j]*exp(a1*(log(t)-lnmu)+b1*(log(1-t)-lnmuc));
     }
   ans=sum*(xu-x)*exp(a1*lnmu-gammln(a)+b1*lnmuc-gammln(b)+gammln(a+b));
   return ans>0.0? 1.0-ans : -ans;
  }
//+------------------------------------------------------------------+
//|          Inverse of incomplete beta function                     |
//+------------------------------------------------------------------+
double CBeta::invbetai(double p,double a,double b)
/*
  Inverse of incomplete beta function. Returns x 
  such that Ix(a,b)=p for argument p between 0 and 1.
*/
  {
   const double Eps_local=1.e-8;
   double pp,t,u,err,x,al,h,w_local,afac,a1=a-1.,b1=b-1.;
   int j;
   if(p<=0.) return 0.;
   else if(p>=1.) return 1.;
   else if(a>=1. && b>=1.)
     {
      pp=(p<0.5)? p : 1.-p;
      t = sqrt(-2.*log(pp));
      x = (2.30753+t*0.27061)/(1.+t*(0.99229+t*0.04481)) - t;
      if(p<0.5) x=-x;
      al=(pow(x,2)-3.)/6.;
      h = 2./(1./(2.*a-1.)+1./(2.*b-1.));
      w_local = (x*sqrt(al+h)/h)-(1./(2.*b-1)-1./(2.*a-1.))*(al+5./6.-2./(3.*h));
      x = a/(a+b*exp(2.*w_local));
        } else {
      double lna=log(a/(a+b)),lnb=log(b/(a+b));
      t = exp(a*lna)/a;
      u = exp(b*lnb)/b;
      w_local = t + u;
      if(p<t/w_local) x=pow(a*w_local*p,1./a);
      else x=1.-pow(b*w_local*(1.-p),1./b);
     }
   afac = -gammln(a)-gammln(b)+gammln(a+b);
   for(j=0;j<10;j++)
     {
      if(x==0. || x==1.) return x;
      err = betai(a,b,x)-p;
      t = exp(a1*log(x)+b1*log(1.-x) + afac);
      u = err/t;
      x-=(t=u/(1.-0.5*fmin(1.,u*(a1/x-b1/(1.-x)))));
      if(x <= 0.) x = 0.5*(x + t);
      if(x >= 1.) x = 0.5*(x + t + 1.);
      if(fabs(t)<Eps_local*x && j>0) break;
     }
   return x;
  }
//+------------------------------------------------------------------+
//|            Incomplete Gamma Function class definition            |
//+------------------------------------------------------------------+
class CGamma : public CGauleg18
  {
private:
   int               ASWITCH;
   double            Eps,
   Fpmin,
   gln;
public:
   //+------------------------------------------------------------------+
   //| CGamma class constructor                                         |
   //+------------------------------------------------------------------+
   void CGamma()
     {
      int aswi=100;
      setCGamma(aswi,EPS,FPMIN);
     };
   void setCGamma(int aswi,double eps,double fpmin) //set-method for CGamma
     {
      ASWITCH=aswi;
      Eps=eps;
      Fpmin=fpmin;
     };
   double            gammp(const double a,const double x);
   double            gammq(const double a,const double x);
   void              gser(double &gamser,double a,double x,double &gln);
   double            gcf(const double a,const double x);
   double            gammpapprox(double a,double x,int psig);
   double            invgammp(double p,double a);
  };
//+------------------------------------------------------------------+
//|          Incomplete gamma function P(a,x)                        |
//+------------------------------------------------------------------+
double CGamma::gammp(const double a,const double x)
//Returns the incomplete gamma function P(a,x)
  {
   double gamser,gln_local;//gammcf
   if(x<0.0 || a<=0.0) Alert("bad args in gammp!");
   if(x==0.0) return 0.0;
   else if((int)a>=ASWITCH) return gammpapprox(a,x,1); //Quadrature
   else if(x<a+1.0)
     {
      gser(gamser,a,x,gln_local);     //Use the series representation
      return gamser;
     }
   else return 1.0-gcf(a,x);          //Use the continued fraction representation
  };
//+------------------------------------------------------------------+
//|          Incomplete gamma function Q(a,x)                        |
//+------------------------------------------------------------------+
double CGamma::gammq(const double a,const double x)
//Returns the incomplete gamma function Q(a,x) ? 1 - P(a,x)
  {
   double gamser,gln_local;//gammcf
   if(x<0.0 || a<=0.0) Alert("bad args in gammq!");
   if(x==0.0) return 1.0;
   else if((int)a>=ASWITCH) return gammpapprox(a,x,0); //Quadrature
   else if(x<a+1.0)
     {
      gser(gamser,a,x,gln_local);     //Use the series representation
      return 1.0-gamser;
     }
   else return gcf(a,x);              //Use the continued fraction representation
  };
//+------------------------------------------------------------------+
//|             Incomplete gamma function P(a,x)                     |
//+------------------------------------------------------------------+
void CGamma::gser(double &gamser,double a,double x,double &gln_local)
/* 
   Returns the incomplete gamma function P(a,x) evaluated
   by its series representation as gamser. Also sets ln Ã(a) as gln.
   User should not call directly.
*/
  {
   int n;
   double sum,del,ap;
   gln_local=gammln(a);
   if(x<=0.0)
     {
      if(x<0.0) nrerror("x less than 0 in routine gser");
      gamser=0.0;
      return;
     }
   else
     {
      ap=a;
      del=sum=1.0/a;
      for(n=1;n<=ITMAX;n++)
        {
         ++ap;
         del *= x/ap;
         sum += del;
         if(fabs(del)<fabs(sum)*EPS)
           {
            gamser=sum*exp(-x+a*log(x)-(gln_local));
            return;
           }
        }
      nrerror("a too large, ITMAX too small in routine gser");
      return;
     }
  };
//+------------------------------------------------------------------+
//|             Incomplete gamma function Q(a,x)                     |
//+------------------------------------------------------------------+
double CGamma::gcf(const double a,const double x)
/*
  Returns the incomplete gamma function Q(a,x) 
  evaluated by its continued fraction representation.
  Also sets ln Ã(a) as gln. User should not call directly.
*/
  {
   int i;
   double an,b,c,d,del,h;
   gln=gammln(a);
   b=x+1.0-a;                     //Set up for evaluating continued fraction
   c=1.0/FPMIN;                   //by modified Lentz’s method
   d=1.0/b;                       //with b0=0
   h=d;
   for(i=1;;i++) //Iterate to convergence
     {
      an = -i*(i-a);
      b += 2.0;
      d=an*d+b;
      if(fabs(d)<FPMIN) d=FPMIN;
      c=b+an/c;
      if(fabs(c)<FPMIN) c=FPMIN;
      d=1.0/d;
      del=d*c;
      h *= del;
      if(fabs(del-1.0)<=EPS) break;
     }
   return exp(-x+a*log(x)-gln)*h;    //Put factors in front
  };
//+------------------------------------------------------------------+
//|             Incomplete gamma by quadrature                       |
//+------------------------------------------------------------------+
double CGamma::gammpapprox(double a,double x,int psig)
/*
  Incomplete gamma by quadrature. Returns P(a,x) or Q(a,x),
  when psig is 1 or 0,respectively. User should not call directly.
*/
  {
   int j;
   double xu,t,sum,ans;
   double a1=a-1.0,lna1=log(a1),sqrta1=sqrt(a1);
   gln=gammln(a);
//Set how far to integrate into the tail:
   if(x>a1) xu=fmax(a1+11.5*sqrta1,x+6.0*sqrta1);
   else xu=fmax(0.,fmin(a1-7.5*sqrta1,x-5.0*sqrta1));
   sum=0;
   for(j=0;j<ngau;j++) //Gauss-Legendre
     {
      t=x+(xu-x)*y[j];
      sum+=w[j]*exp(-(t-a1)+a1*(log(t)-lna1));
     }
   ans=sum*(xu-x)*exp(a1*(lna1-1.)-gln);
   return(psig?(ans>0.0? 1.0-ans:-ans):(ans>=0.0? ans:1.0+ans));
  };
//+------------------------------------------------------------------+
//|         Inverse of incomplete gamma function                     |
//+------------------------------------------------------------------+
double CGamma::invgammp(double p,double a)
//Returns x such that P(a,x)=p for an argument p between 0 and 1.   
  {
   int j;
   double x,err,t,u,pp,lna1=0.,afac=0.,a1=a-1;
   const double Eps_local=1.e-8;    //Accuracy is the square of EPS
   gln=gammln(a);
   if(a <= 0.) Alert("a must be pos in invgammap!");
   if(p >= 1.) return fmax(100.,a + 100.*sqrt(a));
   if(p <= 0.) return 0.0;
   if(a>1.) //Initial guess based on reference [1]
     {
      lna1=log(a1);
      afac= exp(a1*(lna1-1.)-gln);
      pp=(p<0.5)? p : 1.-p;
      t = sqrt(-2.*log(pp));
      x = (2.30753+t*0.27061)/(1.+t*(0.99229+t*0.04481)) - t;
      if(p<0.5) x=-x;
      x=fmax(1.e-3,a*pow(1.-1./(9.*a)-x/(3.*sqrt(a)),3));
        } else {              //Initial guess based on equations (6.2.8)
      t=1.0-a*(0.253+a*0.12); //and (6.2.9)
      if(p<t) x=pow(p/t,1./a);
      else x=1.-log(1.-(p-t)/(1.-t));
     }
   for(j=0;j<12;j++)
     {
      if(x<=0.0) return 0.0;  //x too small to compute accurately
      err=gammp(a,x)-p;
      if(a>1.) t=afac*exp(-(x-a1)+a1*(log(x)-lna1));
      else t=exp(-x+a1*log(x)-gln);
      u=err/t;
      x-=(t=u/(1.-0.5*fmin(1.,u*((a-1.)/x-1)))); //Halley’s method
      if(x<=0.) x=0.5*(x+t);  //Halve old value if x tries to go negative
      if(fabs(t)<Eps_local*x) break;
     }
   return x;
  }
//+------------------------------------------------------------------+
//|                Student-t Distribution class definition           |
//+------------------------------------------------------------------+
class CStudenttdist : CBeta // CBeta class inheritance
  {
public:
   int               nu;    //shape parameter(?)
   double            mu,    //location parameter (?)
   sig,                     //scale parameter (?)
   np,                      // 1/2*(?+1)
   fac;                     // Ã(1/2*(?+1))-Ã(1/2*?)
   //+------------------------------------------------------------------+
   //| CStudenttdist class constructor                                  |
   //+------------------------------------------------------------------+
   void  CStudenttdist()
     {
      int Nu=1;double Mu=0.0,Sig=1.0; //parameters ?, ? and ? by default
      setCStudenttdist(Nu,Mu,Sig);
     }
   void setCStudenttdist(int Nu,double Mu,double Sig)
     {
      nu=Nu;
      mu=Mu;
      sig=Sig;
      if(sig<=0. || nu<=0.) Alert("bad sig,nu in Student-t Distribution!");
      np=0.5*(nu+1.);
      fac=gammln(np)-gammln(0.5*nu);
     }
   //+------------------------------------------------------------------+
   //| Probability density function                                     |
   //| Probability density function (pdf)                               |
   //+------------------------------------------------------------------+
   double pdf(double x)
     {
      return exp(-np*log(1.+pow((x-mu)/sig,2.)/nu)+fac)/(sqrt(M_PI*nu)*sig);
     }
   //+------------------------------------------------------------------+
   //| Cumulative distribution function                                 |
   //| Cumulative distribution function (cdf)                           |
   //+------------------------------------------------------------------+
   double cdf(double t)
     {
      double p=0.5*betai(0.5*nu,0.5,nu/(nu+pow((t-mu)/sig,2)));
      if(t>=mu) return 1.-p;
      else return p;
     }
   //+------------------------------------------------------------------+
   //| Inverse cumulative distribution function (quantile function)     |
   //| Inverse cumulative distribution function (invcdf)                |
   //+------------------------------------------------------------------+
   double invcdf(double p)
     {
      if(p<=0. || p>=1.) Alert("bad p in Student-t Distribution!");
      double x=invbetai(2.*fmin(p,1.-p),0.5*nu,0.5);
      x=sig*sqrt(nu*(1.-x)/x);
      return(p>=0.5? mu+x : mu-x);
     }
   //+------------------------------------------------------------------+
   //| Reliability (survival) function                                  |
   //| Survival function (sf)                                           |
   //+------------------------------------------------------------------+
   double sf(double x)
     {
      return 1-cdf(x);
     }
   //+------------------------------------------------------------------+
   //| The two-tailed cumulative distribution function                  |
   //| The two-tailed cumulative distribution function (aa) A(t|?)      |
   //+------------------------------------------------------------------+
   double aa(double t)
     {
      if(t < 0.) Alert("bad t in Student-t Distribution!");
      return 1.-betai(0.5*nu,0.5,nu/(nu+pow(t,2.)));
     }
   //+------------------------------------------------------------------+
   //| The inverse two-tailed cumulative distribution function          |
   //| The inverse two-tailed cumulative distribution function (invaa)  |
   //| p=A(t|?)                                                         |
   //+------------------------------------------------------------------+
   double invaa(double p)
     {
      if(!(p>=0. && p<1.)) Alert("bad p in Student-t Distribution!");
      double x=invbetai(1.-p,0.5*nu,0.5);
      return sqrt(nu*(1.-x)/x);
     }
  };
//+------------------------------------------------------------------+
double gammln(const double xx)
//Returns the value ln[Ã(xx)] for xx > 0.
  {
   int j;
   double x,tmp,tmp1,y,ser;
   static const double cof[14]=
     {
      57.1562356658629235,-59.5979603554754912,
      14.1360979747417471,-0.491913816097620199,.339946499848118887e-4,
      .465236289270485756e-4,-.983744753048795646e-4,.158088703224912494e-3,
      -.210264441724104883e-3,.217439618115212643e-3,-.164318106536763890e-3,
      .844182239838527433e-4,-.261908384015814087e-4,.368991826595316234e-5
     };
   if(xx<=0)
      Alert("bad arg in gammln!");
   y=x=xx;
   tmp1= x+5.24218750000000000;
   tmp = (x+0.5)*log(tmp1)-tmp1;
   ser = 0.999999999999997092;
   for(j=0;j<14;j++)
      ser+=cof[j]/++y;
   return tmp+log(2.5066282746310005*ser/x);
  }
//+------------------------------------------------------------------+
//|                    Error print function                          |
//+------------------------------------------------------------------+
void nrerror(string error_text)
//Registers errors
  {
   Print("Numerical Recipes run-time error...\n");
   Print("%s\n",error_text);
   Print("...now exiting to system...\n");
   return;
  }
//+------------------------------------------------------------------+
//|                Logistic Distribution class definition            |
//+------------------------------------------------------------------+
class CLogisticdist
  {
public:
   double            alph,//location parameter (?)
   bet;                   //scale parameter (?)
   //+------------------------------------------------------------------+
   //| CLogisticdist class constructor                                  |
   //+------------------------------------------------------------------+
   void  CLogisticdist()
     {
      alph=0.0;bet=1.0; //parameters ? and ? by default
      if(bet<=0.) Alert("bad bet in Logistic Distribution!");
     }
   //+------------------------------------------------------------------+
   //| Probability density function                                     |
   //| Probability density function (pdf)                               |
   //+------------------------------------------------------------------+
   double pdf(double x)
     {
      return exp(-(x-alph)/bet)/(bet*pow(1.+exp(-(x-alph)/bet),2));
     }
   //+------------------------------------------------------------------+
   //| Cumulative distribution function                                 |
   //| Cumulative distribution function (cdf)                           |
   //+------------------------------------------------------------------+
   double cdf(double x)
     {
      double et=exp(-1.*fabs(1.81379936423421785*(x-alph)/bet));
      if(x>=alph) return 1./(1.+et);
      else return et/(1.+et);
     }
   //+------------------------------------------------------------------+
   //| Inverse cumulative distribution function (quantile function)     |
   //| Inverse cumulative distribution function (invcdf)                |
   //+------------------------------------------------------------------+
   double invcdf(double p)
     {
      if(p<=0. || p>=1.) Alert("bad p in Logistic Distribution!");
      return alph+0.551328895421792049*bet*log(p/(1.-p));
     }
   //+------------------------------------------------------------------+
   //| Reliability (survival) function                                  |
   //| Survival function (sf)                                           |
   //+------------------------------------------------------------------+
   double sf(double x)
     {
      return 1-cdf(x);
     }
  };
//+------------------------------------------------------------------+
//|                Exponential Distribution class definition         |
//+------------------------------------------------------------------+
class CExpondist
  {
public:
   double            lambda;   //scale parameter (?)
   //+------------------------------------------------------------------+
   //| CExpondist class constructor                                     |
   //+------------------------------------------------------------------+
   void  CExpondist()
     {
      lambda=1.0;              //parameter ? by default
      if(lambda<=0.) Alert("bad lambda in Exponential Distribution!");
     }
   //+------------------------------------------------------------------+
   //| Probability density function                                     |
   //| Probability density function (pdf)                               |
   //+------------------------------------------------------------------+
   double pdf(double x)
     {
      if(x<0.) Alert("bad x in Exponential Distribution!");
      return lambda*exp(-lambda*x);
     }
   //+------------------------------------------------------------------+
   //| Cumulative distribution function                                 |
   //| Cumulative distribution function (cdf)                           |
   //+------------------------------------------------------------------+
   double cdf(double x)
     {
      if(x < 0.) Alert("bad x in Exponential Distribution!");
      return 1.-exp(-lambda*x);
     }
   //+------------------------------------------------------------------+
   //| Inverse cumulative distribution function (quantile function)     |
   //| Inverse cumulative distribution function (invcdf)                |
   //+------------------------------------------------------------------+
   double invcdf(double p)
     {
      if(p<0. || p>=1.) Alert("bad p in Exponential Distribution!");
      return -log(1.-p)/lambda;
     }
   //+------------------------------------------------------------------+
   //| Reliability (survival) function                                  |
   //| Survival function (sf)                                           |
   //+------------------------------------------------------------------+
   double sf(double x)
     {
      return 1-cdf(x);
     }
  };
//+------------------------------------------------------------------+
//|                Gamma Distribution class definition               |
//+------------------------------------------------------------------+
class CGammadist : CGamma // CGamma class inheritance
  {
public:
   double            alph,//continuous shape parameter (?>0)
   bet,                   //continuous scale parameter  (?>0)
   fac;                   //factor
   //+------------------------------------------------------------------+
   //| CGammaldist class constructor                                    |
   //+------------------------------------------------------------------+
   void  CGammadist()
     {
      setCGammadist();
     }
   void setCGammadist(double Alph=1.0,double Bet=1.0)//parameters ? and ? by default
     {
      alph=Alph; bet=Bet;
      if(alph<=0. || bet<=0.) Alert("bad alph,bet in Gamma Distribution!");
      fac=alph*log(bet)-gammln(alph);
     }
   //+------------------------------------------------------------------+
   //| Probability density function                                     |
   //| Probability density function (pdf)                               |
   //+------------------------------------------------------------------+
   double pdf(double x)
     {
      if(x<=0.) Alert("bad x in Gamma Distribution!");
      return exp(-bet*x+(alph-1.)*log(x)+fac);
     }
   //+------------------------------------------------------------------+
   //| Cumulative distribution function                                 |
   //| Cumulative distribution function (cdf)                           |
   //+------------------------------------------------------------------+
   double cdf(double x)
     {
      if(x<0.) Alert("bad x in Gamma Distribution!");
      return gammp(alph,bet*x);
     }
   //+------------------------------------------------------------------+
   //| Inverse cumulative distribution function (quantile function)     |
   //| Inverse cumulative distribution function (invcdf)                |
   //+------------------------------------------------------------------+
   double invcdf(double p)
     {
      if(p<0. || p>=1.) Alert("bad p in Gamma Distribution!");
      return invgammp(p,alph)/bet;
     }
   //+------------------------------------------------------------------+
   //| Reliability (survival) function                                  |
   //| Survival function (sf)                                           |
   //+------------------------------------------------------------------+
   double sf(double x)
     {
      return 1-cdf(x);
     }
  };
//+------------------------------------------------------------------+
//|                Beta Distribution class definition                |
//+------------------------------------------------------------------+
class CBetadist : CBeta // CBeta class inheritance
  {
public:
   double            alph,//continuous shape parameter (?>0)
   bet,                   //continuous shape parameter (?>0)
   fac;                   //factor
   //+------------------------------------------------------------------+
   //| CBetadist class constructor                                      |
   //+------------------------------------------------------------------+
   void  CBetadist()
     {
      setCBetadist();
     }
   void setCBetadist(double Alph=0.5,double Bet=0.5)//parameters ? and ? by default
     {
      alph=Alph; bet=Bet;
      if(alph<=0. || bet<=0.) Alert("bad alph,bet in Beta Distribution!");
      fac=gammln(alph+bet)-gammln(alph)-gammln(bet);
     }
   //+------------------------------------------------------------------+
   //| Probability density function                                     |
   //| Probability density function (pdf)                               |
   //+------------------------------------------------------------------+
   double pdf(double x)
     {
      if(x<=0. || x>=1.) Alert("bad x in Beta Distribution!");
      return exp((alph-1.)*log(x)+(bet-1.)*log(1.-x)+fac);
     }
   //+------------------------------------------------------------------+
   //| Cumulative distribution function                                 |
   //| Cumulative distribution function (cdf)                           |
   //+------------------------------------------------------------------+
   double cdf(double x)
     {
      if(x<0. || x>1.) Alert("bad x in Beta Distribution");
      return betai(alph,bet,x);
     }
   //+------------------------------------------------------------------+
   //| Inverse cumulative distribution function (quantile function)     |
   //| Inverse cumulative distribution function (invcdf)                |
   //+------------------------------------------------------------------+
   double invcdf(double p)
     {
      if(p<0. || p>1.) Alert("bad p in Beta Distribution!");
      return invbetai(p,alph,bet);
     }
   //+------------------------------------------------------------------+
   //| Reliability (survival) function                                  |
   //| Survival function (sf)                                           |
   //+------------------------------------------------------------------+
   double sf(double x)
     {
      return 1-cdf(x);
     }
  };
//+------------------------------------------------------------------+
//|                 Laplace Distribution class definition            |
//+------------------------------------------------------------------+
class CLaplacedist
  {
public:
   double            alph;   //location parameter (?)
   double            bet;    //scale parameter (?)
   //+------------------------------------------------------------------+
   //| CLaplacedist class constructor                                   |
   //+------------------------------------------------------------------+
   void  CLaplacedist()
     {
      alph=.0;               //parameter ? by default
      bet=1.;                //parameter ? by default
      if(bet<=0.) Alert("bad bet in Laplace Distribution!");
     }
   //+------------------------------------------------------------------+
   //| Probability density function                                     |
   //| Probability density function (pdf)                               |
   //+------------------------------------------------------------------+
   double pdf(double x)
     {
      return exp(-fabs((x-alph)/bet))/2*bet;
     }
   //+------------------------------------------------------------------+
   //| Cumulative distribution function                                 |
   //| Cumulative distribution function (cdf)                           |
   //+------------------------------------------------------------------+
   double cdf(double x)
     {
      double temp;
      if(x<0)
         temp=0.5*exp(-fabs((x-alph)/bet));
      else
         temp=1.-0.5*exp(-fabs((x-alph)/bet));
      return temp;
     }
   //+------------------------------------------------------------------+
   //| Inverse cumulative distribution function (quantile function)     |
   //| Inverse cumulative distribution function (invcdf)                |
   //+------------------------------------------------------------------+
   double invcdf(double p)
     {
      double temp;
      if(p<0. || p>=1.) Alert("bad p in Laplace Distribution!");
      if(p<0.5)
         temp=bet*log(2*p)+alph;
      else
         temp=-1.*(bet*log(2*(1.-p))+alph);
      return temp;
     }
   //+------------------------------------------------------------------+
   //| Reliability (survival) function                                  |
   //| Survival function (sf)                                           |
   //+------------------------------------------------------------------+
   double sf(double x)
     {
      return 1-cdf(x);
     }
  };
//+------------------------------------------------------------------+
//|               Binomial Distribution class definition             |
//+------------------------------------------------------------------+
class CBinomialdist : CBeta // CBeta class inheritance
  {
public:
   int               n;      //number of trials
   double            pe,     //success probability
   fac;                      //factor
   //+------------------------------------------------------------------+
   //| CBinomialdist class constructor                                  |
   //+------------------------------------------------------------------+
   void              CBinomialdist()
     {
      setCBinomialdist();
     }
   void setCBinomialdist(int N=100,double Pe=0.5)//default parameters n and pe
     {
      n=N; pe=Pe;
      if(n<=0 || pe<=0. || pe>=1.) Alert("bad args in Binomial Distribution!");
      fac=gammln(n+1.);
     }
   //+------------------------------------------------------------------+
   //| Probability density function                                     |
   //| Probability density function (pdf)                               |
   //+------------------------------------------------------------------+
   double pdf(int k)
     {
      if(k<0) Alert("bad k in Binomial Distribution!");
      if(k>n) return 0.;
      return exp(k*log(pe)+(n-k)*log(1.-pe)+fac-gammln(k+1.)-gammln(n-k+1.));
     }
   //+------------------------------------------------------------------+
   //| Cumulative distribution function                                 |
   //| Cumulative distribution function (cdf)                           |
   //+------------------------------------------------------------------+
   double cdf(int k)
     {
      if(k<0) Alert("bad k in Binomial Distribution!");
      if(k==0) return 0.;
      if(k>n) return 1.;
      return 1.-betai((double)k,n-k+1.,pe);
     }
   //+------------------------------------------------------------------+
   //| Inverse cumulative distribution function (quantile function)     |
   //| Inverse cumulative distribution function (invcdf)                |
   //+------------------------------------------------------------------+
   int invcdf(double p)
     {
      int k,kl,ku,inc=1;
      if(p<=0. || p>=1.) Alert("bad p in Binomial Distribution!");
      k=(int)fmax(0,fmin(n,(int)(n*pe)));
      if(p<cdf(k))
        {
         do
           {
            k=(int)fmax(k-inc,0);
            inc*=2;
           }
         while(p<cdf(k));
         kl=k; ku=k+inc/2;
           } else {
         do
           {
            k=fmin(k+inc,n+1);
            inc*=2;
           }
         while(p>cdf(k));
         ku=k; kl=k-inc/2;
        }
      while(ku-kl>1)
        {
         k=(kl+ku)/2;
         if(p<cdf(k)) ku=k;
         else kl=k;
        }
      return kl;
     }
   //+------------------------------------------------------------------+
   //| Reliability (survival) function                                  |
   //| Survival function (sf)                                           |
   //+------------------------------------------------------------------+
   double sf(int k)
     {
      return 1.-cdf(k);
     }
  };
//+------------------------------------------------------------------+
//|               Poisson Distribution class definition              |
//+------------------------------------------------------------------+
class CPoissondist : CGamma // CGamma class inheritance
  {
public:
   double            lambda;   //location parameter (?)
   //+------------------------------------------------------------------+
   //| CPoissondist class constructor                                   |
   //+------------------------------------------------------------------+
   void CPoissondist()
     {
      lambda=15.;
      if(lambda<=0.) Alert("bad lambda in Poisson Distribution!");
     }
   //+------------------------------------------------------------------+
   //| Probability density function                                     |
   //| Probability density function (pdf)                               |
   //+------------------------------------------------------------------+
   double pdf(int n)
     {
      if(n<0) Alert("bad n in Poisson Distribution!");
      return exp(-lambda+n*log(lambda)-gammln(n+1.));
     }
   //+------------------------------------------------------------------+
   //| Cumulative distribution function                                 |
   //| Cumulative distribution function (cdf)                           |
   //+------------------------------------------------------------------+
   double cdf(int n)
     {
      if(n<0) Alert("bad n in Poisson Distribution!");
      if(n==0) return 0.;
      return gammq((double)n,lambda);
     }
   //+------------------------------------------------------------------+
   //| Inverse cumulative distribution function (quantile function)     |
   //| Inverse cumulative distribution function (invcdf)                |
   //+------------------------------------------------------------------+
   int invcdf(double p)
     {
      int n,nl,nu,inc=1;
      if(p<=0. || p>=1.) Alert("bad p in Poisson Distribution!");
      if(p<exp(-lambda)) return 0;
      n=(int)fmax(sqrt(lambda),5.);
      if(p<cdf(n))
        {
         do
           {
            n=(int)fmax(n-inc,0);
            inc*=2;
           }
         while(p<cdf(n));
         nl=n; nu=n+inc/2;
           } else {
         do
           {
            n+=inc;
            inc*=2;
           }
         while(p>cdf(n));
         nu=n; nl=n-inc/2;
        }
      while(nu-nl>1)
        {
         n=(nl+nu)/2;
         if(p<cdf(n)) nu=n;
         else nl=n;
        }
      return nl;
     }
   //+------------------------------------------------------------------+
   //| Reliability (survival) function                                  |
   //| Survival function (sf)                                           |
   //+------------------------------------------------------------------+
   double sf(int n)
     {
      return 1.-cdf(n);
     }
  };
//+=====================================================================+
//[EOF]