atmosphere.estimate.generalized

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.generalized_gamma;
9 
10 import core.stdc.tgmath;
11 
12 import std.traits;
13 import std.typecons;
14 
15 import atmosphere.statistic: GeneralizedGammaFixedPowerStatistic;
16 
17 
18 /++
19 Estimates parameters of the generalized gamma distribution.
20 +/
21 Tuple!(T, "shape", T, "scale")
22 generalizedGammaFixedPowerEstimate(T)(in T power, in T[] sample)
23 	if(isFloatingPoint!T)
24 {
25 	return generalizedGammaFixedPowerEstimate(power, GeneralizedGammaFixedPowerStatistic!T(power, sample));
26 }
27 
28 ///
29 unittest
30 {
31 	import std.range;
32 	import std.random;
33 	import atmosphere.random;
34 	import atmosphere.likelihood;
35 	size_t length = 1000;
36 	auto shape = 2.0, power = 1.4, scale = 2.3;
37 	auto rng = Random(1234);
38 	auto sample = GeneralizedGammaSRNG!double(rng, shape, power, scale).take(length).array;
39 	auto params = generalizedGammaFixedPowerEstimate!double(power, sample);
40 	auto lh0 = generalizedGammaLikelihood(shape, power, scale, sample);
41 	auto lh1 = generalizedGammaLikelihood(params.shape, power, params.scale, sample);
42 	assert(lh0 <= lh1);
43 }
44 
45 
46 ///ditto
47 Tuple!(T, "shape", T, "scale")
48 generalizedGammaFixedPowerEstimate(T)(in T power, in T[] sample, in T[] weights)
49 	if(isFloatingPoint!T)
50 {
51 	return generalizedGammaFixedPowerEstimate(power, GeneralizedGammaFixedPowerStatistic!T(power, sample, weights));
52 }
53 
54 ///
55 unittest
56 {
57 	import std.range;
58 	import std.random;
59 	import atmosphere.random;
60 	import atmosphere.likelihood;
61 	size_t length = 1000;
62 	auto shape = 2.0, power = 1.4, scale = 2.3;
63 	auto rng = Random(1234);
64 	auto sample = GeneralizedGammaSRNG!double(rng, shape, power, scale).take(length).array;
65 	auto weights = iota(1.0, length + 1.0).array;
66 	auto params = generalizedGammaFixedPowerEstimate!double(power, sample, weights);
67 	auto lh0 = generalizedGammaLikelihood(shape, power, scale, sample, weights);
68 	auto lh1 = generalizedGammaLikelihood(params.shape, power, params.scale, sample, weights);
69 	assert(lh0 <= lh1);
70 }
71 
72 
73 ///ditto
74 Tuple!(T, "shape", T, "scale")
75 generalizedGammaFixedPowerEstimate(T)(in T power, in GeneralizedGammaFixedPowerStatistic!T stat)
76 	if(isFloatingPoint!T)
77 {
78 	import atmosphere.math: logmdigammaInverse;
79 	with(stat)
80 	{
81 		immutable T shape = logmdigammaInverse(log(meanp) - power * meanl);
82 		immutable T scale = pow(meanp / shape, 1 / power);
83 		return typeof(return)(shape, scale);		
84 	}
85 }