atmosphere.estimate.gamma source code

1 /**
2 Authors: [Ilya Yaroshenko](http://9il.github.io)
3 
4 Copyright: © 2014-2015 [Ilya Yaroshenko](http://9il.github.io)
5 
6 License: MIT
7 */
8 module atmosphere.estimate.gamma;
9 
10 import core.stdc.tgmath;
11 
12 import std.traits;
13 import std.typecons;
14 
15 import atmosphere.statistic: GammaStatistic;
16 
17 /++
18 Estimates parameters of the gamma distribution.
19 +/
20 Tuple!(T, "shape", T, "scale")
21 gammaEstimate(T)(in T[] sample)
22 	if(isFloatingPoint!T)
23 {
24 	return gammaEstimate(GammaStatistic!T(sample));
25 }
26 
27 ///
28 unittest {
29 	import atmosphere.estimate.generalized_gamma;
30 	immutable sample = [1.0, 0.5, 0.75];
31 	immutable p0 = gammaEstimate(sample);
32 	immutable p1 = generalizedGammaFixedPowerEstimate(1.0, sample);
33 	assert(p0.shape == p1.shape);
34 	assert(p0.scale == p1.scale);
35 }
36 
37 
38 ///
39 unittest
40 {
41 	import std.range;
42 	import std.random;
43 	import atmosphere.random;
44 	import atmosphere.likelihood.gamma;
45 	auto length = 1000;
46 	auto shape = 2.0, scale = 3.0;
47 	auto rng = Random(1234);
48 	auto sample = GammaSRNG!double(rng, shape, scale).take(length).array;
49 	auto weights = iota(1.0, length + 1.0).array;
50 	auto params = gammaEstimate!double(sample, weights);
51 	auto lh0 = gammaLikelihood(shape, scale, sample, weights);
52 	auto lh1 = gammaLikelihood(params.shape, params.scale, sample, weights);
53 	assert(lh0 <= lh1);
54 }
55 
56 ///ditto
57 Tuple!(T, "shape", T, "scale")
58 gammaEstimate(T)(in T[] sample, in T[] weights)
59 	if(isFloatingPoint!T)
60 {
61 	return gammaEstimate(GammaStatistic!T(sample, weights));
62 }
63 
64 ///
65 unittest {
66 	import atmosphere.estimate.generalized_gamma;
67 	immutable sample = [1.0, 0.5, 0.75];
68 	immutable weights = [1.0, 4, 3];
69 	immutable p0 = gammaEstimate(sample, weights);
70 	immutable p1 = generalizedGammaFixedPowerEstimate(1.0, sample, weights);
71 	assert(p0.shape == p1.shape);
72 	assert(p0.scale == p1.scale);
73 }
74 
75 ///ditto
76 Tuple!(T, "shape", T, "scale")
77 gammaEstimate(T)(in GammaStatistic!T stat)
78 	if(isFloatingPoint!T)
79 {
80 	import atmosphere.math: logmdigammaInverse;
81 	with(stat)
82 	{
83 		immutable T shape = logmdigammaInverse(log(mean) - meanl);
84 		immutable T scale = mean / shape;
85 		return typeof(return)(shape, scale);
86 	}
87 }