atmosphere.moment source code

1 /++
2 Moment of probability distributions.
3 +/
4 /**
5 Authors: [Ilya Yaroshenko](http://9il.github.io)
6 
7 Copyright: © 2014-2015 [Ilya Yaroshenko](http://9il.github.io)
8 
9 License: MIT
10 */
11 module atmosphere.moment;
12 
13 import atmosphere.math: besselKR;
14 
15 import core.stdc.tgmath;
16 
17 import std.typecons;
18 
19 
20 /++
21 Moments of the generalized gamma distribution
22 +/
23 T generalizedGammaMean(T)(T shape, T scale, T power, uint order = 1)
24 {
25 	import std.math : pow;
26 	return scale.pow(order) * (tgamma(shape+order/power) / tgamma(shape));
27 }
28 
29 ///
30 unittest
31 {
32 	auto expectation     = generalizedGammaMean!double(2.0, 3.0, 4.0);
33 	auto secondRowMoment = generalizedGammaMean!double(2.0, 3.0, 4.0, 2);
34 }
35 
36 /++
37 Moments of the generalized variance-gamma distribution
38 +/
39 T generalizedVarianceGammaMean(T)(T shape, T scale, T power, T beta)
40 {
41 	return beta * generalizedGammaMean(shape, scale, power);
42 }
43 
44 ///
45 unittest
46 {
47 	auto expectation = generalizedVarianceGammaMean!double(2.0, 3.0, 4.0, 5.0);
48 }
49 
50 
51 /++
52 Moments of the generalized inverse Gaussian distribution
53 +/
54 T properGeneralizedInverseGaussianMean(T)(T lambda, T eta, T omega, uint order = 1)
55 in {
56 	assert(order);
57 }
58 body {
59 	import std.math : pow;
60 	return eta.pow(order) * besselKR(lambda, omega, order-1);
61 }
62 
63 ///
64 unittest
65 {
66 	auto expectation     = properGeneralizedInverseGaussianMean!double(2.0, 3.0, 4.0);
67 	auto secondRowMoment = properGeneralizedInverseGaussianMean!double(2.0, 3.0, 4.0, 2);
68 }
69 
70 
71 /++
72 Variance of the generalized inverse Gaussian distribution
73 +/
74 //(besselK[x-1, y] * besselK[x+1, y] - besselK[x, y]^2) / besselK[x-1, y]^2
75 T properGeneralizedInverseGaussianVariance(T)(T lambda, T eta, T omega)
76 {
77 	immutable r = besselKR(lambda, omega);
78 	return eta^^2 * (r * (2 * (lambda + 1) / omega - r) + 1);
79 }
80 
81 ///
82 unittest
83 {
84 	import std.math : approxEqual;
85 	auto expectation     = properGeneralizedInverseGaussianMean!double(2.0, 3.0, 4.0);
86 	auto secondRowMoment = properGeneralizedInverseGaussianMean!double(2.0, 3.0, 4.0, 2);
87 	auto variance        = properGeneralizedInverseGaussianVariance!double(2.0, 3.0, 4.0);
88 	assert(approxEqual(secondRowMoment - expectation^^2, variance));
89 }
90 
91 
92 /++
93 Mean of the generalized Hyperbolic distribution
94 +/
95 T generalizedHyperbolicMean(T)(T lambda, T beta, T eta, T omega)
96 {
97 	return beta * properGeneralizedInverseGaussianMean(lambda, eta, omega);
98 }
99 
100 ///
101 unittest
102 {
103 	auto expectation = generalizedHyperbolicMean!double(2.0, 5.0, 3.0, 4.0);
104 }
105 
106 
107 /++
108 Variance of the generalized Hyperbolic distribution
109 +/
110 T generalizedHyperbolicVariance(T)(T lambda, T beta, T eta, T omega)
111 {
112 	immutable r = besselKR(lambda, omega);
113 	return eta * r + (beta * eta) ^^ 2 * (r * (2 * (lambda + 1) / omega - r) + 1);
114 }
115 
116 ///
117 unittest
118 {
119 	auto variance = generalizedHyperbolicVariance!double(2.0, 5.0, 3.0, 4.0);
120 }