Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "Beta.h"
11 : #include "math.h"
12 : #include "libmesh/utility.h"
13 :
14 : registerMooseObject("StochasticToolsApp", Beta);
15 :
16 : InputParameters
17 14 : Beta::validParams()
18 : {
19 14 : InputParameters params = Distribution::validParams();
20 14 : params.addClassDescription("Beta distribution");
21 28 : params.addRequiredRangeCheckedParam<Real>("alpha", "alpha > 0", "Shape parameter 1.");
22 28 : params.addRequiredRangeCheckedParam<Real>("beta", "beta > 0", "Shape parameter 2.");
23 14 : return params;
24 0 : }
25 :
26 7 : Beta::Beta(const InputParameters & parameters)
27 21 : : Distribution(parameters), _alpha(getParam<Real>("alpha")), _beta(getParam<Real>("beta"))
28 : {
29 7 : }
30 :
31 : Real
32 21 : Beta::pdf(const Real & x, const Real & alpha, const Real & beta)
33 : {
34 21 : if (x <= 0.0 || x > 1.0)
35 : return 0.0;
36 : else
37 21 : return std::pow(x, alpha - 1.0) * std::pow(1.0 - x, beta - 1.0) / betaFunction(alpha, beta);
38 : }
39 :
40 : Real
41 7 : Beta::cdf(const Real & x, const Real & alpha, const Real & beta)
42 : {
43 7 : if (x <= 0.0)
44 : return 0.0;
45 7 : else if (x >= 1.0)
46 : return 1.0;
47 : else
48 7 : return incompleteBeta(alpha, beta, x);
49 : }
50 :
51 : Real
52 7 : Beta::quantile(const Real & p, const Real & alpha, const Real & beta)
53 : {
54 7 : return incompleteBetaInv(alpha, beta, p);
55 : }
56 :
57 : Real
58 7 : Beta::pdf(const Real & x) const
59 : {
60 7 : return pdf(x, _alpha, _beta);
61 : }
62 :
63 : Real
64 7 : Beta::cdf(const Real & x) const
65 : {
66 7 : return cdf(x, _alpha, _beta);
67 : }
68 :
69 : Real
70 7 : Beta::quantile(const Real & p) const
71 : {
72 7 : return quantile(p, _alpha, _beta);
73 : }
74 :
75 : Real
76 1127 : Beta::betaFunction(const Real & a, const Real & b)
77 : {
78 1127 : return std::tgamma(a) * std::tgamma(b) / std::tgamma(a + b);
79 : }
80 :
81 : Real
82 1736 : Beta::incompleteBeta(const Real & a, const Real & b, const Real & x)
83 : {
84 1736 : if (x > ((a + 1.0) / (a + b + 2.0)))
85 672 : return 1.0 - incompleteBeta(b, a, 1.0 - x);
86 :
87 : const Real tol = 1e-14;
88 : const unsigned int max_iter = 1e3;
89 1064 : const Real coef = std::pow(x, a) * std::pow(1.0 - x, b) / a / betaFunction(a, b);
90 : Real fn = 1.0;
91 : Real cn = 1.0;
92 : Real dn = 0.0;
93 : Real num;
94 : Real m = 0.0;
95 12278 : for (unsigned int i = 0; i < max_iter; ++i)
96 : {
97 12278 : if (i == 0)
98 : num = 1.0;
99 11214 : else if (i % 2 == 0)
100 : {
101 5320 : m = i / 2.0;
102 5320 : num = m * (b - m) * x / (a + 2.0 * m) / (a + 2.0 * m - 1.0);
103 : }
104 : else
105 : {
106 5894 : m = (i - 1) / 2.0;
107 5894 : num = -(a + m) * (a + b + m) * x / (a + 2.0 * m) / (a + 2.0 * m + 1.0);
108 : }
109 :
110 12278 : dn = 1.0 / (1.0 + num * dn);
111 12278 : cn = 1.0 + num / cn;
112 12278 : const Real cd = cn * dn;
113 12278 : fn *= cd;
114 :
115 12278 : if (std::abs(1.0 - cd) < tol)
116 1064 : return coef * (fn - 1.0);
117 : }
118 :
119 : mooseAssert(false, "Could not compute incomplete beta function.");
120 0 : return coef * (fn - 1.0);
121 : }
122 :
123 : Real
124 42 : Beta::incompleteBetaInv(const Real & a, const Real & b, const Real & p)
125 : {
126 42 : if (a < b)
127 14 : return 1.0 - incompleteBetaInv(b, a, 1.0 - p);
128 :
129 : // Try with Newton method
130 28 : Real x = 0.5;
131 : const Real tol = 1e-14;
132 : const unsigned int max_iter = 1e3;
133 28 : const Real scale = betaFunction(a, b);
134 : Real f;
135 : Real df;
136 56 : for (unsigned int i = 0; i < max_iter; ++i)
137 : {
138 56 : f = incompleteBeta(a, b, x) - p;
139 56 : if (std::abs(f) < tol)
140 7 : return x;
141 49 : df = std::pow(x, a - 1.0) * std::pow(1.0 - x, b - 1.0) / scale;
142 49 : x -= f / df;
143 : // There's a divergence so try out bisection
144 49 : if (x < 0 || x > 1.0)
145 : break;
146 : }
147 :
148 : // Try with bisection method
149 : Real x1 = 0; // f(0) = -p
150 : Real x2 = 1; // f(1) = 1 - p
151 980 : for (unsigned int i = 0; i < max_iter; ++i)
152 : {
153 980 : x = (x2 + x1) / 2.0;
154 980 : f = incompleteBeta(a, b, x) - p;
155 980 : if (std::abs(f) < tol)
156 21 : return x;
157 959 : else if (f < 0)
158 504 : x1 = x;
159 : else
160 455 : x2 = x;
161 : }
162 :
163 : mooseAssert(false, "Could not find inverse of incomplete gamma function.");
164 0 : return x;
165 : }
|
