1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // QUESO - a library to support the Quantification of Uncertainty
5 // for Estimation, Simulation and Optimization
6 //
7 // Copyright (C) 2008-2017 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24
25 #include <queso/BasicPdfsCXX11.h>
26
27 #ifdef QUESO_HAVE_CXX11
28
29 #include <cmath>
30
31 namespace QUESO {
32
BasicPdfsCXX11(int worldRank)33 BasicPdfsCXX11::BasicPdfsCXX11(int worldRank)
34 :
35 BasicPdfsBase(worldRank)
36 {
37 }
38
~BasicPdfsCXX11()39 BasicPdfsCXX11::~BasicPdfsCXX11()
40 {
41 }
42
43 double
betaPdfActualValue(double x,double alpha,double beta)44 BasicPdfsCXX11::betaPdfActualValue(double x, double alpha, double beta) const
45 {
46 // C++11 doesn't have a beta distribution, so we just implement one.
47
48 double ln_gamma_apb;
49 double ln_gamma_a;
50 double ln_gamma_b;
51
52 if (x < 0 || x > 1) { // Outside the support; return zero.
53 return 0.0;
54 }
55 else {
56 if (x == 0.0 || x == 1.0) { // On the boundary
57 if (alpha > 1.0 && beta > 1.0) { // Not blowing up
58 return 0.0;
59 }
60 else { // Might be blowing up, so let it.
61 ln_gamma_apb = std::lgamma(alpha + beta);
62 ln_gamma_a = std::lgamma(alpha);
63 ln_gamma_b = std::lgamma(beta);
64
65 return std::exp(ln_gamma_apb - ln_gamma_a - ln_gamma_b) *
66 std::pow(x, alpha - 1) * std::pow(1 - x, beta - 1);
67 }
68 }
69 else { // Compute beta pdf in log space.
70 ln_gamma_apb = std::lgamma(alpha + beta);
71 ln_gamma_a = std::lgamma(alpha);
72 ln_gamma_b = std::lgamma(beta);
73
74 return std::exp(ln_gamma_apb - ln_gamma_a - ln_gamma_b +
75 (alpha - 1) * std::log(x) + (beta - 1) * std::log1p(-x));
76 }
77 }
78 }
79
80 double
gammaPdfActualValue(double x,double a,double b)81 BasicPdfsCXX11::gammaPdfActualValue(double x, double a, double b) const
82 {
83 if (x < 0) { // Out of bounds.
84 return 0;
85 }
86 else {
87 if (x == 0.0) { // On the boundary.
88 if (a > 1.0) { // Not blowing up.
89 return 0.0;
90 }
91 else {
92 double ln_gamma_a = std::lgamma(a);
93 return std::exp(-ln_gamma_a - a * std::log(b)) * std::pow(x, a - 1) *
94 std::exp(- x / b);
95 }
96 }
97 else {
98 double ln_gamma_a = std::lgamma(a);
99 return std::exp(-ln_gamma_a - a * std::log(b) + (a - 1) * std::log(x) -
100 (x / b));
101 }
102 }
103 }
104
105 } // End namespace QUESO
106
107 #endif // QUESO_HAVE_CXX11
108