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