1 /*
2 * Copyright (C) 2002, 2003 Laird Breyer
3 *
4 * This program is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program; if not, write to the Free Software
16 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
17 *
18 * Author: Laird Breyer <laird@lbreyer.com>
19 */
20
21 #include <math.h>
22 #include "util.h"
23
24 #ifndef M_PI
25 #define M_PI 3.14159265358979323846
26 #define M_SQRT2 1.41421356237309504880
27 #define M_2_SQRTPI 1.12837916709551257390
28 #endif
29
30
sample_mean(double x,double n)31 double sample_mean(double x, double n) {
32 return x/n;
33 }
34
sample_variance(double s,double x,double n)35 double sample_variance(double s, double x, double n) {
36 return (n * s - (x * x))/(n * (n - 1.0));
37 }
38
39 /* returns the logarithm of a Poisson distribution
40 (calculations based on Stirling's formula) */
log_poisson(int k,double lambda)41 double log_poisson(int k, double lambda) {
42
43 return (double)k * (log(lambda) -log((double)k) + 1.0) -
44 log(2.0 * M_PI * (double)k ) / 2.0 - lambda;
45
46 }
47
48 /* returns the probability that a chi squared with df
49 degrees of freedom is less than or equal to x */
chi2_cdf(double df,double x)50 double chi2_cdf(double df, double x) {
51 /* to be implemented */
52 return 1.0 - igamc(df/2.0, x/2.0);
53 }
54
gamma_tail(double a,double b,double x)55 double gamma_tail(double a, double b, double x) {
56 /* don't call igamc with extreme numbers, it can exit with an error */
57 return igamc(a, x/b);
58 }
59
normal_cdf(double x)60 double normal_cdf(double x) {
61 return ndtr(x);
62 }
63
64 /* Given X_1,...,X_n independent Gaussians, computes the probability
65 that X_k = min_j X_j.
66 eg:
67
68 Prob(X_1 <= min_j X_j) = \int \phi(x) \prod_{j\neq 1} Prob(X_j > x)
69
70 Call this function once for each X_j. The computation uses Gaussian
71 quadrature with Hermite polynomials, 10 points from Abramovitz and
72 Stegun, p.924. It doesn't need to be very accurate, but an error of
73 more than 1% looks bad when displayed to the user.
74 */
min_prob(int k,int n,double mu[],double sigma[])75 double min_prob(int k, int n, double mu[], double sigma[]) {
76 const double gauss_hermite[10][2] = {
77 {-3.4361591188, 1.0254516913},
78 {-2.5327316742, 0.8206661264},
79 {-1.7566836492, 0.7414419319},
80 {-1.0366108297 , 0.7032963231},
81 {-0.3429013272, 0.6870818539},
82 {0.3429013272, 0.6870818539},
83 {1.0366108297 , 0.7032963231},
84 {1.7566836492, 0.7414419319},
85 {2.5327316742, 0.8206661264},
86 {3.4361591188, 1.0254516913},
87 };
88 double x, p, srt, f, g;
89 int i,j;
90
91 srt = sigma[k] * M_SQRT2;
92 p = 0.0;
93 for(i = 0; i < 10; i++) {
94 x = gauss_hermite[i][0];
95 f = exp(-x*x);
96 g = f;
97 /* don't compute full product if too small to matter (less than 1%) */
98 for(j = 0; (j < n) && (g > 0.01); j++) {
99 if( j == k ) continue;
100 g *= (1.0 - ndtr( ((mu[k] + x * srt) - mu[j])/sigma[j] ));
101 }
102 p += gauss_hermite[i][1] * g;
103 }
104 return p * M_2_SQRTPI / 2.0;
105 }
106