Belle II Software light-2406-ragdoll
CalcMeanCov.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8#include <framework/utilities/TestHelpers.h>
9#include <framework/utilities/CalcMeanCov.h>
10#include <TRandom.h>
11#include <TMath.h>
12
13#include <gtest/gtest.h>
14
15using namespace std;
16using namespace Belle2;
17
18namespace {
21 TEST(CalcMeanCov, Normal)
22 {
23 constexpr int N = 10000;
24 std::vector<double> values(N, 0);
25 CalcMeanCov<> meancov;
26
27 for (double gmean : { -5.0, 0.0, 1.0, 3.5, 150.}) {
28 for (double gsigma : {0.1, 1.0, 4.2, 25.}) {
29 double sum1(0), sum2(0), sum3(0);
30 for (int i = 0; i < N; ++i) {
31 const double x = gRandom->Gaus(gmean, gsigma);
32 sum1 += x;
33 values[i] = x;
34 meancov.add(x);
35 }
36 const double mean = sum1 / N;
37 for (auto x : values) {
38 const double d = x - mean;
39 sum2 += d * d;
40 sum3 += d;
41 }
42 const double variance = (sum2 - sum3 * sum3 / N) / N;
43 const double sigma = std::sqrt(variance);
44
45 EXPECT_FLOAT_EQ(N, meancov.getEntries());
46 EXPECT_FLOAT_EQ(mean, meancov.getMean());
47 EXPECT_FLOAT_EQ(sigma, meancov.getStddev());
48 EXPECT_FLOAT_EQ(sum1, meancov.getSum());
49 meancov.clear();
50 }
51 }
52 }
53
58 TEST(CalcMeanCov, ManyUniform)
59 {
60 CalcMeanCov<3> meancov;
61 constexpr int N = 10001;
62 //Let's use a districte, uniform distribution because sum,
63 //mean and variance are very easy to calculate
64 //Result should not depend on the order of entries and a constant offset
65 //should only influence mean and sum, not the variance/standard deviation
66 for (int i = 0; i < N; ++i) {
67 meancov.add(i + 1, i + 13.5, N - i);
68 }
69 const double sum1 = (N * N + N) / 2.;
70 const double sum2 = sum1 + 12.5 * N;
71 const double mean1 = (N + 1.0) / 2;
72 const double mean2 = (13.5 + N + 12.5) / 2;
73 const double varia = (N * N - 1.) / 12;
74 const double sigma = std::sqrt(varia);
75
76 EXPECT_DOUBLE_EQ(N, meancov.getEntries());
77 EXPECT_DOUBLE_EQ(mean1, meancov.getMean<0>());
78 EXPECT_DOUBLE_EQ(mean2, meancov.getMean<1>());
79 EXPECT_DOUBLE_EQ(mean1, meancov.getMean<2>());
80 EXPECT_DOUBLE_EQ(sigma, meancov.getStddev<0>());
81 EXPECT_DOUBLE_EQ(sigma, meancov.getStddev<1>());
82 EXPECT_DOUBLE_EQ(sigma, meancov.getStddev<2>());
83 EXPECT_DOUBLE_EQ(sum1, meancov.getSum<0>());
84 EXPECT_DOUBLE_EQ(sum2, meancov.getSum<1>());
85 EXPECT_DOUBLE_EQ(sum1, meancov.getSum<2>());
86 //Check that correlation is correct
87 EXPECT_DOUBLE_EQ(+1, (meancov.getCorrelation<0, 1>()));
88 EXPECT_DOUBLE_EQ(-1, (meancov.getCorrelation<0, 2>()));
89 EXPECT_DOUBLE_EQ(-1, (meancov.getCorrelation<1, 2>()));
90 //Check that the covariance matrix is symmetric
91 EXPECT_EQ((meancov.getCovariance<1, 0>()), (meancov.getCovariance<0, 1>()));
92 EXPECT_EQ((meancov.getCovariance<2, 0>()), (meancov.getCovariance<0, 2>()));
93 EXPECT_EQ((meancov.getCovariance<2, 1>()), (meancov.getCovariance<1, 2>()));
94 }
95
100 TEST(CalcMeanCov, Binomial)
101 {
102 CalcMeanCov<> meancov;
103 for (double n : {10, 13, 25}) {
104 for (double p : {0.1, 0.5, 0.75}) {
105 for (int k = 0; k <= n; ++k) {
106 const double b = TMath::Binomial(n, k) * std::pow(p, k) * std::pow(1 - p, n - k);
107 meancov.addWeighted(b, k);
108 }
109 EXPECT_DOUBLE_EQ(1.0, meancov.getEntries()) << "n: " << n << ", p: " << p;
110 EXPECT_DOUBLE_EQ(n * p, meancov.getMean()) << "n: " << n << ", p: " << p;
111 EXPECT_DOUBLE_EQ(n * p * (1 - p), meancov.getVariance()) << "n: " << n << ", p: " << p;
112 meancov.clear();
113 }
114 }
115 }
116
118 TEST(CalcMeanCov, TemplateVsArray)
119 {
120 CalcMeanCov<3> templated;
121 CalcMeanCov<3> arrayed;
122 double values[3];
123 for (int i = 0; i < 10000; ++i) {
124 const double weight = gRandom->Uniform(0, 10);
125 values[0] = gRandom->Gaus(0, 1);
126 values[1] = gRandom->Gaus(values[0], 1);
127 values[2] = gRandom->Uniform(-5, 5);
128 templated.addWeighted(weight, values[0], values[1], values[2]);
129 arrayed.addWeightedArray(weight, values);
130
131 ASSERT_EQ(templated.getEntries(), arrayed.getEntries());
132 ASSERT_EQ(templated.getSum<0>(), arrayed.getSum(0));
133 ASSERT_EQ(templated.getSum<1>(), arrayed.getSum(1));
134 ASSERT_EQ(templated.getSum<2>(), arrayed.getSum(2));
135 ASSERT_EQ(templated.getMean<0>(), arrayed.getMean(0));
136 ASSERT_EQ(templated.getMean<1>(), arrayed.getMean(1));
137 ASSERT_EQ(templated.getMean<2>(), arrayed.getMean(2));
138 ASSERT_EQ(templated.getStddev<0>(), arrayed.getStddev(0));
139 ASSERT_EQ(templated.getStddev<1>(), arrayed.getStddev(1));
140 ASSERT_EQ(templated.getStddev<2>(), arrayed.getStddev(2));
141 ASSERT_EQ((templated.getCovariance<0, 1>()), arrayed.getCovariance(0, 1));
142 ASSERT_EQ((templated.getCovariance<0, 2>()), arrayed.getCovariance(0, 2));
143 ASSERT_EQ((templated.getCovariance<1, 2>()), arrayed.getCovariance(1, 2));
144 }
145 }
146
148 TEST(CalcMeanCov, AddingTwoSets)
149 {
150 constexpr int N = 10000;
151 CalcMeanCov<> meancov; //adding one by one
152 CalcMeanCov<> sum;
153 CalcMeanCov<> otherhalf; //added to sum
154
155 for (double gmean : { -5.0, 0.0, 1.0, 3.5, 150.}) {
156 for (double gsigma : {0.1, 1.0, 4.2, 25.}) {
157 for (int i = 0; i < N; ++i) {
158 const double x = gRandom->Gaus(gmean, gsigma);
159 meancov.add(x);
160 if (i % 2 == 0)
161 sum.add(x);
162 else
163 otherhalf.add(x);
164 }
165 sum.add(otherhalf);
166
167 EXPECT_FLOAT_EQ(sum.getEntries(), meancov.getEntries());
168 EXPECT_FLOAT_EQ(sum.getMean(), meancov.getMean());
169 EXPECT_FLOAT_EQ(sum.getStddev(), meancov.getStddev());
170 EXPECT_FLOAT_EQ(sum.getSum(), meancov.getSum());
171 meancov.clear();
172 sum.clear();
173 otherhalf.clear();
174 }
175 }
176 }
177} // namespace
Class to calculate mean and and covariance between a number of parameters on running data without sto...
Definition: CalcMeanCov.h:35
void addWeighted(value_type weight, T... values)
Update mean and covariance by adding a new, weighted entry.
Definition: CalcMeanCov.h:59
value_type getVariance(int i) const
Return the variance for paramter i.
Definition: CalcMeanCov.h:137
value_type getCorrelation(int i, int j) const
Return the correlation coefficient between parameters i and j.
Definition: CalcMeanCov.h:130
value_type getStddev(int i) const
Return the standard deviation for parameter i.
Definition: CalcMeanCov.h:139
value_type getMean(int i) const
Return the mean for parameter i.
Definition: CalcMeanCov.h:122
value_type getCovariance(int i, int j) const
Return the covariance between parameters i and j.
Definition: CalcMeanCov.h:124
void addWeightedArray(value_type weight, value_type *values)
Update mean and covarianced by adding a new weighted entry.
Definition: CalcMeanCov.h:100
value_type getSum(int i) const
Return the weighted sum values for parameter i.
Definition: CalcMeanCov.h:141
void add(T... values)
Update mean and covariance by adding a new entry.
Definition: CalcMeanCov.h:71
void clear()
Clear all values.
Definition: CalcMeanCov.h:45
value_type getEntries() const
Return the number of entries.
Definition: CalcMeanCov.h:120
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.