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 34 : Beta::validParams()
18 : {
19 34 : InputParameters params = Distribution::validParams();
20 34 : params.addClassDescription("Beta distribution");
21 68 : params.addRequiredRangeCheckedParam<Real>("alpha", "alpha > 0", "Shape parameter 1.");
22 68 : params.addRequiredRangeCheckedParam<Real>("beta", "beta > 0", "Shape parameter 2.");
23 34 : return params;
24 0 : }
25 :
26 17 : Beta::Beta(const InputParameters & parameters)
27 51 : : Distribution(parameters), _alpha(getParam<Real>("alpha")), _beta(getParam<Real>("beta"))
28 : {
29 17 : }
30 :
31 : Real
32 51 : Beta::pdf(const Real & x, const Real & alpha, const Real & beta)
33 : {
34 51 : if (x <= 0.0 || x > 1.0)
35 : return 0.0;
36 : else
37 51 : return std::pow(x, alpha - 1.0) * std::pow(1.0 - x, beta - 1.0) / betaFunction(alpha, beta);
38 : }
39 :
40 : Real
41 17 : Beta::cdf(const Real & x, const Real & alpha, const Real & beta)
42 : {
43 17 : if (x <= 0.0)
44 : return 0.0;
45 17 : else if (x >= 1.0)
46 : return 1.0;
47 : else
48 17 : return incompleteBeta(alpha, beta, x);
49 : }
50 :
51 : Real
52 17 : Beta::quantile(const Real & p, const Real & alpha, const Real & beta)
53 : {
54 17 : return incompleteBetaInv(alpha, beta, p);
55 : }
56 :
57 : Real
58 17 : Beta::pdf(const Real & x) const
59 : {
60 17 : return pdf(x, _alpha, _beta);
61 : }
62 :
63 : Real
64 17 : Beta::cdf(const Real & x) const
65 : {
66 17 : return cdf(x, _alpha, _beta);
67 : }
68 :
69 : Real
70 17 : Beta::quantile(const Real & p) const
71 : {
72 17 : return quantile(p, _alpha, _beta);
73 : }
74 :
75 : Real
76 2737 : Beta::betaFunction(const Real & a, const Real & b)
77 : {
78 2737 : return std::tgamma(a) * std::tgamma(b) / std::tgamma(a + b);
79 : }
80 :
81 : Real
82 4216 : Beta::incompleteBeta(const Real & a, const Real & b, const Real & x)
83 : {
84 4216 : if (x > ((a + 1.0) / (a + b + 2.0)))
85 1632 : return 1.0 - incompleteBeta(b, a, 1.0 - x);
86 :
87 : const Real tol = 1e-14;
88 : const unsigned int max_iter = 1e3;
89 2584 : 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 29818 : for (unsigned int i = 0; i < max_iter; ++i)
96 : {
97 29818 : if (i == 0)
98 : num = 1.0;
99 27234 : else if (i % 2 == 0)
100 : {
101 12920 : m = i / 2.0;
102 12920 : num = m * (b - m) * x / (a + 2.0 * m) / (a + 2.0 * m - 1.0);
103 : }
104 : else
105 : {
106 14314 : m = (i - 1) / 2.0;
107 14314 : num = -(a + m) * (a + b + m) * x / (a + 2.0 * m) / (a + 2.0 * m + 1.0);
108 : }
109 :
110 29818 : dn = 1.0 / (1.0 + num * dn);
111 29818 : cn = 1.0 + num / cn;
112 29818 : const Real cd = cn * dn;
113 29818 : fn *= cd;
114 :
115 29818 : if (std::abs(1.0 - cd) < tol)
116 2584 : 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 102 : Beta::incompleteBetaInv(const Real & a, const Real & b, const Real & p)
125 : {
126 102 : if (a < b)
127 34 : return 1.0 - incompleteBetaInv(b, a, 1.0 - p);
128 :
129 : // Try with Newton method
130 68 : Real x = 0.5;
131 : const Real tol = 1e-14;
132 : const unsigned int max_iter = 1e3;
133 68 : const Real scale = betaFunction(a, b);
134 : Real f;
135 : Real df;
136 136 : for (unsigned int i = 0; i < max_iter; ++i)
137 : {
138 136 : f = incompleteBeta(a, b, x) - p;
139 136 : if (std::abs(f) < tol)
140 17 : return x;
141 119 : df = std::pow(x, a - 1.0) * std::pow(1.0 - x, b - 1.0) / scale;
142 119 : x -= f / df;
143 : // There's a divergence so try out bisection
144 119 : 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 2380 : for (unsigned int i = 0; i < max_iter; ++i)
152 : {
153 2380 : x = (x2 + x1) / 2.0;
154 2380 : f = incompleteBeta(a, b, x) - p;
155 2380 : if (std::abs(f) < tol)
156 51 : return x;
157 2329 : else if (f < 0)
158 1224 : x1 = x;
159 : else
160 1105 : x2 = x;
161 : }
162 :
163 : mooseAssert(false, "Could not find inverse of incomplete gamma function.");
164 0 : return x;
165 : }
|
