1 /++
2 Probability density functions
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.pdf;
12 
13 import core.stdc.tgmath;
14 
15 import std.traits;
16 import std.numeric;
17 import std.math : isNormal, isFinite, M_2_SQRTPI, SQRT2, PI, approxEqual;
18 import std.typecons;
19 import std.algorithm;
20 import std.array;
21 
22 import atmosphere.params;
23 import atmosphere.utilities;
24 import atmosphere.math: besselK;
25 
26 
27 /++
28 Normal PDF
29 +/
30 struct NormalSPDF(T)
31 {
32 	private T c, mu, sigma2;
33 
34 	/++
35 	Params:
36 		mu = location
37 		sigma2 = sigma^2
38 	+/
39 	this(T mu, T sigma2)
40 	in {
41 		assert(sigma2 > 0);
42 		assert(mu.isFinite);
43 	}
44 	body {
45 		c = 1 / sqrt(2*PI*sigma2);
46 		this.mu = mu;
47 		this.sigma2 = sigma2;
48 	}
49 
50 	///
51 	T opCall(T x) const
52 	{
53 		import std.math : exp;
54 		return c * exp((x-mu)^^2 / (-2*sigma2));
55 	}
56 }
57 
58 ///
59 unittest
60 {
61 	auto pdf = NormalSPDF!double(1.9, 2.3);
62 }
63 
64 /++
65 Gamma PDF
66 +/
67 struct GammaSPDF(T)
68 	if(isFloatingPoint!T)
69 {
70 	import std.mathspecial : gamma;
71 	private T shapem1, scale, c;
72 
73 	/++
74 	Params:
75 		shape = gamma shape parameter
76 		scale = gamma scale parameter
77 	+/
78 	this(T shape, T scale)
79 	in {
80 		assert(shape.isNormal);
81 		assert(shape > 0);
82 		assert(scale.isNormal);
83 		assert(scale > 0);
84 	}
85 	body {
86 		this.c = 1 / (gamma(shape) * scale);
87 		assert(c.isNormal);
88 		this.shapem1 = shape - 1;
89 		this.scale = scale;
90 	}
91 
92 	///
93 	T opCall(T x) const
94 	{
95 		if(x < 0)
96 			return 0;
97 		x /= scale;
98 		return c 
99 			* pow(x, shapem1) 
100 			* exp(-x);
101 	}
102 }
103 
104 ///
105 unittest 
106 {
107 	auto pdf = GammaSPDF!double(3, 2);
108 	auto x = pdf(0.1);
109 	assert(x.isNormal);
110 
111 	import scid.calculus : integrate;
112 	auto result = pdf.integrate(-double.infinity, double.infinity);
113 	assert(fabs(result.value - 1) < result.error);
114 }
115 
116 
117 /++
118 Inverse-gamma PDF
119 +/
120 struct InverseGammaSPDF(T)
121 	if(isFloatingPoint!T)
122 {
123 	import std.mathspecial : gamma;
124 	private T shapem1, scale, c;
125 
126 	/++
127 	Params:
128 		shape = gamma shape parameter
129 		scale = gamma scale parameter
130 	+/
131 	this(T shape, T scale)
132 	in {
133 		assert(shape.isNormal);
134 		assert(shape > 0);
135 		assert(scale.isNormal);
136 		assert(scale > 0);
137 	}
138 	body {
139 		this.c = 1 / (gamma(shape) * scale);
140 		assert(c.isNormal);
141 		this.shapem1 = shape + 1;
142 		this.scale = scale;
143 	}
144 
145 	///
146 	T opCall(T x) const
147 	{
148 		if(x < 0)
149 			return 0;
150 		x /= scale;
151 		return c 
152 			* pow(x, -shapem1) 
153 			* exp(-1/x);
154 	}
155 }
156 
157 ///
158 unittest 
159 {
160 	auto pdf = InverseGammaSPDF!double(3, 4);
161 	auto x = pdf(0.1);
162 	assert(x.isNormal);
163 
164 	import scid.calculus : integrate;
165 	auto result = pdf.integrate(-double.infinity, double.infinity);
166 	assert(fabs(result.value - 1) < result.error);
167 }
168 
169 
170 /++
171 Generalized gamma PDF
172 +/
173 struct GeneralizedGammaSPDF(T)
174 	if(isFloatingPoint!T)
175 {
176 	import std.mathspecial : gamma;
177 	private T shapem1, power, scale, c, e;
178 
179 	/++
180 	Params:
181 		shape = shape parameter
182 		power = power parameter
183 		scale = scale parameter
184 	+/
185 	this(T shape, T power, T scale)
186 	in {
187 		assert(shape.isNormal);
188 		assert(shape > 0);
189 		assert(power.isFinite);
190 		assert(scale.isNormal);
191 		assert(scale > 0);
192 	}
193 	body {
194 		this.power = power;
195 		this.scale = scale;
196 		this.c = fabs(power) / gamma(shape) / scale;
197 		this.e = power * shape - 1;
198 	}
199 
200 	///
201 	T opCall(T x) const
202 	{
203 		if(x < 0)
204 			return 0;
205 		x /= scale;
206 		return c
207 			* pow(x, e)
208 			* exp(-pow(x, power));
209 	}
210 }
211 
212 ///
213 unittest 
214 {
215 	auto pdf = GeneralizedGammaSPDF!double(3, 2, 0.5);
216 	auto x = pdf(0.1);
217 	assert(x.isNormal);
218 
219 	import scid.calculus : integrate;
220 	auto result = pdf.integrate(-double.infinity, double.infinity);
221 	assert(fabs(result.value - 1) < result.error);
222 }
223 
224 
225 /++
226 Inverse Gaussian PDF
227 
228 See_Also: [distribution.params](distribution/params.html)
229 +/
230 struct InverseGaussianSPDF(T)
231 	if(isFloatingPoint!T)
232 {
233 	private T omega, chi, psi;
234 
235 	///
236 	this(T chi, T psi)
237 	in {
238 		assert(chi.isNormal);
239 		assert(chi > 0);
240 		assert(psi.isNormal);
241 		assert(psi > 0);
242 	}
243 	body {
244 		this.chi = chi;
245 		this.psi = psi;
246 		this.omega = GIGChiPsi!T(chi, psi).omega;
247 	}
248 
249 	///
250 	T opCall(T x) const
251 	{
252 		return x < 0 ? 0 : sqrt(chi / (2*PI*x^^3)) * exp(omega - (chi / x + psi * x) / 2);
253 	}
254 }
255 
256 ///
257 unittest 
258 {
259 	auto pdf = InverseGaussianSPDF!double(3, 2);
260 	auto x = pdf(0.1);
261 	assert(x.isNormal);
262 
263 	import scid.calculus : integrate;
264 	auto result = pdf.integrate(-double.infinity, double.infinity);
265 	assert(fabs(result.value - 1) < result.error);
266 }
267 
268 
269 /++
270 Proper generalized inverse Gaussian PDF
271 
272 See_Also: [distribution.params](distribution/params.html)
273 +/
274 struct ProperGeneralizedInverseGaussianSPDF(T)
275 	if(isFloatingPoint!T)
276 {
277 	private T omega, eta, c, lambdam1;
278 
279 	///
280 	this(T lambda, T eta, T omega)
281 	in {
282 		assert(lambda.isFinite);
283 		assert(eta.isNormal);
284 		assert(eta > 0);
285 		assert(omega.isNormal);
286 		assert(omega > 0);
287 
288 	}
289 	body {
290 		this.c = 2 * eta * besselK!(Flag!"ExponentiallyScaled".yes)(lambda, omega);
291 		this.eta = eta;
292 		this.omega = omega;
293 		this.lambdam1 = lambda - 1;
294 	}
295 
296 	///
297 	T opCall(T x) const
298 	{
299 		if(x <= 0)
300 			return 0;
301 		x /= eta;
302 		return x < 0 ? 0 : pow(x, lambdam1) * exp(omega * (1 - (1/x + x) / 2)) / c;
303 	}
304 }
305 
306 ///
307 unittest 
308 {
309 	auto pdf = ProperGeneralizedInverseGaussianSPDF!double(4, 3, 2);
310 	auto x = pdf(0.1);
311 	assert(x.isNormal);
312 
313 	import scid.calculus : integrate;
314 	auto result = pdf.integrate(-double.infinity, double.infinity);
315 	assert(fabs(result.value - 1) < result.error);
316 }
317 
318 
319 /++
320 Variance gamma (gamma mixture of normals) PDF
321 
322 See_Also: [distribution.params](distribution/params.html)
323 +/
324 struct VarianceGammaSPDF(T)
325 	if(isFloatingPoint!T)
326 {
327 	import std.mathspecial : gamma;
328 	private T lambdamh, alpha, beta, mu, c;
329 
330 	/++
331 	Params:
332 		lambda = 
333 		alpha = 
334 		beta = 
335 		delta = 0
336 		mu = location
337 	+/
338 	this(T lambda, T alpha, T beta, T mu)
339 	in {
340 		assert(lambda.isNormal);
341 		assert(lambda > 0);
342 		assert(alpha.isNormal);
343 		assert(alpha > 0);
344 		assert(beta.isFinite);
345 		assert(beta > -alpha);
346 		assert(beta < +alpha);
347 		assert(mu.isFinite);
348 	}
349 	body {
350 		immutable params = GHypAlphaDelta!T(alpha, beta, 0);
351 		this.c = M_2_SQRTPI / SQRT2 * pow(params.psi / 2, lambda) / gamma(lambda);
352 		assert(c.isNormal);
353 		this.lambdamh = lambda - 0.5f;
354 		this.alpha = alpha;
355 		this.beta = beta;
356 		this.mu = mu;
357 	}
358 
359 	///
360 	T opCall(T x) const
361 	{
362 		immutable y = x - mu;
363 		immutable z = fabs(y);
364 		immutable a = alpha * z;
365 		return c 
366 			* pow(z / alpha, lambdamh)
367 			* besselK!(Flag!"ExponentiallyScaled".yes)(lambdamh, a)
368 			* exp(beta * y - a);
369 	}
370 }
371 
372 ///
373 unittest 
374 {
375 	auto pdf = VarianceGammaSPDF!double(1.1, 1.1, 1.0, 1.1);
376 	auto x = pdf(0.1);
377 	assert(x.isNormal);
378 
379 	import scid.calculus : integrate;
380 	auto result = pdf.integrate(-double.infinity, double.infinity);
381 	assert(fabs(result.value - 1) < result.error);
382 }
383 
384 
385 /++
386 Hyperbolic asymmetric T (inverse gamma mixture of normals) PDF
387 
388 See_Also: [distribution.params](distribution/params.html)
389 +/
390 struct HyperbolicAsymmetricTSPDF(T)
391 	if(isFloatingPoint!T)
392 {
393 	import std.mathspecial : gamma;
394 	private T lambdamh, beta, delta, mu, c;
395 
396 	/++
397 	Params:
398 		lambda = 
399 		alpha = `|beta|`
400 		beta = 
401 		delta = 
402 		mu = location
403 	+/
404 	this(T lambda, T beta, T delta, T mu)
405 	in {
406 		assert(lambda.isNormal);
407 		assert(lambda < 0);
408 		assert(beta.isFinite);
409 		assert(delta.isNormal);
410 		assert(delta > 0);
411 		assert(mu.isFinite);
412 	}
413 	body {
414 		immutable params = GHypAlphaDelta!T(fabs(beta), beta, delta);
415 		this.c = M_2_SQRTPI / SQRT2 * pow(params.chi / 2, -lambda) / gamma(-lambda);
416 		assert(c.isNormal);
417 		this.lambdamh = lambda - 0.5f;
418 		this.beta = beta;
419 		this.delta = delta;
420 		this.mu = mu;
421 	}
422 
423 	///
424 	T opCall(T x) const
425 	{
426 		immutable y = x - mu;
427 		immutable z = hypot(delta, y);
428 		immutable alpha = fabs(beta);
429 		immutable a = z * alpha;
430 		return c
431 			* pow(z / alpha, lambdamh)
432 			* besselK!(Flag!"ExponentiallyScaled".yes)(lambdamh, a)
433 			* exp(beta * y - a);
434 	}
435 }
436 
437 ///
438 unittest 
439 {
440 	auto pdf = HyperbolicAsymmetricTSPDF!double(-1.1, 1.1, 1.1, 1.1);
441 	auto x = pdf(0.1);
442 	import std.conv;
443 	assert(x.isNormal, text(x));
444 }
445 
446 
447 /++
448 Normal-inverse Gaussian (inverse Gaussian mixture of normals) PDF
449 
450 See_Also: [distribution.params](distribution/params.html)
451 +/
452 struct NormalInverseGaussianSPDF(T)
453 	if(isFloatingPoint!T)
454 {
455 	private T alpha, beta, delta, mu, c;
456 
457 	/++
458 	Params:
459 		lambda = -1/2
460 		alpha = 
461 		beta = 
462 		delta = 
463 		mu = location
464 	+/
465 	this(T alpha, T beta, T delta, T mu)
466 	in {
467 		assert(alpha.isNormal);
468 		assert(alpha > 0);
469 		assert(beta.isFinite);
470 		assert(beta > -alpha);
471 		assert(beta < +alpha);
472 		assert(delta.isNormal);
473 		assert(delta > 0);
474 		assert(mu.isFinite);
475 	}
476 	body {
477 		immutable params = GHypAlphaDelta!T(alpha, beta, delta);
478 		this.c = exp(params.omega) * alpha * delta / PI;
479 		assert(c.isNormal);
480 		this.alpha = alpha;
481 		this.beta = beta;
482 		this.delta = delta;
483 		this.mu = mu;
484 	}
485 
486 	///
487 	T opCall(T x) const
488 	{
489 		immutable y = x - mu;
490 		immutable z = hypot(delta, y);
491 		immutable a = z * alpha;
492 		return c
493 			/ z
494 			* besselK!(Flag!"ExponentiallyScaled".yes)(T(1), a)
495 			* exp(beta * y - a);
496 	}
497 }
498 
499 ///
500 unittest 
501 {
502 	auto pdf = NormalInverseGaussianSPDF!double(1.1, 0.8, 1.1, 1.1);
503 	auto x = pdf(0.1);
504 	assert(x.isNormal);
505 
506 	import scid.calculus : integrate;
507 	auto result = pdf.integrate(-double.infinity, double.infinity);
508 	assert(fabs(result.value - 1) < result.error);
509 }
510 
511 
512 /++
513 Proper generalized hyperbolic (generalized inverse Gaussian mixture of normals) PDF
514 
515 See_Also: [distribution.params](distribution/params.html)
516 +/
517 struct ProperGeneralizedHyperbolicSPDF(T)
518 	if(isFloatingPoint!T)
519 {
520 	private T lambdamh, alpha, beta, delta, mu, omega, c;
521 
522 	/++
523 	Params:
524 		lambda =
525 		alpha = 
526 		beta = 
527 		delta = 
528 		mu = location
529 	+/
530 	this(T lambda, T alpha, T beta, T delta, T mu)
531 	in {
532 		assert(lambda.isFinite);
533 		assert(alpha.isNormal);
534 		assert(alpha > 0);
535 		assert(beta.isFinite);
536 		assert(beta > -alpha);
537 		assert(beta < +alpha);
538 		assert(delta.isNormal);
539 		assert(delta > 0);
540 		assert(mu.isFinite);
541 	}
542 	body {
543 		immutable params = GHypAlphaDelta!T(alpha, beta, delta);
544 		this.omega = params.omega;
545 		this.c = M_2_SQRTPI / (SQRT2 * 2)
546 			* pow(params.eta, -lambda)
547 			/ besselK!(Flag!"ExponentiallyScaled".yes)(lambda, omega);
548 		assert(c.isNormal);
549 		this.lambdamh = lambda - 0.5f;
550 		this.alpha = alpha;
551 		this.beta = beta;
552 		this.delta = delta;
553 		this.mu = mu;
554 	}
555 
556 	///
557 	T opCall(T x) const
558 	{
559 		immutable y = x - mu;
560 		immutable z = hypot(delta, y);
561 		immutable a = z * alpha;
562 		return c
563 			* pow(z / alpha, lambdamh)
564 			* besselK!(Flag!"ExponentiallyScaled".yes)(lambdamh, a)
565 			* exp(beta * y - (a - omega));
566 	}
567 }
568 
569 ///
570 unittest 
571 {
572 	auto pdf = ProperGeneralizedHyperbolicSPDF!double(1.1, 1.1, 0.8, 1.1, 1.1);
573 	auto x = pdf(0.1);
574 	assert(x.isNormal);
575 
576 	import scid.calculus : integrate;
577 	auto result = pdf.integrate(-double.infinity, double.infinity);
578 	assert(fabs(result.value - 1) < result.error);
579 }
580 
581 
582 /++
583 Probability density function interface
584 +/
585 interface PDF(T)
586 {
587 	/**
588 	Call operator
589 	*/
590 	T opCall(T x);
591 }
592 
593 
594 ///
595 unittest
596 {
597 	import std.traits, std.mathspecial;
598 
599 	class NormalPDF : PDF!real
600 	{
601 		real opCall(real x)
602 		{
603 			// 1/sqrt(2 PI)
604 			enum c = 0.398942280401432677939946L;
605 			return c * exp(-0.5f * x * x);
606 		}
607 	}
608 
609 	auto pdf = new NormalPDF;
610 	auto x = pdf(0.1);
611 	assert(x.isNormal);
612 
613 	import scid.calculus : integrate;
614 	auto result = pdf.integrate(-double.infinity, double.infinity);
615 	assert(fabs(result.value - 1) < result.error);
616 }
617 
618 ///
619 alias toPDF = convertTo!PDF;
620 
621 ///
622 unittest
623 {
624 	PDF!double pdf = GammaSPDF!double(1, 3).toPDF;
625 }
626 
627 /++
628 Variance-mean mixture of normals. $(RED Unstable) algorithm.
629 +/
630 abstract class NormalVarianceMeanMixturePDF(T) : PDF!T
631 	if(isFloatingPoint!T)
632 {
633 	private PDF!T pdf;
634 	private T beta, mu;
635 	private T epsRel, epsAbs;
636 	private const(T)[] subdivisions;
637 
638 	/++
639     Params:
640 		pdf     = The PDF to _integrate.
641 		beta = NVMM scale
642 		mu = NVMM location
643 		subdivisions     = TODO.
644 		epsRel  = (optional) The requested relative accuracy.
645 		epsAbs  = (optional) The requested fabsolute accuracy.
646 	See_also: [struct Result](https://github.com/kyllingstad/scid/blob/a9f3916526e4bf9a4da35d14a969e1abfa17a496/source/scid/types.d)
647 	+/
648 	this(PDF!T pdf, T beta, T mu, const(T)[] subdivisions = null, T epsRel = 1e-6, T epsAbs = 0)
649 	in {
650 		assert(subdivisions.all!isFinite);
651 		assert(subdivisions.all!(s => s > 0));
652 		assert(subdivisions.isSorted);
653 		assert(subdivisions.findAdjacent.empty);
654 	}
655 	body {
656 		this.pdf = pdf;
657 		this.beta = beta;
658 		this.mu = mu;
659 		this.epsRel = epsRel;
660 		this.epsAbs = epsAbs;
661 		this.subdivisions = subdivisions;
662 	}
663 
664 	T opCall(T x)
665 	{
666 		import scid.calculus : integrate;
667 		T f(T z) {
668 			return exp(-0.5f * (x - (mu + beta * z)) ^^ 2 / z) / sqrt(2*PI*z) * pdf(z);
669 		}
670 		T sum = 0;
671 		T a = 0;
672 		foreach(s; subdivisions)
673 		{
674 			sum += integrate(&f, a, s, epsRel, epsAbs);
675 			a = s;
676 		}
677 		sum += integrate(&f, a, T.infinity);
678 		return sum;
679 	}
680 }
681 
682 ///
683 unittest
684 {
685 	class MyGHypPDF(T) : NormalVarianceMeanMixturePDF!T ///$(RED Unstable) algorithm.
686 	{
687 		import atmosphere.moment;
688 		this(T lambda, GHypEtaOmega!T params, T mu)
689 		{
690 			with(params)
691 			{
692 				auto pgig = ProperGeneralizedInverseGaussianSPDF!T(lambda, eta, omega);
693 				auto e = mu + properGeneralizedInverseGaussianMean(lambda, eta, omega);
694 				super(pgig.toPDF, params.beta, mu, [e]);				
695 			}
696 		}
697 	}
698 
699 	import atmosphere.params;
700 	immutable double lambda = 2;
701 	immutable double mu = 0.3;
702 	immutable params = GHypEtaOmega!double(2, 3, 4);
703 	auto pghyp = new GeneralizedHyperbolicPDF!double(lambda, params.alpha, params.beta, params.delta, mu);
704 	auto pnvmm = new MyGHypPDF!double(lambda, params, mu);
705 	foreach(i; 0..5)
706 		assert(approxEqual(pghyp(i), pnvmm(i)));
707 }
708 
709 
710 /++
711 Generalized variance-gamma (generalized gamma mixture of normals) PDF. $(RED Unstable) algorithm.
712 +/
713 final class GeneralizedVarianceGammaPDF(T) : NormalVarianceMeanMixturePDF!T
714 	if(isFloatingPoint!T)
715 {
716 	/++
717 	Params:
718 		shape = shape parameter (generalized gamma)
719 		power = power parameter (generalized gamma)
720 		scale = scale parameter (generalized gamma)
721 		beta = NVMM scale
722 		mu = NVMM location
723 	+/
724 	this(T shape, T power, T scale, T beta, T mu)
725 	in {
726 		assert(shape.isNormal);
727 		assert(shape > 0);
728 		assert(power.isFinite);
729 		assert(scale.isNormal);
730 		assert(scale > 0);
731 		assert(beta.isFinite);
732 		assert(mu.isFinite);
733 	}
734 	body {
735 		import atmosphere.moment : generalizedGammaMean;
736 		auto pdf  = GeneralizedGammaSPDF!double(shape, power, scale);
737 		auto e = generalizedGammaMean(shape, power, scale);
738 		assert(e > 0);
739 		assert(e.isFinite);
740 		super(pdf.toPDF, beta, mu, [e]);
741 	}
742 }
743 
744 ///
745 unittest 
746 {
747 	auto pdf = new GeneralizedVarianceGammaPDF!double(1.1, 1.1, 0.9, 1.1, 1.1);
748 	auto x = pdf(0.1);
749 	assert(x.isNormal);
750 
751 	import scid.calculus : integrate;
752 	auto result = pdf.integrate(-double.infinity, double.infinity);
753 	assert(fabs(result.value - 1) < result.error);
754 }
755 
756 
757 /++
758 Generalized inverse Gaussian PDF
759 
760 See_Also: [distribution.params](distribution/params.html)
761 +/
762 final class GeneralizedInverseGaussianPDF(T) : PDF!T
763 	if(isFloatingPoint!T)
764 {
765 	private PDF!T pdf;
766 
767 	///
768 	this(T lambda, T chi, T psi)
769 	in {
770 		assert(lambda.isFinite);
771 		assert(chi.isNormal);
772 		assert(chi >= 0);
773 		assert(psi.isNormal);
774 		assert(psi >= 0);
775 	}
776 	body {
777 		immutable params = GIGChiPsi!T(chi, psi);
778 		if (chi < T.min_normal)
779 			this.pdf = GammaSPDF!T(lambda, 2 / psi).toPDF;
780 		else if (psi < T.min_normal)
781 			this.pdf = InverseGammaSPDF!T(-lambda, chi / 2).toPDF;
782 		else if (lambda == -0.5f)
783 			this.pdf = InverseGaussianSPDF!T(params.eta, chi).toPDF;
784 		else
785 			this.pdf = ProperGeneralizedInverseGaussianSPDF!T(lambda, params.eta, params.omega).toPDF;
786 	}
787 
788 	T opCall(T x)
789 	{
790 		return pdf(x);
791 	}
792 }
793 
794 ///
795 unittest 
796 {
797 	auto pdf = new GeneralizedInverseGaussianPDF!double(3, 2, 1);
798 	auto x = pdf(0.1);
799 	assert(x.isNormal);
800 
801 	import scid.calculus : integrate;
802 	auto result = pdf.integrate(-double.infinity, double.infinity);
803 	assert(fabs(result.value - 1) < result.error);
804 }
805 
806 
807 /++
808 Generalized hyperbolic (generalized inverse Gaussian mixture of normals) PDF
809 
810 See_Also: [distribution.params](distribution/params.html)
811 +/
812 final class GeneralizedHyperbolicPDF(T) : PDF!T
813 	if(isFloatingPoint!T)
814 {
815 	private PDF!T pdf;
816 
817 	/++
818 	+/
819 	this(T lambda, T alpha, T beta, T delta, T mu)
820 	in {
821 		assert(lambda.isFinite);
822 		assert(alpha.isNormal || alpha == 0);
823 		assert(alpha >= 0);
824 		assert(beta.isFinite);
825 		assert(beta >= -alpha);
826 		assert(beta <= +alpha);
827 		assert(delta.isNormal || delta == 0);
828 		assert(delta >= 0);
829 		assert(mu.isFinite);
830 	}
831 	body {
832 		if (delta == 0)
833 			this.pdf = VarianceGammaSPDF!T(lambda, alpha, beta, mu).toPDF;
834 		else if (alpha == fabs(beta))
835 			this.pdf = HyperbolicAsymmetricTSPDF!T(lambda, beta, delta, mu).toPDF;
836 		else if (lambda == -0.5f)
837 			this.pdf = NormalInverseGaussianSPDF!T(alpha, beta, delta, mu).toPDF;
838 		else
839 			this.pdf = ProperGeneralizedHyperbolicSPDF!T(lambda, alpha, beta, delta, mu).toPDF;
840 	}
841 
842 	T opCall(T x)
843 	{
844 		return pdf(x);
845 	}
846 }
847 
848 ///
849 unittest 
850 {
851 	auto pdf = new GeneralizedHyperbolicPDF!double(1.1, 1.1, 0.9, 1.1, 1.1);
852 	auto x = pdf(0.1);
853 	assert(x.isNormal);
854 
855 	import scid.calculus : integrate;
856 	auto result = pdf.integrate(-double.infinity, double.infinity);
857 	assert(fabs(result.value - 1) < result.error);
858 }