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