/* file: probability-distribution.h (version 0.0)
   last modification: august 17, 2010
   use: academic purpose, limited accuracy
   author: Anthony Phan
*/

#include <math.h>
#define infinity 1e12 /* NaN;*/
static double eps = 1e-12;
static double accuracy = 1e-12;
static double normlimit = 10.0;
static double PI = 3.141592653589793;

static double min(double x, double y){if (x < y) return x; else return y;}

static double max(double x, double y){if (x < y) return y; else return x;}

double lnGamma(double x){
  return log(x+4.5)*(x-.5)-x-4.5
    +log(2.50662827465*(1+76.18009173/x
			-86.50532033/(x+1)
			+24.01409822/(x+2)
			-1.231739516/(x+3)
			+0.00120858003/(x+4)
			-0.00000536382/(x+5)));
}

double Gamma(double x){
  return exp(lnGamma(x));
}

double Beta(double x, double y){
  return exp(lnGamma(x)+lnGamma(y)-lnGamma(x+y));
}

/* Gamma distributions */

double gammapdf(double x, double a, double lambda){
  if (x <= 0) return 0;
  else return lambda*exp((a-1)*log(lambda*x)-lambda*x-lnGamma(a));
}
  
double gammacdf(double x, double a, double lambda){
  float x_; x_ = lambda*x;
  if (x_ <= 0 || a <= 0) return 0;
  else {
    if (x_ <= a+1){
      double sum, del; int n;
      sum = 1/a; del = sum;
      for (n = 1; n <= 100; n++){
	del = del*x_/(a+n); sum = sum+del;
	if (fabs(del) < fabs(sum)*accuracy) break;
      }
      return sum*exp(-x_+a*log(x_)-lnGamma(a));
    }
    else { 
      double g, gold, aa, aaa, bb, bbb, fac; int n;
      gold = 0.0; aa = 1; aaa = x_; bb = 0.0; bbb = 1; fac = 1;
      for (n = 1; n <= 100; n++){
	aa = (aaa+aa*(n-a))*fac;
	bb = (bbb+bb*(n-a))*fac;
	aaa = x_*aa+n*fac*aaa;
	bbb = x_*bb+n*fac*bbb;
	if (aaa != 0){fac = 1/aaa; g = bbb*fac;}
	if (fabs((g-gold)/g) < accuracy) break;
	gold = g;
      }
      return 1-exp(-x_+a*log(x_)-lnGamma(a))*g;
    }
  }
}

double gammaicdf(double p, double a, double lambda){
  if (p < accuracy) return 0;
  else if (p > 1-accuracy) return infinity;
    else {double lb, ub; lb = 0.0; ub = 1.0;
      while (gammacdf(ub, a, lambda) <= p){lb = ub; ub = ub*2;}
      while (ub-lb > accuracy)
	if (gammacdf((lb+ub)/2, a, lambda) < p)
	  lb = (lb+ub)/2; else ub = (lb+ub)/2;
      return (lb+ub)/2;
    }
}

/* N(0, 1) distribution */

double normalpdf(double x){return .3989423*exp(-x*x/2);}

double normalcdf(double x){
  if (fabs(x) < normlimit){
    if (x == 0) return 0.5;
    else {double u, p; int n;
      u = x/sqrt(2*PI); p = 0.5+u; n = 0;
      while (n <= x*x || fabs(u/(2*n+1)) >= accuracy){
	n++;
	u = -u*x*x/(2*n);
	p += u/(2*n+1);
      }
      return max(min(p, 1), 0);
    }
  }
  else {
    if (x < 0) return 0; else return 1;
  }
}

double normalcdf_(double x){
  double p; p = 0.5*gammacdf(0.5*x*x, 0.5, 1)+0.5;
  if (x < 0) p = -p; return p;
}

double normalicdf(double p){
  if (p > 1-accuracy) return infinity;
  else if (p < accuracy) return -infinity;
  else if (p == 0.5) return 0;
  else {double x, lb, ub;
    lb = -normlimit; ub = normlimit;
    while (ub-lb > accuracy){
      if (normalcdf((lb+ub)/2) < p) lb = (lb+ub)/2; else ub = (lb+ub)/2;
    }
    return (lb+ub)/2;
  }
}

double normalicdf_(double p){
  if (p > 1-accuracy) return infinity;
  else if (p < accuracy) return -infinity;
  else if (p == 0.5) return 0;
  else {double x, lb, ub;
    lb = -normlimit; ub = normlimit;
    while (ub-lb > accuracy){
      if (normalcdf_((lb+ub)/2) < p) lb = (lb+ub)/2; else ub = (lb+ub)/2;
    }
    return (lb+ub)/2;
  }
}

/* N(m, sigma) normal or Gaussian distributions
double gaussianpdf(double x, double m, double sigma){
  return normalpdf((x-m)/sigma)/sigma;}
double gaussiancdf(double x, double m, double sigma){
  return normalcdf((x-m)/sigma);}
double gaussianicdf(double p, double m, double sigma){
  return sigma*normalicdf(p)+m;}
*/

/* $\chi^2$ distributions (Pearson) */

double chisquarepdf(double x, double df){
  if (x > 0) return exp((0.5*df-1)*log(x)-x/2-0.5*df*log(2)-lnGamma(.5*df));
  else return 0;
}

double chisquarecdf(double x, double df){
  if (x <= 0) return 0; else return gammacdf(0.5*x, 0.5*df, 1);
}

double chisquareicdf(double p, double df){
  if (p < accuracy) return 0;
  else if (p > 1-accuracy) return infinity;
  else return 2*gammaicdf(p, 0.5*df, 1);
}

/* Beta distributions */

double betapdf(double x, double a, double b){
  if (x <= 0 || x >= 1) return  0;
  else return pow(x, a-1)*pow(1-x, b-1)/Beta(a, b); 
}

double _betacdf_(double x, double a, double b){
  double am, bm, az, bz, d, ap, bp, app, bpp, aold; int m;
  am = 1; bm = 1; az = 1; bz = 1-((a+b)/(a+1))*x;
  for (m = 1; m <= 100; m++){
    d = (m/(a-1+2*m))*((b-m)/(a+2*m))*x;
    ap = az+d*am;
    bp = bz+d*bm;
    d = -((a+m)/(a+2*m))*((a+b+m)/(a+1+2*m))*x;
    app = ap+d*az;
    bpp = bp+d*bz;
    aold = az;
    am = ap/bpp;
    bm = bp/bpp;
    az = app/bpp;
    bz = 1;
    if (fabs(az-aold) < accuracy*fabs(az)) break;
  }
  return az;
}

double betacdf(double x, double a, double b){
  if (x <= 0) return 0;
  else if (x >= 1) return 1;
  else {double bt;
    bt = exp(lnGamma(a+b)-lnGamma(a)-lnGamma(b)+a*log(x)+b*log(1-x));
    if (x <= (a+1)/(a+b+2)) return bt*_betacdf_(x, a, b)/a;
    else return 1-bt*_betacdf_(1-x, b, a)/b;
  }
}

double betaicdf(double p, double a, double b){
  if (p < accuracy) return 0;
  else if (p > 1-accuracy) return 1;
  else {double lb, ub; lb = 0.0; ub = 1.0;
    while (ub-lb > accuracy){
      if (betacdf((lb+ub)/2, a, b) < p) lb = (lb+ub)/2; else ub = (lb+ub)/2;
    }
    return (lb+ub)/2;
  }
}

/* Student's (T) distributions */

double studentpdf(double x, double df){
  return pow(1+x*x/df, -0.5*(df+1))/sqrt(df)/Beta(0.5*df, 0.5);
}

double studentcdf(double x, double df){
  double p;
  p = 0.5*betacdf(df/(df+x*x), 0.5*df, 0.5);
  if (x >= 0) return 1-p; else return p;
}

double studenticdf(double p, double df){
  double x;
  if (p > 1-accuracy) return infinity;
  else if (p < accuracy) return -infinity;
  else if (p == 0.5) return 0;
  else if (p < 0.5) return -sqrt(df/betaicdf(2*p, 0.5*df, 0.5)-df);
  else  return sqrt(df/betaicdf(2-2*p, 0.5*df, 0.5)-df);
    /* {double lb, ub;
       if (p > 0.5){lb = 0.0; ub = 1.0;
       while (studentcdf(ub, df) < p){lb = ub; ub = 2*ub;}
       }
       else{lb = -1.0; ub = 0.0;
       while (studentcdf(lb, df) > p){ub = lb; lb = 2*lb;}
       }
       while (ub-lb > accuracy){
       if (studentcdf((lb+ub)/2, df) < p) lb =(lb+ub)/2; else ub = (lb+ub)/2;
       }
       return (lb+ub)/2;
       } */
}

/* Fisher's (F) distributions */

double fisherpdf(double x, double ndf, double ddf){
  if (x <= 0) return 0;
  else return exp(0.5*ndf*log(ndf/ddf)+(0.5*ndf-1)*log(x)
	       -0.5*(ndf+ddf)*log(1+(ndf/ddf)*x))/Beta(0.5*ndf, 0.5*ddf);
}

double fishercdf(double x, double ndf, double ddf){
  if (x <= 0) return 0;
  else return 1-betacdf(ddf/(ddf+ndf*x), 0.5*ddf, 0.5*ndf);
}

double fishericdf(double p, double ndf, double ddf){
  if (p > 1-accuracy) return infinity;
  else if (p < accuracy) return 0;
  else return ddf/ndf*(1/betaicdf(1-p, 0.5*ddf, 0.5*ndf)-1);
     /* {
       double lb, ub;
       lb = 0.0; ub = 1.0;
       while (fishercdf(ub, ndf, ddf) < p){lb = ub; ub = 2*ub;}
       while (ub-lb > accuracy){
       if (fishercdf((lb+ub)/2, ndf, ddf) < p) lb = (lb+ub)/2;
       else ub = (lb+ub)/2; 
       }
       return (lb+ub)/2;
       } */
}

/* Poisson's distributions */

double poissonpdf(double x, double lambda){
  if (floor(x) == x && x >= 0){
    if (lambda > 0) return exp(-lambda+x*log(lambda)-lnGamma(x+1));
    else if (x == 0) return 1; else return 0;
  }
  else return 0;
}

double poissoncdf(double x, double lambda){
  if (x < 0) return 0;
  else
    if (lambda > 0){double p; int k;
      p = 1.0;
      for (k = floor(x); k >= 1; k--) p = 1+p/k*lambda;
      return p*exp(-lambda);
    }
    else return 1;
}

double poissonicdf(double p, double lambda){
  if (p <= 0) return 0;
  else if (p >= 1) return infinity;
  else {double x, s, t;
      x = 0.0; s = exp(-lambda); t = s;
      while (s+eps < p){x = x+1; t = t*lambda/x; s = s+t;}
      if (fabs(s-p) < eps) x = x+0.5;
      return x;
    }
}

/* Binomial distributions */

double binomialpdf(double x, int n, double pi){
  if (floor(x) == x && x >= 0 && x <= n)
    if (pi*(1-pi) > 0)
      return exp(lnGamma(n+1)-lnGamma(x+1)-lnGamma(n-x+1)
		 +x*log(pi)+(n-x)*log(1-pi));
    else if ((pi == 0 && x == 0) || (pi == 1 && x == n)) return 1;
      else return 0;
  else return 0;
}

double binomialcdf(double x, int n, double pi){
  if (x < 0) return 0;
  else if (x >= n) return 1;
  else {
      if (pi*(1-pi) > 0){double p, c, l; int k;
	p = 0.0; c = n*log(1-pi); l = log(pi)-log(1-pi);
	for (k = 0; k <= x; k++){
	  p = p+exp(c);
	  c = c+log(n-k)-log(k+1)+l;
	}
	  return p;
      }
	  else if (pi == 0) return 1; else return 0;
  }
}

double binomialicdf(double p, int n, double pi){
  if (p <= 0) return 0;
  else if (p >= 1) return n;
  else {double x, s, t;
    x = 0.0; s = binomialpdf(x, n, pi);
    while (s+eps < p && x < n){
      x = x+1; s = s+binomialpdf(x, n, pi);
    }
    if (fabs(s-p) < eps) return x+0.5; else return x;
  }
}

/* Geometric distributions */

double geometricpdf(double x, double pi){
  if (x >= 1 && floor(x) == x && pi > 0 && pi <= 1)
    return pi*pow(1-pi, x-1);
  else return 0;
}

double geometriccdf(double x, double pi){
    if (x >= 1 && pi > 0 && pi <= 1)
      return 1-pow(1-pi, floor(x));
    else return 0;
  }

double geometricicdf(double p, double pi){
  if (p <= 0) return 0;
  else if (p >= 1) return infinity;
  else
    if (pi > 0 && pi < 1){double x;
      x = log(1-p)/log(1-pi);
      if (floor(x) == x) return floor(x)+0.5; else return ceil(x);
    }
    else if (pi >= 1) return 1; else return infinity;
}
  
/* Negative binomial distribution */

double negativebinomialpdf(double x, int n, double pi){
  if (floor(x) == x && x >= n)
    if (pi*(1-pi) > 0) return exp(lnGamma(x)-lnGamma(n)-lnGamma(x-n+1)
				+n*log(pi)+(x-n)*log(1-pi));
    else if (pi == 1 && x == n) return 1;
    else return 0;
  else return 0; 
}

double negativebinomialcdf(double x, int n, double pi){
  return 1-binomialcdf(n-1, floor(x), pi);
}

double negativebinomialicdf(double p, int n, double pi){
  if (p <= 0) return n;
  else if (p >= 1) return infinity;
  else {double x, s, t;
    x = 1.0*n; s = negativebinomialpdf(x, n, pi);
    while (s+eps < p){
      x = x+1; s = s+negativebinomialpdf(x, n, pi);
    }
    if (fabs(s-p) < eps) return x+0.5; else return x;
  }
}

/* Hypergeometric distributions */

double hypergeometricpdf(double x, int N, int n , double pi){
  /*  if (floor(N*pi) <> N*pi)
    print("Warning: arguments of hypergeometric distribution are wrong.\n");
  */
  if (floor(x) == x && x >= max(0, n-N*(1-pi)) && x <= min(n, N*pi))
    return exp(lnGamma(N*pi+1)+lnGamma(N*(1-pi)+1)-lnGamma(N+1)
	+lnGamma(n+1)-lnGamma(x+1)-lnGamma(n-x+1)
	       +lnGamma(N-n+1)-lnGamma(N*pi-x+1)-lnGamma(N*(1-pi)-n+x+1));
  else return 0;
}
  
double hypergeometriccdf(double x, int N, int n, double pi){
  if (x < max(0, n-N*(1-pi))) return 0;
      else if (x >= min(n, N*pi)) return 1;
      else {double p; int k; p = 0.0;
	for (k = max(0, n-N*(1-pi)); k <= x; k++)
	  {
	    p = p+hypergeometricpdf(k, N, n, pi);
	  }
	return min(p, 1);
      }
}

double hypergeometricicdf(double p, int N, int n, double pi){
  if (p <= 0) return max(0, n-N*(1-pi));
  else if (p >= 1) return min(n, N*pi);
  else {double x, s;
    x = max(0, n-N*(1-pi));
    s = hypergeometricpdf(x, N, n, pi);
    while (s+eps < p && x < min(n, N*pi)){
      x = x+1; s = s+hypergeometricpdf(x, N, n, pi);
    }
    if (fabs(s-p) < eps) return x+0.5; else return x;
  }
}

/* Kolmogorov distributions
   d'apr\`es George Marsaglia and Wai Wan Tsang,
   ``Evaluating Kolmogorov's Distribution''
*/

/* Les matrices sont trait\'ees comme des colonnes
   unidimensionnelles de flottants en double pr\'ecision. */
 
void mMultiply(double *A, double *B, double *C, int m){
  int i, j, k;
  double s;
  for (i = 0; i < m; i++)
    for (j = 0; j < m; j++)
      {s = 0.; /* initialisation d'un flottant */
      for (k = 0; k < m; k++) s += A[i*m+k]*B[k*m+j];
      C[i*m+j] = s;}
}

void mPower(double *A, int eA, double *V, int *eV, int m, int n){
  double *B;
  int eB, i;
  if (n == 1){for (i = 0; i < m*m; i++) V[i] = A[i]; *eV = eA; return;}
  mPower(A, eA, V, eV, m, n/2); /* n/2 est en fait \'egal \`a floor(n/2) */
  B = (double*)malloc((m*m)*sizeof(double));
  mMultiply(V, V, B, m);
  eB = 2*(*eV);
  if (n % 2 == 0)
    {for (i = 0; i < m*m; i++) V[i] = B[i];
    *eV = eB;}
  else
    {mMultiply(A, B, V, m);
    *eV = eA+eB;}
  if (V[(m/2)*m+(m/2)] > 1e140)
    {for (i = 0; i < m*m; i++) V[i] = V[i]*1e-140;
    *eV += 140;}
  free(B);
}

/* to be defined? */
double kolmogorovpdf(double x, int n){return 0;}

double kolmogorovcdf(double x, int n){
  int k, m, i, j, g, eH, eQ;
  double h, s, *H, *Q;
  /* Supprimer les 3 prochaines lignes pour une
     pr\'ecision sup\'erieure \`a 7 d\'ecimales */
  s = x*x*n;
  if (s > 7.24 || (s > 3.76 && n > 99))
    return 1-2*exp(-(2.000071+.331/sqrt(n)+1.409/n)*s);
  k = (int)(n*x)+1; m = 2*k-1; h = k-n*x;
  H = (double*)malloc((m*m)*sizeof(double));
  Q = (double*)malloc((m*m)*sizeof(double));
  for (i = 0; i < m; i++)
    for (j = 0; j < m; j++)
      if (i-j+1 < 0) H[i*m+j] = 0;
      else H[i*m+j] = 1;
  for (i = 0; i < m; i++)
    {H[i*m] -= pow(h, i+1);
    H[(m-1)*m+i] -= pow(h, (m-i));}
  H[(m-1)*m] += (2*h-1 > 0 ? pow(2*h-1, m) : 0);
  for (i = 0; i < m; i++)
    for (j = 0; j < m; j++)
      if (i-j+1 > 0) for (g = 1; g <= i-j+1; g++) H[i*m+j] /= g;
  eH = 0; mPower(H, eH, Q, &eQ, m, n);
  s = Q[(k-1)*m+k-1];
  for (i = 1; i <= n; i++)
    {s = s*i/n; if (s < 1e-140){s *= 1e140; eQ -= 140;}}
  s *= pow(10., eQ); /* attention au flottant 10.*/
  free(H);
  free(Q);
  return s;
}

double kolmogorovicdf(double p, int n){
  if (p > 1-accuracy) return 1;
  else if (p < accuracy) return 0;
  else{
    double lb, ub; lb = 0.0; ub = 1.0;
    while (ub-lb > accuracy){
      if (kolmogorovcdf((lb+ub)/2, n) < p) lb = (lb+ub)/2;
      else ub = (lb+ub)/2;
    }
    return (lb+ub)/2;
  }
}

/* Kolmogorov Limiting distribution (Dudley, 1964),
   with numerical adjustments (Stephens M.A., 1970) 
   $$
   \lim_{n\to\infty}\Pr\bigl\{K_n\leq u/\!\sqrt n\bigr\}
   =1+2\sum_{k=1}^\infty(-1)^k\exp\bigl(-2k^2u^2\bigr).
   $$
*/
/* to be defined? */
double klmpdf(double x, int n){return 0;}

double klmcdf(double x, int n){
  if (x <= 0) return 0;
  else if (x >= 1) return 1;
  else {
    double z, p, term; int k, pm;
    z = -2*pow(sqrt(n)+0.12+0.11/sqrt(n), 2)*x*x;
    p = 1.0; k = 1; pm = -1; term = 1;
    while (fabs(term) > accuracy && k < 100){
      term = 2*pm*exp(z*k*k);
      p += term; pm = -pm; k++;
  }
  return max(min(p, 1), 0);
  }
}

double klmicdf(double p, int n){
  if (p > 1-accuracy) return 1;
  else if (p < accuracy) return 0;
  else {
    double lb, ub; lb = 0.0; ub = 1.0;
    while (ub-lb > accuracy){
      if (klmcdf((lb+ub)/2, n) < p) lb = (lb+ub)/2;
      else ub = (lb+ub)/2;
    }
    return (lb+ub)/2;
  }
}

double kpmcdf(double x, int n) {
  if (x <= 0) return 0;
  else if (x >= 1) return 1;
  else {
    double p; int k; p = 0.0;
    for (k = n; k > n*x; k--)
      p = p+exp(lnGamma(n+1)-lnGamma(k+1)-lnGamma(n-k+1)
        +k*log(1.0*k/n-x)+(n-k-1)*log(x+1-1.0*k/n));
     return 1-x*p;
  }
}

double kpmicdf(double p, int n){
  if (p > 1-accuracy) return 1;
  else if (p < accuracy) return 0;
  else {
    double lb, ub; lb = 0.0; ub = 1.0;
    while (ub-lb > accuracy){
      if (kpmcdf((lb+ub)/2, n) < p) lb = (lb+ub)/2;
      else ub = (lb+ub)/2;
    }
    return (lb+ub)/2;
  }
}