samtools/bcftools/kfunc.c.pysam.c

   1 #include "pysam.h"
   2
   3 #include <math.h>
   4
   5
   6 /* Log gamma function
   7  * \log{\Gamma(z)}
   8  * AS245, 2nd algorithm, http://lib.stat.cmu.edu/apstat/245
   9  */
  10 double kf_lgamma(double z)
  11 {
  12         double x = 0;
  13         x += 0.1659470187408462e-06 / (z+7);
  14         x += 0.9934937113930748e-05 / (z+6);
  15         x -= 0.1385710331296526     / (z+5);
  16         x += 12.50734324009056      / (z+4);
  17         x -= 176.6150291498386      / (z+3);
  18         x += 771.3234287757674      / (z+2);
  19         x -= 1259.139216722289      / (z+1);
  20         x += 676.5203681218835      / z;
  21         x += 0.9999999999995183;
  22         return log(x) - 5.58106146679532777 - z + (z-0.5) * log(z+6.5);
  23 }
  24
  25 /* complementary error function
  26  * \frac{2}{\sqrt{\pi}} \int_x^{\infty} e^{-t^2} dt
  27  * AS66, 2nd algorithm, http://lib.stat.cmu.edu/apstat/66
  28  */
  29 double kf_erfc(double x)
  30 {
  31         const double p0 = 220.2068679123761;
  32         const double p1 = 221.2135961699311;
  33         const double p2 = 112.0792914978709;
  34         const double p3 = 33.912866078383;
  35         const double p4 = 6.37396220353165;
  36         const double p5 = .7003830644436881;
  37         const double p6 = .03526249659989109;
  38         const double q0 = 440.4137358247522;
  39         const double q1 = 793.8265125199484;
  40         const double q2 = 637.3336333788311;
  41         const double q3 = 296.5642487796737;
  42         const double q4 = 86.78073220294608;
  43         const double q5 = 16.06417757920695;
  44         const double q6 = 1.755667163182642;
  45         const double q7 = .08838834764831844;
  46         double expntl, z, p;
  47         z = fabs(x) * M_SQRT2;
  48         if (z > 37.) return x > 0.? 0. : 2.;
  49         expntl = exp(z * z * - .5);
  50         if (z < 10. / M_SQRT2) // for small z
  51             p = expntl * ((((((p6 * z + p5) * z + p4) * z + p3) * z + p2) * z + p1) * z + p0)
  52                         / (((((((q7 * z + q6) * z + q5) * z + q4) * z + q3) * z + q2) * z + q1) * z + q0);
  53         else p = expntl / 2.506628274631001 / (z + 1. / (z + 2. / (z + 3. / (z + 4. / (z + .65)))));
  54         return x > 0.? 2. * p : 2. * (1. - p);
  55 }
  56
  57 /* The following computes regularized incomplete gamma functions.
  58  * Formulas are taken from Wiki, with additional input from Numerical
  59  * Recipes in C (for modified Lentz's algorithm) and AS245
  60  * (http://lib.stat.cmu.edu/apstat/245).
  61  *
  62  * A good online calculator is available at:
  63  *
  64  *   http://www.danielsoper.com/statcalc/calc23.aspx
  65  *
  66  * It calculates upper incomplete gamma function, which equals
  67  * kf_gammaq(s,z)*tgamma(s).
  68  */
  69
  70 #define KF_GAMMA_EPS 1e-14
  71 #define KF_TINY 1e-290
  72
  73 // regularized lower incomplete gamma function, by series expansion
  74 static double _kf_gammap(double s, double z)
  75 {
  76         double sum, x;
  77         int k;
  78         for (k = 1, sum = x = 1.; k < 100; ++k) {
  79                 sum += (x *= z / (s + k));
  80                 if (x / sum < KF_GAMMA_EPS) break;
  81         }
  82         return exp(s * log(z) - z - kf_lgamma(s + 1.) + log(sum));
  83 }
  84 // regularized upper incomplete gamma function, by continued fraction
  85 static double _kf_gammaq(double s, double z)
  86 {
  87         int j;
  88         double C, D, f;
  89         f = 1. + z - s; C = f; D = 0.;
  90         // Modified Lentz's algorithm for computing continued fraction
  91         // See Numerical Recipes in C, 2nd edition, section 5.2
  92         for (j = 1; j < 100; ++j) {
  93                 double a = j * (s - j), b = (j<<1) + 1 + z - s, d;
  94                 D = b + a * D;
  95                 if (D < KF_TINY) D = KF_TINY;
  96                 C = b + a / C;
  97                 if (C < KF_TINY) C = KF_TINY;
  98                 D = 1. / D;
  99                 d = C * D;
 100                 f *= d;
 101                 if (fabs(d - 1.) < KF_GAMMA_EPS) break;
 102         }
 103         return exp(s * log(z) - z - kf_lgamma(s) - log(f));
 104 }
 105
 106 double kf_gammap(double s, double z)
 107 {
 108         return z <= 1. || z < s? _kf_gammap(s, z) : 1. - _kf_gammaq(s, z);
 109 }
 110
 111 double kf_gammaq(double s, double z)
 112 {
 113         return z <= 1. || z < s? 1. - _kf_gammap(s, z) : _kf_gammaq(s, z);
 114 }
 115
 116 /* Regularized incomplete beta function. The method is taken from
 117  * Numerical Recipe in C, 2nd edition, section 6.4. The following web
 118  * page calculates the incomplete beta function, which equals
 119  * kf_betai(a,b,x) * gamma(a) * gamma(b) / gamma(a+b):
 120  *
 121  *   http://www.danielsoper.com/statcalc/calc36.aspx
 122  */
 123 static double kf_betai_aux(double a, double b, double x)
 124 {
 125         double C, D, f;
 126         int j;
 127         if (x == 0.) return 0.;
 128         if (x == 1.) return 1.;
 129         f = 1.; C = f; D = 0.;
 130         // Modified Lentz's algorithm for computing continued fraction
 131         for (j = 1; j < 200; ++j) {
 132                 double aa, d;
 133                 int m = j>>1;
 134                 aa = (j&1)? -(a + m) * (a + b + m) * x / ((a + 2*m) * (a + 2*m + 1))
 135                         : m * (b - m) * x / ((a + 2*m - 1) * (a + 2*m));
 136                 D = 1. + aa * D;
 137                 if (D < KF_TINY) D = KF_TINY;
 138                 C = 1. + aa / C;
 139                 if (C < KF_TINY) C = KF_TINY;
 140                 D = 1. / D;
 141                 d = C * D;
 142                 f *= d;
 143                 if (fabs(d - 1.) < KF_GAMMA_EPS) break;
 144         }
 145         return exp(kf_lgamma(a+b) - kf_lgamma(a) - kf_lgamma(b) + a * log(x) + b * log(1.-x)) / a / f;
 146 }
 147 double kf_betai(double a, double b, double x)
 148 {
 149         return x < (a + 1.) / (a + b + 2.)? kf_betai_aux(a, b, x) : 1. - kf_betai_aux(b, a, 1. - x);
 150 }
 151
 152 #ifdef KF_MAIN
 153 #include <stdio.h>
 154 int main(int argc, char *argv[])
 155 {
 156         double x = 5.5, y = 3;
 157         double a, b;
 158         printf("erfc(%lg): %lg, %lg\n", x, erfc(x), kf_erfc(x));
 159         printf("upper-gamma(%lg,%lg): %lg\n", x, y, kf_gammaq(y, x)*tgamma(y));
 160         a = 2; b = 2; x = 0.5;
 161         printf("incomplete-beta(%lg,%lg,%lg): %lg\n", a, b, x, kf_betai(a, b, x) / exp(kf_lgamma(a+b) - kf_lgamma(a) - kf_lgamma(b)));
 162         return 0;
 163 }
 164 #endif