Line data Source code
1 : /*************************************************************************
2 : *
3 : * Project
4 : * _____ _____ __ __ _____
5 : * / ____| __ \| \/ | __ \
6 : * ___ _ __ ___ _ __ | | __| |__) | \ / | |__) |
7 : * / _ \| '_ \ / _ \ '_ \| | |_ | ___/| |\/| | ___/
8 : *| (_) | |_) | __/ | | | |__| | | | | | | |
9 : * \___/| .__/ \___|_| |_|\_____|_| |_| |_|_|
10 : * | |
11 : * |_|
12 : *
13 : * Copyright (C) Akiel Aries, <akiel@akiel.org>, et al.
14 : *
15 : * This software is licensed as described in the file LICENSE, which
16 : * you should have received as part of this distribution. The terms
17 : * among other details are referenced in the official documentation
18 : * seen here : https://akielaries.github.io/openGPMP/ along with
19 : * important files seen in this project.
20 : *
21 : * You may opt to use, copy, modify, merge, publish, distribute
22 : * and/or sell copies of the Software, and permit persons to whom
23 : * the Software is furnished to do so, under the terms of the
24 : * LICENSE file. As this is an Open Source effort, all implementations
25 : * must be of the same methodology.
26 : *
27 : *
28 : *
29 : * This software is distributed on an AS IS basis, WITHOUT
30 : * WARRANTY OF ANY KIND, either express or implied.
31 : *
32 : ************************************************************************/
33 : #include <algorithm>
34 : #include <cmath>
35 : #include <limits>
36 : #include <numeric>
37 : #include <openGPMP/stats/cdfs.hpp>
38 : #include <random>
39 : #include <stdexcept>
40 : #include <vector>
41 :
42 : // Bernoulli CDF
43 0 : double gpmp::stats::CDF::bernoulli(double x, double p) {
44 0 : if (x < 0)
45 0 : return 0.0;
46 0 : else if (x < 1)
47 0 : return 1 - (1 - p);
48 : else
49 0 : return 1.0;
50 : }
51 :
52 : // Beta CDF using the incomplete beta function
53 0 : double gpmp::stats::CDF::beta(double x, double alpha, double beta) {
54 0 : if (x <= 0)
55 0 : return 0.0;
56 0 : else if (x >= 1)
57 0 : return 1.0;
58 : else
59 0 : return incomplete_beta(alpha, beta, x);
60 : }
61 :
62 : // Binomial CDF
63 0 : double gpmp::stats::CDF::binomial(int k, int n, double p) {
64 0 : if (k < 0)
65 0 : return 0.0;
66 0 : else if (k >= n)
67 0 : return 1.0;
68 : else
69 0 : return incomplete_beta(1.0 - p, n - k, k + 1);
70 : }
71 :
72 : // Cauchy CDF
73 0 : double gpmp::stats::CDF::cauchy(double x, double x0, double gamma) {
74 0 : return 0.5 + atan((x - x0) / gamma) / M_PI;
75 : }
76 :
77 : // Chi-squared CDF
78 0 : double gpmp::stats::CDF::chi_squared(double x, double k) {
79 0 : if (x < 0)
80 0 : return 0.0;
81 : else
82 0 : return incomplete_gamma(k / 2.0, x / 2.0);
83 : }
84 :
85 : // Exponential CDF
86 0 : double gpmp::stats::CDF::exponential(double x, double lambda) {
87 0 : if (x < 0)
88 0 : return 0.0;
89 : else
90 0 : return 1.0 - exp(-lambda * x);
91 : }
92 :
93 : // F CDF
94 0 : double gpmp::stats::CDF::f(double x, double df1, double df2) {
95 0 : if (x <= 0)
96 0 : return 0.0;
97 : else
98 0 : return incomplete_beta(df1 / 2.0, df2 / 2.0, df1 / (df1 + df2 * x));
99 : }
100 :
101 : // Gamma CDF
102 0 : double gpmp::stats::CDF::gamma(double x, double shape, double scale) {
103 0 : if (x < 0)
104 0 : return 0.0;
105 : else
106 0 : return incomplete_gamma(shape, x / scale);
107 : }
108 :
109 : // Inverse-Gamma CDF
110 0 : double gpmp::stats::CDF::inverse_gamma(double x, double shape, double scale) {
111 0 : if (x <= 0)
112 0 : return 0.0;
113 : else
114 0 : return 1.0 - incomplete_gamma(shape, scale / x);
115 : }
116 :
117 : // Inverse-Gaussian CDF
118 0 : double gpmp::stats::CDF::inverse_gaussian(double x, double mu, double lambda) {
119 0 : if (x <= 0)
120 0 : return 0.0;
121 : else
122 0 : return normal_cdf(sqrt(lambda / x) * (x / mu - 1.0));
123 : }
124 :
125 : // Laplace CDF
126 0 : double gpmp::stats::CDF::laplace(double x, double mu, double b) {
127 0 : if (x < mu)
128 0 : return 0.5 * exp((x - mu) / b);
129 : else
130 0 : return 1.0 - 0.5 * exp(-(x - mu) / b);
131 : }
132 :
133 : // Logistic CDF
134 0 : double gpmp::stats::CDF::logistic(double x, double mu, double s) {
135 0 : return 1.0 / (1.0 + exp(-(x - mu) / s));
136 : }
137 :
138 : // Log-Normal CDF
139 0 : double gpmp::stats::CDF::log_normal(double x, double mu, double sigma) {
140 0 : if (x <= 0)
141 0 : return 0.0;
142 : else
143 0 : return 0.5 + 0.5 * erf((log(x) - mu) / (sqrt(2.0) * sigma));
144 : }
145 :
146 : // Normal (Gaussian) CDF
147 0 : double gpmp::stats::CDF::gaussian(double x, double mu, double sigma) {
148 0 : return 0.5 * (1 + erf((x - mu) / (sigma * sqrt(2))));
149 : }
150 :
151 : // Poisson CDF
152 0 : double gpmp::stats::CDF::poisson(int k, double lambda) {
153 0 : if (k < 0)
154 0 : return 0.0;
155 : else
156 0 : return incomplete_gamma(k + 1, lambda);
157 : }
158 :
159 : // Rademacher CDF
160 0 : double gpmp::stats::CDF::rademacher(double x) {
161 0 : if (x < 0)
162 0 : return 0.0;
163 0 : else if (x < 0.5)
164 0 : return 0.0;
165 : else
166 0 : return 1.0;
167 : }
168 :
169 : // Student's t CDF
170 0 : double gpmp::stats::CDF::student_t(double x, double df) {
171 0 : if (df <= 0)
172 0 : return NAN;
173 0 : return 0.5 + 0.5 * std::tgamma((df + 1) / 2) * std::hypot(1, x / sqrt(df)) /
174 0 : (sqrt(df) * std::tgamma(df / 2));
175 : }
176 :
177 : // Uniform CDF
178 0 : double gpmp::stats::CDF::uniform(double x, double a, double b) {
179 0 : if (x < a)
180 0 : return 0.0;
181 0 : else if (x < b)
182 0 : return (x - a) / (b - a);
183 : else
184 0 : return 1.0;
185 : }
186 :
187 : // Weibull CDF
188 0 : double gpmp::stats::CDF::weibull(double x, double shape, double scale) {
189 0 : if (x < 0)
190 0 : return 0.0;
191 : else
192 0 : return 1.0 - exp(-pow(x / scale, shape));
193 : }
194 :
195 : // Normal cumulative distribution function
196 0 : double gpmp::stats::CDF::normal_cdf(double x) {
197 0 : return 0.5 * (1 + erf(x / sqrt(2.0)));
198 : }
199 :
200 : // Incomplete beta function for beta and binomial CDF
201 0 : double gpmp::stats::CDF::incomplete_beta(double a, double b, double x) {
202 0 : const int maxIterations = 1000;
203 0 : const double epsilon = 1e-12;
204 0 : double result = 0.0;
205 0 : double term = 1.0;
206 0 : for (int k = 0; k < maxIterations; ++k) {
207 0 : term *=
208 0 : (k == 0 ? 1.0
209 0 : : (a + k - 1) * (b - k) / (a + b + k - 1) * x / (k + 1));
210 0 : result += term;
211 0 : if (std::abs(term) < epsilon * std::abs(result))
212 0 : break;
213 : }
214 0 : return result * std::pow(x, a) * std::pow(1 - x, b) / (a * std::beta(a, b));
215 : }
216 :
217 : // Incomplete gamma function for chi-squared and gamma CDFs
218 0 : double gpmp::stats::CDF::incomplete_gamma(double a, double x) {
219 0 : const int maxIterations = 1000;
220 0 : const double epsilon = 1e-12;
221 0 : double result = 0.0;
222 0 : double term = 1.0;
223 0 : for (int k = 0; k < maxIterations; ++k) {
224 0 : term *= x / (a + k);
225 0 : result += term;
226 0 : if (std::abs(term) < epsilon * std::abs(result))
227 0 : break;
228 : }
229 0 : return exp(-x) * std::pow(x, a) * result / std::tgamma(a);
230 : }
|
