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 : // Bootstrap CI
10 :
11 : #include "BootstrapCalculators.h"
12 :
13 : namespace StochasticTools
14 : {
15 :
16 : MooseEnum
17 2236 : makeBootstrapCalculatorEnum()
18 : {
19 4472 : return MooseEnum("percentile=0 bca=1");
20 : }
21 :
22 : // PERCENTILE //////////////////////////////////////////////////////////////////////////////////////
23 : template <typename InType, typename OutType>
24 : std::vector<OutType>
25 3440 : Percentile<InType, OutType>::compute(const InType & data, const bool is_distributed)
26 : {
27 : // Bootstrap estimates
28 3440 : const std::vector<OutType> values = this->computeBootstrapEstimates(data, is_distributed);
29 :
30 : // Extract percentiles
31 : std::vector<OutType> output;
32 3440 : if (this->processor_id() == 0)
33 14730 : for (const Real & level : this->_levels)
34 : {
35 11994 : long unsigned int index = std::lrint(level * (this->_replicates - 1));
36 11994 : output.push_back(values[index]);
37 : }
38 :
39 3440 : return output;
40 3440 : }
41 :
42 : // BIASCORRECTEDACCELERATED ////////////////////////////////////////////////////////////////////////
43 : template <typename InType, typename OutType>
44 : std::vector<OutType>
45 101 : BiasCorrectedAccelerated<InType, OutType>::compute(const InType & data, const bool is_distributed)
46 : {
47 : // Bootstrap estimates
48 101 : const std::vector<OutType> values = this->computeBootstrapEstimates(data, is_distributed);
49 :
50 : // Compute bias-correction, Efron and Tibshirani (2003), Eq. 14.14, p. 186
51 101 : const OutType value = this->_calc.compute(data, is_distributed);
52 101 : const Real count = std::count_if(
53 : values.begin(),
54 : values.end(),
55 : [&value](Real v) { return v < value; }); // use Real for non-integer division below
56 101 : const Real bias = NormalDistribution::quantile(count / this->_replicates, 0, 1);
57 :
58 : // Compute Acceleration, Efron and Tibshirani (2003), Eq. 14.15, p. 186
59 101 : const OutType acc = data.empty() ? OutType() : acceleration(data, is_distributed);
60 :
61 : // Compute intervals, Efron and Tibshirani (2003), Eq. 14.10, p. 185
62 : std::vector<OutType> output;
63 1004 : for (const Real & level : this->_levels)
64 : {
65 903 : const Real z = NormalDistribution::quantile(level, 0, 1);
66 903 : const Real x = bias + (bias + (bias + z) / (1 - acc * (bias + z)));
67 903 : const Real alpha = NormalDistribution::cdf(x, 0, 1);
68 :
69 903 : long unsigned int index = std::lrint(alpha * (this->_replicates - 1));
70 903 : output.push_back(values[index]);
71 : }
72 101 : return output;
73 101 : }
74 :
75 : template <typename InType, typename OutType>
76 : OutType
77 101 : BiasCorrectedAccelerated<InType, OutType>::acceleration(const InType & data,
78 : const bool is_distributed)
79 : {
80 : const std::size_t local_size = data.size();
81 101 : std::vector<std::size_t> local_sizes = {local_size};
82 101 : if (is_distributed)
83 33 : this->_communicator.allgather(local_sizes);
84 101 : const std::size_t count = std::accumulate(local_sizes.begin(), local_sizes.end(), 0);
85 101 : const processor_id_type rank = is_distributed ? this->processor_id() : 0;
86 :
87 : // Jackknife statistics
88 101 : std::vector<OutType> theta_i(local_size);
89 :
90 : // Compute jackknife estimates, Ch. 11, Eq. 11.2, p. 141
91 268 : for (processor_id_type r = 0; r < local_sizes.size(); ++r)
92 4961 : for (std::size_t i = 0; i < local_sizes[r]; ++i)
93 : {
94 4794 : this->_calc.initializeCalculator();
95 2064238 : for (std::size_t il = 0; il < local_size; ++il)
96 2059444 : if (i != il || r != rank)
97 2055310 : this->_calc.updateCalculator(data[il]);
98 4794 : this->_calc.finalizeCalculator(is_distributed);
99 4794 : if (r == rank)
100 4134 : theta_i[i] = this->_calc.getValue();
101 : }
102 :
103 : // Compute jackknife sum, Ch. 11, Eq. 11.4, p. 141
104 101 : OutType theta_dot = std::accumulate(theta_i.begin(), theta_i.end(), OutType());
105 101 : if (is_distributed)
106 33 : this->_communicator.sum(theta_dot);
107 101 : theta_dot /= count;
108 :
109 : // Acceleration, Ch. 14, Eq. 14.15, p. 185
110 101 : std::vector<OutType> num_den(2);
111 4235 : for (const auto & jk : theta_i)
112 : {
113 4134 : num_den[0] += MathUtils::pow(theta_dot - jk, 3);
114 4134 : num_den[1] += MathUtils::pow(theta_dot - jk, 2);
115 : }
116 101 : if (is_distributed)
117 33 : this->_communicator.sum(num_den);
118 :
119 : mooseAssert(num_den[1] != OutType(), "The acceleration denomenator must not be zero.");
120 101 : return num_den[0] / (6 * std::pow(num_den[1], 3. / 2.));
121 101 : }
122 :
123 : template <typename InType, typename OutType>
124 : std::unique_ptr<BootstrapCalculator<InType, OutType>>
125 2667 : BootstrapCalculatorBuilder<InType, OutType>::build(
126 : const MooseEnum & item,
127 : const libMesh::ParallelObject & other,
128 : const std::vector<Real> & levels,
129 : unsigned int replicates,
130 : unsigned int seed,
131 : StochasticTools::Calculator<InType, OutType> & calc)
132 : {
133 : std::unique_ptr<BootstrapCalculator<InType, OutType>> ptr = nullptr;
134 2667 : if (item == "percentile")
135 2532 : ptr =
136 : std::make_unique<Percentile<InType, OutType>>(other, item, levels, replicates, seed, calc);
137 135 : else if (item == "bca")
138 135 : ptr = std::make_unique<BiasCorrectedAccelerated<InType, OutType>>(
139 : other, item, levels, replicates, seed, calc);
140 :
141 2667 : if (!ptr)
142 0 : ::mooseError("Failed to create Statistics::BootstrapCalculator object for ", item);
143 :
144 2667 : return ptr;
145 0 : }
146 :
147 : #define createBootstrapCalculators(InType, OutType) \
148 : template class Percentile<InType, OutType>; \
149 : template class BiasCorrectedAccelerated<InType, OutType>; \
150 : template struct BootstrapCalculatorBuilder<InType, OutType>
151 :
152 : createBootstrapCalculators(std::vector<Real>, Real);
153 : createBootstrapCalculators(std::vector<int>, Real);
154 :
155 : } // namespace
|
