https://mooseframework.inl.gov
Loading...
Searching...
No Matches
TestBootstrapCalculators.C
Go to the documentation of this file.
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#include <vector>
10#include <cmath>
11
12#include "gtest/gtest.h"
13#include "Calculators.h"
14#include "BootstrapCalculators.h"
15#include "MooseRandom.h"
16#include "Normal.h"
17#include "libmesh/communicator.h"
18#include "libmesh/parallel_object.h"
19
20using namespace StochasticTools;
21
35class NormalSampler
36{
37public:
38 NormalSampler(Real mean, Real std, unsigned int seed) : _mean(mean), _std(std)
39 {
40 _generator.seed(seed);
41 }
42
43 Real sample() const { return _generator.randNormal(_mean, _std); }
44 std::vector<Real> sample(std::size_t n) const
45 {
46 std::vector<Real> data(n);
47 std::for_each(
48 data.begin(), data.end(), [&](Real & v) { v = _generator.randNormal(_mean, _std); });
49 return data;
50 }
51
52 Real meanConfidence(Real level, std::size_t n) const
53 {
54 const Real z = computeZ(level);
55 const Real mean_std = _std / std::sqrt(n);
56 return z * mean_std;
57 }
58
59 Real stdConfidence(Real level, std::size_t n) const
60 {
61 const Real z = computeZ(level);
62 const Real nr = n;
63 const Real gg = std::exp(std::lgamma(nr / 2.) - std::lgamma((nr - 1.) / 2.));
64 const Real std_std = _std * std::sqrt(1.0 - 2. / (nr - 1.) * gg * gg);
65 return z * std_std;
66 }
67
68private:
69 static Real computeZ(Real level)
70 {
71 const Real alpha = level < 0.5 ? level : 1. - level;
72 const Real z = -Normal::quantile(alpha / 2., 0, 1);
73 return level < 0.5 ? -z : z;
74 }
75
76 const Real _mean;
77 const Real _std;
78 MooseRandom _generator;
79};
80
81TEST(BootstrapCalculators, Percentile)
82{
83 // Sampling quantities
84 const Real mean_dist = 1993;
85 const Real std_dist = 27;
86 const std::size_t nsamp = 1000;
87
88 // CI quantities
89 const unsigned int replicates = 1e4;
90 const std::vector<Real> levels = {0.05, 0.1, 0.2, 0.8, 0.9, 0.95};
91
92 // Parallel object to give to calculators
93 Parallel::Communicator comm;
94 libMesh::ParallelObject po(comm);
95
96 // Construct mean and standard-deviation calculators
97 MultiMooseEnum calc("mean stddev", "mean stddev", true);
98 auto mean_calc = makeCalculator(calc[0], po);
99 auto std_calc = makeCalculator(calc[1], po);
100
101 // Construct bootstrap calculators
102 MooseEnum boot("percentile", "percentile", true);
103 auto mean_boot_calc = makeBootstrapCalculator(boot, po, levels, replicates, 2613, *mean_calc);
104 auto std_boot_calc = makeBootstrapCalculator(boot, po, levels, replicates, 2613, *std_calc);
105
106 // Construct data and run bootstrapping
107 NormalSampler sampler(mean_dist, std_dist, 1945);
108 const auto data = sampler.sample(nsamp);
109 const Real mean_samp = mean_calc->compute(data, false);
110 const Real std_samp = std_calc->compute(data, false);
111 const std::vector<Real> mean_ci = mean_boot_calc->compute(data, false);
112 const std::vector<Real> std_ci = std_boot_calc->compute(data, false);
113
114 // Compare with reference values
115 const Real tol = 5e-1;
116 for (const auto & l : index_range(levels))
117 {
118 const Real mean_ref = sampler.meanConfidence(levels[l], nsamp);
119 const Real std_ref = sampler.stdConfidence(levels[l], nsamp);
120 EXPECT_NEAR(mean_ci[l] - mean_samp, mean_ref, std::abs(mean_ref * tol));
121 EXPECT_NEAR(std_ci[l] - std_samp, std_ref, std::abs(std_ref * tol));
122 }
123}
124
125TEST(BootstrapCalculators, BiasCorrectedAccelerated)
126{
127 // Sampling quantities
128 const Real mean_dist = 1993;
129 const Real std_dist = 27;
130 const std::size_t nsamp = 1000;
131
132 // CI quantities
133 const unsigned int replicates = 1e4;
134 const std::vector<Real> levels = {0.05, 0.1, 0.2, 0.8, 0.9, 0.95};
135
136 // Parallel object to give to calculators
137 Parallel::Communicator comm;
138 libMesh::ParallelObject po(comm);
139
140 // Construct mean and standard-deviation calculators
141 MultiMooseEnum calc("mean stddev", "mean stddev", true);
142 auto mean_calc = makeCalculator(calc[0], po);
143 auto std_calc = makeCalculator(calc[1], po);
144
145 // Construct bootstrap calculators
146 MooseEnum boot("bca", "bca", true);
147 auto mean_boot_calc = makeBootstrapCalculator(boot, po, levels, replicates, 2613, *mean_calc);
148 auto std_boot_calc = makeBootstrapCalculator(boot, po, levels, replicates, 2613, *std_calc);
149
150 // Construct data and run bootstrapping
151 NormalSampler sampler(mean_dist, std_dist, 1945);
152 const auto data = sampler.sample(nsamp);
153 const Real mean_samp = mean_calc->compute(data, false);
154 const Real std_samp = std_calc->compute(data, false);
155 const std::vector<Real> mean_ci = mean_boot_calc->compute(data, false);
156 const std::vector<Real> std_ci = std_boot_calc->compute(data, false);
157
158 // Compare with reference values
159 const Real tol = 5e-1;
160 for (const auto & l : index_range(levels))
161 {
162 const Real mean_ref = sampler.meanConfidence(levels[l], nsamp);
163 const Real std_ref = sampler.stdConfidence(levels[l], nsamp);
164 EXPECT_NEAR(mean_ci[l] - mean_samp, mean_ref, std::abs(mean_ref * tol));
165 EXPECT_NEAR(std_ci[l] - std_samp, std_ref, std::abs(std_ref * tol));
166 }
167}
168
169TEST(BootstrapCalculators, Percentile_Vec)
170{
171 // Sampling quantities
172 const std::size_t nsamp = 1000;
173 const std::size_t nval = 26;
174 std::vector<NormalSampler> samplers;
175 for (const auto & k : make_range(nval))
176 samplers.emplace_back(/*mean = */ 1993 + 42 * k, /*std = */ 27 + 7 * k, 1945);
177
178 // CI quantities
179 const unsigned int replicates = 1e4;
180 const std::vector<Real> levels = {0.05, 0.1, 0.2, 0.8, 0.9, 0.95};
181
182 // Parallel object to give to calculators
183 Parallel::Communicator comm;
184 libMesh::ParallelObject po(comm);
185
186 // Construct mean and standard-deviation calculators
187 MultiMooseEnum calc("mean stddev", "mean stddev", true);
188 auto mean_calc = makeCalculator<std::vector<std::vector<Real>>, std::vector<Real>>(calc[0], po);
189 auto std_calc = makeCalculator<std::vector<std::vector<Real>>, std::vector<Real>>(calc[1], po);
190
191 // Construct bootstrap calculators
192 MooseEnum boot("percentile", "percentile", true);
193 auto mean_boot_calc = makeBootstrapCalculator<std::vector<std::vector<Real>>, std::vector<Real>>(
194 boot, po, levels, replicates, 2613, *mean_calc);
195 auto std_boot_calc = makeBootstrapCalculator<std::vector<std::vector<Real>>, std::vector<Real>>(
196 boot, po, levels, replicates, 2613, *std_calc);
197
198 // Construct data and run bootstrapping
199 std::vector<std::vector<Real>> data(nsamp);
200 for (auto & dt : data)
201 for (const auto & samp : samplers)
202 dt.push_back(samp.sample());
203 const std::vector<Real> mean_samp = mean_calc->compute(data, false);
204 const std::vector<Real> std_samp = std_calc->compute(data, false);
205 const std::vector<std::vector<Real>> mean_ci = mean_boot_calc->compute(data, false);
206 const std::vector<std::vector<Real>> std_ci = std_boot_calc->compute(data, false);
207
208 // Compare with reference values
209 const Real tol = 5e-1;
210 for (const auto & l : index_range(levels))
211 for (const auto & k : make_range(nval))
212 {
213 const Real mean_ref = samplers[k].meanConfidence(levels[l], nsamp);
214 const Real std_ref = samplers[k].stdConfidence(levels[l], nsamp);
215 EXPECT_NEAR(mean_ci[l][k] - mean_samp[k], mean_ref, std::abs(mean_ref * tol));
216 EXPECT_NEAR(std_ci[l][k] - std_samp[k], std_ref, std::abs(std_ref * tol));
217 }
218}
219
220TEST(BootstrapCalculators, BiasCorrectedAccelerated_Vec)
221{
222 // Sampling quantities
223 const std::size_t nsamp = 1000;
224 const std::size_t nval = 26;
225 std::vector<NormalSampler> samplers;
226 for (const auto & k : make_range(nval))
227 samplers.emplace_back(/*mean = */ 1993 + 42 * k, /*std = */ 27 + 7 * k, 1945);
228
229 // CI quantities
230 const unsigned int replicates = 1e4;
231 const std::vector<Real> levels = {0.05, 0.1, 0.2, 0.8, 0.9, 0.95};
232
233 // Parallel object to give to calculators
234 Parallel::Communicator comm;
235 libMesh::ParallelObject po(comm);
236
237 // Construct mean and standard-deviation calculators
238 MultiMooseEnum calc("mean stddev", "mean stddev", true);
239 auto mean_calc = makeCalculator<std::vector<std::vector<Real>>, std::vector<Real>>(calc[0], po);
240 auto std_calc = makeCalculator<std::vector<std::vector<Real>>, std::vector<Real>>(calc[1], po);
241
242 // Construct bootstrap calculators
243 MooseEnum boot("bca", "bca", true);
244 auto mean_boot_calc = makeBootstrapCalculator<std::vector<std::vector<Real>>, std::vector<Real>>(
245 boot, po, levels, replicates, 2613, *mean_calc);
246 auto std_boot_calc = makeBootstrapCalculator<std::vector<std::vector<Real>>, std::vector<Real>>(
247 boot, po, levels, replicates, 2613, *std_calc);
248
249 // Construct data and run bootstrapping
250 std::vector<std::vector<Real>> data(nsamp);
251 for (auto & dt : data)
252 for (const auto & samp : samplers)
253 dt.push_back(samp.sample());
254 const std::vector<Real> mean_samp = mean_calc->compute(data, false);
255 const std::vector<Real> std_samp = std_calc->compute(data, false);
256 const std::vector<std::vector<Real>> mean_ci = mean_boot_calc->compute(data, false);
257 const std::vector<std::vector<Real>> std_ci = std_boot_calc->compute(data, false);
258
259 // Compare with reference values
260 const Real tol = 5e-1;
261 for (const auto & l : index_range(levels))
262 for (const auto & k : make_range(nval))
263 {
264 const Real mean_ref = samplers[k].meanConfidence(levels[l], nsamp);
265 const Real std_ref = samplers[k].stdConfidence(levels[l], nsamp);
266 EXPECT_NEAR(mean_ci[l][k] - mean_samp[k], mean_ref, std::abs(mean_ref * tol));
267 EXPECT_NEAR(std_ci[l][k] - std_samp[k], std_ref, std::abs(std_ref * tol));
268 }
269}
tol
const double tol
Definition EquilibriumConstantFitTest.C:16
v
const double v
Definition GasLiquidMassTransferTest.C:23
TEST
TEST(BootstrapCalculators, Percentile)
Definition TestBootstrapCalculators.C:81
MooseRandom::seed
void seed(std::size_t i, unsigned int seed)
MooseRandom::randNormal
Real randNormal(std::size_t i, Real mean, Real sigma)
NormalSampler
These tests are meant to use the bootstrap calculators and test against analytical confidence interva...
Definition TestBootstrapCalculators.C:36
NormalSampler::NormalSampler
NormalSampler(Real mean, Real std, unsigned int seed)
Definition TestBootstrapCalculators.C:38
NormalSampler::computeZ
static Real computeZ(Real level)
Definition TestBootstrapCalculators.C:69
NormalSampler::meanConfidence
Real meanConfidence(Real level, std::size_t n) const
Definition TestBootstrapCalculators.C:52
NormalSampler::sample
std::vector< Real > sample(std::size_t n) const
Definition TestBootstrapCalculators.C:44
NormalSampler::stdConfidence
Real stdConfidence(Real level, std::size_t n) const
Definition TestBootstrapCalculators.C:59
Normal::quantile
virtual Real quantile(const Real &p) const override
Definition Normal.C:80
StochasticTools
Enum for batch type in stochastic tools MultiApp.
Definition StochasticToolsTypes.h:14
StochasticTools::makeCalculator
std::unique_ptr< Calculator< InType, OutType > > makeCalculator(const MooseEnumItem &item, const libMesh::ParallelObject &other)
Definition Calculators.h:335
StochasticTools::makeBootstrapCalculator
std::unique_ptr< BootstrapCalculator< InType, OutType > > makeBootstrapCalculator(const MooseEnum &, const libMesh::ParallelObject &, const std::vector< Real > &, unsigned int, unsigned int, StochasticTools::Calculator< InType, OutType > &)
Definition BootstrapCalculators.h:183
Generated on Fri Sep 11 2026 20:37:02 for https://mooseframework.inl.gov by doxygen 1.9.8