atmosphere.likelihood.generalized_inverse

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.likelihood.generalized_inverse_gaussian;
9 
10 import std.traits;
11 import std.typecons;
12 import std.math;
13 
14 import atmosphere.statistic: GeneralizedInverseGaussinStatistic;
15 
16 /++
17 Normalized log-likelihood function of the generalized inverse Gaussian distribution.
18 See_Also: `distribution.params.GIGEtaOmega`
19 +/
20 T properGeneralizedInverseGaussianLikelihood(T)(in T lambda, in T eta, in T omega, in T[] sample)
21 	if(isFloatingPoint!T)
22 {
23 	return properGeneralizedInverseGaussianLikelihood(lambda, eta, omega, GeneralizedInverseGaussinStatistic!T(sample));
24 }
25 
26 ///ditto
27 T properGeneralizedInverseGaussianLikelihood(T)(in T lambda, in T eta, in T omega, in T[] sample, in T[] weights)
28 	if(isFloatingPoint!T)
29 {
30 	return properGeneralizedInverseGaussianLikelihood(lambda, eta, omega, GeneralizedInverseGaussinStatistic!T(sample, weights));
31 }
32 
33 ///
34 unittest {
35 	immutable l = properGeneralizedInverseGaussianLikelihood!double(1,2,3,[1,2,2]);
36 	immutable m = properGeneralizedInverseGaussianLikelihood!double(1,2,3,[1,2],[2,4]);
37 	assert(l == m);
38 }
39 
40 ///ditto
41 T properGeneralizedInverseGaussianLikelihood(T)(in T lambda, in T eta, in T omega, in GeneralizedInverseGaussinStatistic!T stat)
42 	if(isFloatingPoint!T)
43 in {
44 	assert(omega >= T.min_normal);
45 }
46 body {
47 	version(LDC)
48 	{
49 		import ldc.intrinsics: log = llvm_log;
50 	}
51 	import std.math : LN2;
52 	import atmosphere.math: logBesselK;
53 	with(stat) return
54 		- cast(T)LN2 - lambda * log(eta) - logBesselK(lambda, omega)
55 		+ (lambda - 1) * stat.meanl
56 		- omega * (mean / eta + meani * eta) / 2;
57 }