Theoretica
A C++ numerical and automatic mathematical library
Loading...
Searching...
No Matches
distributions.h
Go to the documentation of this file.
1
5
6#ifndef THEORETICA_DISTRIBUTIONS
7#define THEORETICA_DISTRIBUTIONS
8
9#include "../core/real_analysis.h"
10#include "../algebra/vec.h"
11#include "./statistics.h"
12#include "../core/function.h"
13#include "../core/special.h"
14
15
16namespace theoretica {
17
18
19 namespace stats {
20
27 inline real likelihood(const vec<real>& X, const vec<real>& theta, stat_function f) {
28
29 real res = 1;
30
31 for (unsigned int i = 0; i < X.size(); ++i)
32 res *= f(X[i], theta);
33
34 return res;
35 }
36
37
44 inline real log_likelihood(const vec<real>& X, const vec<real>& theta, stat_function f) {
45
46 real res = 0;
47
48 for (unsigned int i = 0; i < X.size(); ++i)
49 res += ln(f(X[i], theta));
50
51 return res;
52 }
53 }
54
55
57 namespace distribution {
58
59
62
63 return (1.0 / (sigma * SQRT2 * SQRTPI))
64 * exp(-square(x - X) / (2 * square(sigma)));
65 }
66
67
69 inline real gaussian(real x, const vec<real>& theta) {
70
71 if(theta.size() != 2) {
72 TH_MATH_ERROR("distribution::gaussian", theta.size(), INVALID_ARGUMENT);
73 return nan();
74 }
75
76 return gaussian(x, theta[0], theta[1]);
77 }
78
79
81 inline real bernoulli(unsigned int k, real p) {
82
83 return pow(p, k) * pow(1 - p, 1 - k);
84 }
85
86
88 inline real bernoulli(real k, const vec<real>& theta) {
89
90 if(theta.size() != 1) {
91 TH_MATH_ERROR("distribution::bernoulli", theta.size(), INVALID_ARGUMENT);
92 return nan();
93 }
94
95 return bernoulli(static_cast<unsigned int>(k), theta[0]);
96
97 }
98
99
101 inline real poisson(unsigned int k, real lambda) {
102
103 return exp(-lambda) * pow(lambda, k) / (real) fact(k);
104 }
105
106
108 inline real poisson(real k, const vec<real>& theta) {
109
110 if(theta.size() != 1) {
111 TH_MATH_ERROR("distribution::poisson", theta.size(), INVALID_ARGUMENT);
112 return nan();
113 }
114
115 return poisson(static_cast<unsigned int>(k), theta[0]);
116 }
117
118
120 inline real binomial(unsigned int nu, unsigned int n, real p) {
121
122 return binomial_coeff(n, nu) *
123 pow(p, nu) * pow(1 - p, n - nu);
124 }
125
126
128 inline real binomial(real nu, const vec<real>& theta) {
129
130 if(theta.size() != 2) {
131 TH_MATH_ERROR("distribution::binomial", theta.size(), INVALID_ARGUMENT);
132 return nan();
133 }
134
135 return binomial(static_cast<unsigned int>(nu),
136 static_cast<unsigned int>(theta[0]), theta[1]);
137 }
138
139
142 std::vector<unsigned int> x,
143 unsigned int n,
144 unsigned int k,
145 std::vector<real> p) {
146
147 if(x.size() != p.size() || x.size() != k) {
148 TH_MATH_ERROR("distribution::multinomial", x.size(), INVALID_ARGUMENT);
149 return nan();
150 }
151
152 real res = fact(n);
153
154 for (unsigned int i = 0; i < k; ++i)
155 res *= pow(p[i], x[i]) / fact(x[i]);
156
157 return res;
158 }
159
160
165 inline real chi_squared(real x, unsigned int k) {
166
167 if(x < 0)
168 return 0;
169
170 const real coeff = special::half_gamma(k);
171
172 if(k % 2 == 0)
173 return pow(x, int(k / 2) - 1) * exp(-x / 2.0)
174 / (pow(SQRT2, k) * coeff);
175 else
176 return pow(sqrt(x), k - 2) * exp(-x / 2.0)
177 / (pow(SQRT2, k) * coeff);
178 }
179
180
190 inline real chi_squared(real x, unsigned int k, real half_gamma_k) {
191
192 if(x < 0)
193 return 0;
194
195 if(k % 2 == 0)
196 return pow(x, int(k / 2) - 1) * exp(-x / 2.0)
197 / (pow(SQRT2, k) * half_gamma_k);
198 else
199 return pow(sqrt(x), k - 2) * exp(-x / 2.0)
200 / (pow(SQRT2, k) * half_gamma_k);
201 }
202
203
205 inline real chi_squared(real x, const vec<real>& theta) {
206
207 if(theta.size() != 1) {
208 TH_MATH_ERROR("distribution::chi_squared", theta.size(), INVALID_ARGUMENT);
209 return nan();
210 }
211
212 return chi_squared(x, static_cast<unsigned int>(theta[0]));
213 }
214
215
219
220 return powf(beta, alpha) * powf(x, alpha - 1)
221 * exp(-beta * x) / special::gamma(alpha);
222 }
223
224
226 inline real gamma(real x, const vec<real>& theta) {
227
228 if(theta.size() != 2) {
229 TH_MATH_ERROR("distribution::gamma", theta.size(), INVALID_ARGUMENT);
230 return nan();
231 }
232
233 return gamma(x, theta[0], theta[1]);
234 }
235
236
240
241 return powf(x, alpha - 1) * powf(1 - x, beta - 1)
243 }
244
245
247 inline real beta(real x, const vec<real>& theta) {
248
249 if(theta.size() != 2) {
250 TH_MATH_ERROR("distribution::beta", theta.size(), INVALID_ARGUMENT);
251 return nan();
252 }
253
254 return beta(x, theta[0], theta[1]);
255 }
256
257
259 inline real student(real x, unsigned int nu) {
260
261 const real a = 1 + (x * x / nu);
262
264 * pow(sqrt(a), -nu - 1) / (sqrt(nu) * SQRTPI);
265 }
266
267
269 inline real student(real x, const vec<real>& theta) {
270
271 if(theta.size() != 1) {
272 TH_MATH_ERROR("distribution::student", theta.size(), INVALID_ARGUMENT);
273 return nan();
274 }
275
276 return student(x, static_cast<unsigned int>(theta[0]));
277 }
278
279
282
283 return 1.0 / (2.50662827463 * sigma * x) *
284 exp(-square(ln(x) - mu) / (2 * square(sigma)));
285 }
286
287
289 inline real log_normal(real x, const vec<real>& theta) {
290
291 if(theta.size() != 2) {
292 TH_MATH_ERROR("distribution::log_normal", theta.size(), INVALID_ARGUMENT);
293 return nan();
294 }
295
296 return log_normal(x, theta[0], theta[1]);
297 }
298
299
302
303 if(x < 0)
304 return 0;
305
306 return lambda * exp(-lambda * x);
307 }
308
309
311 inline real exponential(real x, const vec<real>& theta) {
312
313 if(theta.size() != 1) {
314 TH_MATH_ERROR("distribution::exponential", theta.size(), INVALID_ARGUMENT);
315 return nan();
316 }
317
318 return exponential(x, theta[0]);
319 }
320
321
324
325 if(x < 0)
326 return 0;
327
328 if(sigma < MACH_EPSILON) {
329 TH_MATH_ERROR("distribution::rayleigh", sigma, DIV_BY_ZERO);
330 return nan();
331 }
332
333 return x * exp(-square(x / sigma) / 2) / square(sigma);
334 }
335
336
338 inline real rayleigh(real x, const vec<real>& theta) {
339
340 if(theta.size() != 1) {
341 TH_MATH_ERROR("distribution::rayleigh", theta.size(), INVALID_ARGUMENT);
342 return nan();
343 }
344
345 return rayleigh(x, theta[0]);
346 }
347
348
351
352 return 1.0 / (PI * alpha * (1 + (square(x - mu) / square(alpha))));
353 }
354
355
357 inline real cauchy(real x, const vec<real>& theta) {
358
359 if(theta.size() != 2) {
360 TH_MATH_ERROR("distribution::cauchy", theta.size(), INVALID_ARGUMENT);
361 return nan();
362 }
363
364 return cauchy(x, theta[0], theta[1]);
365 }
366
367
370
371 return Gamma / (2 * PI * (square(x - M) + square(Gamma / 2.0)));
372 }
373
374
376 inline real breit_wigner(real x, const vec<real>& theta) {
377
378 if(theta.size() != 2) {
379 TH_MATH_ERROR("distribution::breit_wigner", theta.size(), INVALID_ARGUMENT);
380 return nan();
381 }
382
383 return breit_wigner(x, theta[0], theta[1]);
384 }
385
386
388 inline real maxwell(real x, real a) {
389
390 return (SQRT2 / SQRTPI) * square(x) * exp(-square(x / a) / 2.0) / cube(a);
391 }
392
393
395 inline real maxwell(real x, const vec<real>& theta) {
396
397 if(theta.size() != 1) {
398 TH_MATH_ERROR("distribution::maxwell", theta.size(), INVALID_ARGUMENT);
399 return nan();
400 }
401
402 return maxwell(x, theta[0]);
403 }
404
405
407 inline real laplace(real x, real mu, real b) {
408
409 return (1.0 / (2.0 * b)) * exp(-abs(x - mu) / b);
410 }
411
412
414 inline real laplace(real x, const vec<real>& theta) {
415
416 if(theta.size() != 2) {
417 TH_MATH_ERROR("distribution::laplace", theta.size(), INVALID_ARGUMENT);
418 return nan();
419 }
420
421 return (1.0 / (2.0 * theta[1])) * exp(-abs(x - theta[0]) / theta[1]);
422 }
423
424
427
428 if(alpha <= 0) {
429 TH_MATH_ERROR("distribution::pareto", alpha, OUT_OF_DOMAIN);
430 return nan();
431 }
432
433 if(x < x_m)
434 return 0;
435
436 return alpha * powf(x_m, alpha) / powf(x, alpha + 1);
437 }
438
439
441 inline real pareto(real x, const vec<real>& theta) {
442
443 if(theta.size() != 2) {
444 TH_MATH_ERROR("distribution::pareto", theta.size(), INVALID_ARGUMENT);
445 return nan();
446 }
447
448 return pareto(x, theta[0], theta[1]);
449 }
450
451
453 inline real erlang(real x, unsigned int k, real lambda) {
454
455 if(!k) {
456 TH_MATH_ERROR("distribution::erlang", k, INVALID_ARGUMENT);
457 return nan();
458 }
459
460 return pow(lambda, k) * pow(x, k - 1)
461 * exp(-lambda * x) / fact(k - 1);
462 }
463
464
466 inline real erlang(real x, const vec<real>& theta) {
467
468 if(theta.size() != 2) {
469 TH_MATH_ERROR("distribution::erlang", theta.size(), INVALID_ARGUMENT);
470 return nan();
471 }
472
473 return erlang(x, static_cast<unsigned int>(theta[0]), theta[1]);
474 }
475
476 }
477}
478
479#endif
TH_CONSTEXPR unsigned int size() const
Returns the size of the vector (N)
Definition vec.h:406
#define TH_MATH_ERROR(F_NAME, VALUE, EXCEPTION)
TH_MATH_ERROR is a macro which throws exceptions or modifies errno (depending on which compiling opti...
Definition error.h:225
real maxwell(real x, real a)
Maxwell-Boltzmann distribution.
Definition distributions.h:388
real binomial(unsigned int nu, unsigned int n, real p)
Binomial distribution function.
Definition distributions.h:120
real exponential(real x, real lambda)
Exponential distribution.
Definition distributions.h:301
real poisson(unsigned int k, real lambda)
Poisson distribution.
Definition distributions.h:101
real pareto(real x, real x_m, real alpha)
Pareto distribution.
Definition distributions.h:426
real bernoulli(unsigned int k, real p)
Bernoulli distribution.
Definition distributions.h:81
real beta(real x, real alpha, real beta)
Beta distribution density function with parameters alpha (shape) and beta (shape)
Definition distributions.h:239
real breit_wigner(real x, real M, real Gamma)
Breit Wigner distribution.
Definition distributions.h:369
real chi_squared(real x, unsigned int k)
Chi-squared distribution.
Definition distributions.h:165
real erlang(real x, unsigned int k, real lambda)
Erlang distribution.
Definition distributions.h:453
real laplace(real x, real mu, real b)
Laplace distribution.
Definition distributions.h:407
real student(real x, unsigned int nu)
Student's t distribution.
Definition distributions.h:259
real log_normal(real x, real mu, real sigma)
Log-normal distribution.
Definition distributions.h:281
real multinomial(std::vector< unsigned int > x, unsigned int n, unsigned int k, std::vector< real > p)
Multinomial distribution.
Definition distributions.h:141
real gamma(real x, real alpha, real beta)
Gamma distribution density function with parameters alpha (shape) and beta (rate)
Definition distributions.h:218
real cauchy(real x, real mu, real alpha)
Cauchy distribution.
Definition distributions.h:350
real gaussian(real x, real X, real sigma)
Gaussian distribution function.
Definition distributions.h:61
real rayleigh(real x, real sigma)
Rayleigh distribution.
Definition distributions.h:323
real beta(real x1, real x2)
Beta special function of real argument.
Definition special.h:141
real half_gamma(unsigned int k)
Half Gamma special function, defined as HG(n) = Gamma(n / 2) for any positive integer n.
Definition special.h:41
real gamma(unsigned int k)
Gamma special function of positive integer argument.
Definition special.h:23
real likelihood(const vec< real > &X, const vec< real > &theta, stat_function f)
Compute the likelihood of a distribution <f> with the given parameters <theta> and measures <X>
Definition distributions.h:27
real log_likelihood(const vec< real > &X, const vec< real > &theta, stat_function f)
Compute the log likelihood of a distribution <f> with the given parameters <theta> and measures <X>
Definition distributions.h:44
Main namespace of the library which contains all functions and objects.
Definition algebra.h:27
double real
A real number, defined as a floating point type.
Definition constants.h:198
TH_CONSTEXPR IntType binomial_coeff(unsigned int n, unsigned int m)
Compute the binomial coefficient.
Definition real_analysis.h:1229
dual2 sqrt(dual2 x)
Compute the square root of a second order dual number.
Definition dual2_functions.h:54
dual2 ln(dual2 x)
Compute the natural logarithm of a second order dual number.
Definition dual2_functions.h:151
dual2 abs(dual2 x)
Compute the absolute value of a second order dual number.
Definition dual2_functions.h:198
std::remove_reference_t< decltype(std::declval< Structure >()[0])> vector_element_t
Extract the type of a vector (or any indexable container) from its operator[].
Definition core_traits.h:134
constexpr real SQRTPI
The square root of Pi.
Definition constants.h:234
dual2 exp(dual2 x)
Compute the exponential of a second order dual number.
Definition dual2_functions.h:138
real powf(real x, real a)
Approximate x elevated to a real exponent.
Definition real_analysis.h:795
TH_CONSTEXPR real nan()
Return a quiet NaN number in floating point representation.
Definition error.h:54
TH_CONSTEXPR IntType fact(unsigned int n)
Compute the factorial of n.
Definition real_analysis.h:670
constexpr real MACH_EPSILON
Machine epsilon for the real type.
Definition constants.h:207
std::function< real(real, const vec< real > &)> stat_function
Function pointer to a probability distribution function where the first argument is the variable and ...
Definition function.h:30
constexpr real SQRT2
The square root of 2.
Definition constants.h:261
dual2 square(dual2 x)
Return the square of a second order dual number.
Definition dual2_functions.h:23
constexpr real PI
The Pi mathematical constant.
Definition constants.h:216
dual2 pow(dual2 x, int n)
Compute the n-th power of a second order dual number.
Definition dual2_functions.h:41
dual2 cube(dual2 x)
Return the cube of a second order dual number.
Definition dual2_functions.h:29
Statistical functions.