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.
25 : *
26 : *
27 : *
28 : * This software is distributed on an AS IS basis, WITHOUT
29 : * WARRANTY OF ANY KIND, either express or implied.
30 : *
31 : ************************************************************************/
32 : #include <algorithm>
33 : #include <cmath>
34 : #include <limits>
35 : #include <numeric>
36 : #include <openGPMP/core/utils.hpp>
37 : #include <openGPMP/stats/pdfs.hpp>
38 : #include <random>
39 : #include <stdexcept>
40 : #include <vector>
41 :
42 : // Bernoulli distribution
43 : // Bernoulli distribution
44 0 : double gpmp::stats::PDF::bernoulli(double x, double p) {
45 0 : if (std::abs(x - 0) < std::numeric_limits<double>::epsilon())
46 0 : return 1 - p;
47 0 : else if (std::abs(x - 1) < std::numeric_limits<double>::epsilon())
48 0 : return p;
49 : else
50 0 : return 0;
51 : }
52 :
53 : // Beta distribution
54 0 : double gpmp::stats::PDF::beta(double x, double alpha, double beta) {
55 0 : if (x < 0 || x > 1)
56 0 : return 0;
57 : else
58 0 : return (std::pow(x, alpha - 1) * std::pow(1 - x, beta - 1)) /
59 0 : beta_function(alpha, beta);
60 : }
61 :
62 : // Binomial distribution
63 0 : double gpmp::stats::PDF::binomial(int k, int n, double p) {
64 0 : if (k < 0 || k > n)
65 0 : return 0;
66 : else
67 0 : return binomial_coefficient(n, k) * std::pow(p, k) *
68 0 : std::pow(1 - p, n - k);
69 : }
70 :
71 : // Cauchy distribution
72 0 : double gpmp::stats::PDF::cauchy(double x, double x0, double gamma) {
73 0 : return (1 / (M_PI * gamma)) *
74 0 : (gamma * gamma / ((x - x0) * (x - x0) + gamma * gamma));
75 : }
76 :
77 : // Chi-squared distribution
78 0 : double gpmp::stats::PDF::chi_squared(double x, int k) {
79 0 : if (x < 0)
80 0 : return 0;
81 : else
82 0 : return (std::pow(x, k / 2.0 - 1) * std::exp(-x / 2.0)) /
83 0 : (std::pow(2, k / 2.0) * std::tgamma(k / 2.0));
84 : }
85 :
86 : // Exponential distribution
87 0 : double gpmp::stats::PDF::exponential(double x, double lambda) {
88 0 : if (x < 0)
89 0 : return 0;
90 : else
91 0 : return lambda * std::exp(-lambda * x);
92 : }
93 :
94 : // F distribution
95 0 : double gpmp::stats::PDF::f_dist(double x, int df1, int df2) {
96 0 : if (x < 0)
97 0 : return 0;
98 : else
99 0 : return std::pow(df1, df1 / 2.0) * std::pow(df2, df2 / 2.0) *
100 0 : std::pow(x, df1 / 2.0 - 1) /
101 0 : (std::pow(df2 + df1 * x, (df1 + df2) / 2.0) *
102 0 : std::tgamma(df1 / 2.0) * std::tgamma(df2 / 2.0));
103 : }
104 :
105 : // Gamma distribution
106 0 : double gpmp::stats::PDF::gamma(double x, double alpha, double beta) {
107 0 : if (x < 0)
108 0 : return 0;
109 : else
110 0 : return (std::pow(x, alpha - 1) * std::exp(-x / beta)) /
111 0 : (std::pow(beta, alpha) * std::tgamma(alpha));
112 : }
113 :
114 : // Inverse-Gamma distribution
115 0 : double gpmp::stats::PDF::inverse_gamma(double x, double alpha, double beta) {
116 0 : if (x <= 0)
117 0 : return 0;
118 : else
119 0 : return (std::pow(beta, alpha) * std::pow(x, -alpha - 1) *
120 0 : std::exp(-beta / x)) /
121 0 : std::tgamma(alpha);
122 : }
123 :
124 : // Inverse-Gaussian distribution
125 0 : double gpmp::stats::PDF::inverse_gaussian(double x, double mu, double lambda) {
126 0 : if (x <= 0)
127 0 : return 0;
128 : else
129 0 : return std::sqrt(lambda / (2 * M_PI * x * x * x)) *
130 0 : std::exp(-lambda * (x - mu) * (x - mu) / (2 * mu * mu * x));
131 : }
132 :
133 : // Laplace distribution
134 0 : double gpmp::stats::PDF::laplace(double x, double mu, double b) {
135 0 : return 0.5 * std::exp(-std::abs(x - mu) / b) / b;
136 : }
137 :
138 : // Logistic distribution
139 0 : double gpmp::stats::PDF::logistic(double x, double mu, double s) {
140 0 : double z = (x - mu) / s;
141 0 : return std::exp(-z) / (s * std::pow(1 + std::exp(-z), 2));
142 : }
143 :
144 : // Log-Normal distribution
145 0 : double gpmp::stats::PDF::log_normal(double x, double mu, double sigma) {
146 0 : if (x <= 0)
147 0 : return 0;
148 : else
149 0 : return (1 / (x * sigma * std::sqrt(2 * M_PI))) *
150 0 : std::exp(-0.5 * std::pow((std::log(x) - mu) / sigma, 2));
151 : }
152 :
153 : // Normal (Gaussian) distribution
154 0 : double gpmp::stats::PDF::gaussian(double x, double mu, double sigma) {
155 0 : return (1 / (sigma * std::sqrt(2 * M_PI))) *
156 0 : std::exp(-0.5 * std::pow((x - mu) / sigma, 2));
157 : }
158 :
159 : // Poisson distribution
160 0 : double gpmp::stats::PDF::poisson(int k, double lambda) {
161 0 : if (k < 0)
162 0 : return 0;
163 : else
164 0 : return std::exp(-lambda) * std::pow(lambda, k) /
165 0 : gpmp::core::Misc::factorial(k);
166 : }
167 :
168 : // Rademacher distribution
169 0 : double gpmp::stats::PDF::rademacher(int k) {
170 0 : if (k == -1)
171 0 : return 0.5;
172 0 : else if (k == 1)
173 0 : return 0.5;
174 : else
175 0 : return 0;
176 : }
177 :
178 : // Student's t distribution
179 0 : double gpmp::stats::PDF::student_t(double x, int df) {
180 0 : double numerator = std::tgamma((df + 1) / 2.0);
181 0 : double denominator = std::sqrt(df * M_PI) * std::tgamma(df / 2.0);
182 0 : return std::pow(1 + x * x / df, -(df + 1) / 2.0) * numerator / denominator;
183 : }
184 :
185 : // Uniform distribution
186 0 : double gpmp::stats::PDF::uniform(double x, double a, double b) {
187 0 : if (x >= a && x <= b)
188 0 : return 1 / (b - a);
189 : else
190 0 : return 0;
191 : }
192 :
193 : // Weibull distribution
194 0 : double gpmp::stats::PDF::weibull(double x, double k, double lambda) {
195 0 : if (x < 0)
196 0 : return 0;
197 : else
198 0 : return (k / lambda) * std::pow(x / lambda, k - 1) *
199 0 : std::exp(-std::pow(x / lambda, k));
200 : }
201 :
202 : // Function to calculate binomial coefficient
203 0 : double gpmp::stats::PDF::binomial_coefficient(int n, int k) {
204 0 : double result = 1.0;
205 0 : for (int i = 1; i <= k; ++i) {
206 0 : result *= static_cast<double>(n - (k - i)) / i;
207 : }
208 0 : return result;
209 : }
210 :
211 : // Function to calculate beta function
212 0 : double gpmp::stats::PDF::beta_function(double alpha, double beta) {
213 0 : return std::tgamma(alpha) * std::tgamma(beta) / std::tgamma(alpha + beta);
214 : }
|
