-
Notifications
You must be signed in to change notification settings - Fork 7
/
StatCalculate.hpp
177 lines (145 loc) · 6.99 KB
/
StatCalculate.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
// Flow Vector Correction Framework
//
// Copyright (C) 2020 Lukas Kreis Ilya Selyuzhenkov
// Contact: l.kreis@gsi.de; ilya.selyuzhenkov@gmail.com
// For a full list of contributors please see docs/Credits
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifndef QNTOOLS_STATBIN_H_
#define QNTOOLS_STATBIN_H_
#include <cmath>
#include <vector>
#include "Rtypes.h"
#include "Stat.hpp"
#include "StatCollect.hpp"
#include "Statistics.hpp"
namespace Qn {
class StatCalculate : public Stat {
public:
static StatCalculate PreferObservable(const StatCalculate &lhs, const StatCalculate &rhs);
/// Default constructor.
StatCalculate() = default;
/// Construct StatCalculate from Stats.
explicit StatCalculate(StatCollect &stats) {
mean_ = stats.GetStatistics().Mean();
variance_ = stats.GetStatistics().Variance();
sum_weight_ = stats.GetStatistics().SumWeights();
sum_weight2_ = stats.GetStatistics().SumWeights2();
const auto& samples = stats.GetBootStrap();
sample_means_ = samples.GetMeans();
sample_weights_ = samples.GetWeights();
if (stats.GetWeightType() == Stat::WeightType::OBSERVABLE) {
SetWeightType(WeightType::OBSERVABLE);
} else {
SetWeightType(WeightType::REFERENCE);
}
}
/// Construct StatCalculate from two Stats, prefering OBSERVABLE or REFERENCE.
StatCalculate(const StatCalculate &lhs, const StatCalculate &rhs) :
StatCalculate(PreferObservable(lhs, rhs)) {}
virtual ~StatCalculate();
/// Returns the sample Variance from error propagation
[[nodiscard]] double VarianceFromPropagation() const { return variance_; }
/// Returns variance of the sample mean from error propagation
[[nodiscard]] double VarianceOfMeanFromPropagation() const {
if (sum_weight_ > 0.) {
return VarianceFromPropagation() * NEffectiveInverse();
} else {
return 0.;
}
}
/// Returns variance of the sample mean from bootstrapping using the variance statistic.
[[nodiscard]] double VarianceOfMeanFromBootstrap() const {
Statistics stats;
for (std::size_t i = 0; i < sample_means_.size(); ++i) {
stats.Fill(sample_means_[i], sample_weights_[i]);
}
return stats.Variance();
}
/// Returns standard error of the sample mean from bootstrapping using the percentile method.
[[nodiscard]] double StdDevOfMeanFromBootstrapPercentile() const {
auto means = sample_means_;
const double alpha = 0.3173;
std::sort(means.begin(), means.end());
auto lower = means.at(std::floor(means.size() * alpha / 2.));
auto upper = means.at(std::ceil(means.size() * (1 - alpha / 2.)));
return (upper - lower) / 2;
}
/// Returns standard error from error propagation
[[nodiscard]] double StdDevFromPropagation() const {
return std::sqrt(VarianceFromPropagation());
}
/// Returns the standard error of the mean from error propagation
[[nodiscard]] double StdDevOfMeanFromPropagation() const {
return std::sqrt(VarianceOfMeanFromPropagation());
}
/// Returns the standard error of the mean from bootstrapping using the variance statistic.
[[nodiscard]] double StdDevOfMeanFromBootstrapVariance() const {
return std::sqrt(VarianceOfMeanFromBootstrap());
}
/// Returns the sample mean.
[[nodiscard]] double Mean() const { return mean_; }
/// Returns the sample variance.
[[nodiscard]] double Variance() const { return variance_; }
/// Returns the sum of weights.
[[nodiscard]] double SumWeights() const { return sum_weight_; }
/// Returns the sum of weights squared.
[[nodiscard]] double SumWeights2() const { return sum_weight2_; }
/// Retunrs the effective number of events.
/// See For more information:
/// https://en.wikipedia.org/w/index.php?title=Effective_sample_size&oldid=976981114#Weighted_samples
[[nodiscard]] double NEffectiveInverse() const { return SumWeights2() / (SumWeights() * SumWeights()); }
/// Returns standard error of the mean using the chosen method.
[[nodiscard]] double StandardErrorOfMean() const {
if (GetErrorType() == ErrorType::PROPAGATION) return StdDevOfMeanFromPropagation();
else return StdDevOfMeanFromBootstrapVariance();
}
/// Merge using pooled statistics
/// For merging of pooled variance see:
/// https://en.wikipedia.org/w/index.php?title=Pooled_variance&oldid=953883248#Sample-based_statistics
friend StatCalculate Merge(const StatCalculate &, const StatCalculate &);
friend StatCalculate MergeBins(const StatCalculate &, const StatCalculate &);
/// Implements algebra of random variables
/// For algebra of random variables used here see:
/// https://en.wikipedia.org/w/index.php?title=Algebra_of_random_variables&oldid=945632545
friend StatCalculate operator+(const StatCalculate &, const StatCalculate &);
friend StatCalculate operator-(const StatCalculate &, const StatCalculate &);
friend StatCalculate operator*(const StatCalculate &, const StatCalculate &);
friend StatCalculate operator/(const StatCalculate &, const StatCalculate &);
friend StatCalculate operator/(const StatCalculate &, double);
friend StatCalculate operator*(const StatCalculate &, double);
friend StatCalculate operator*(double, const StatCalculate &);
friend StatCalculate Pow(const StatCalculate &, double);
friend StatCalculate Sqrt(const StatCalculate &);
private:
std::vector<double> sample_means_; /// means of bootstrap samples
std::vector<double> sample_weights_; /// weights of bootstrap samples
double mean_ = 0.; /// mean of the distribution
double variance_ = 0.; /// variance of the distribution
double sum_weight_ = 0.; /// sum of weights
double sum_weight2_ = 0.; /// sum of {weights}^{2}
/// \cond CLASSIMP
ClassDef(StatCalculate, 3);
/// \endcond
};
StatCalculate Merge(const StatCalculate &lhs, const StatCalculate &rhs);
StatCalculate operator+(const StatCalculate &lhs, const StatCalculate &rhs);
StatCalculate operator-(const StatCalculate &lhs, const StatCalculate &rhs);
StatCalculate operator*(const StatCalculate &lhs, const StatCalculate &rhs);
StatCalculate operator/(const StatCalculate &num, const StatCalculate &den);
StatCalculate operator*(const StatCalculate &operand, double scale);
StatCalculate operator*(double operand, const StatCalculate &rhs);
StatCalculate operator/(const StatCalculate &operand, double scale);
StatCalculate Pow(const StatCalculate &base, double exp);
StatCalculate Sqrt(const StatCalculate &operand);
}
#endif // QNTOOLS_INCLUDE_BASE_QNTOOLS_STATBIN_H_