test/multivariate_gaussian.cc

Code
Comments
Other
Rev Date Author Line
4281 29 Jan 23 peter 1 // $Id$
4281 29 Jan 23 peter 2
4281 29 Jan 23 peter 3 /*
4281 29 Jan 23 peter 4   Copyright (C) 2023 Peter Johansson
4281 29 Jan 23 peter 5
4281 29 Jan 23 peter 6   This file is part of the yat library, https://dev.thep.lu.se/yat
4281 29 Jan 23 peter 7
4281 29 Jan 23 peter 8   The yat library is free software; you can redistribute it and/or
4281 29 Jan 23 peter 9   modify it under the terms of the GNU General Public License as
4281 29 Jan 23 peter 10   published by the Free Software Foundation; either version 3 of the
4281 29 Jan 23 peter 11   License, or (at your option) any later version.
4281 29 Jan 23 peter 12
4281 29 Jan 23 peter 13   The yat library is distributed in the hope that it will be useful,
4281 29 Jan 23 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
4281 29 Jan 23 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4281 29 Jan 23 peter 16   General Public License for more details.
4281 29 Jan 23 peter 17
4281 29 Jan 23 peter 18   You should have received a copy of the GNU General Public License
4281 29 Jan 23 peter 19   along with yat. If not, see <https://www.gnu.org/licenses/>.
4281 29 Jan 23 peter 20 */
4281 29 Jan 23 peter 21
4281 29 Jan 23 peter 22 #include <config.h>
4281 29 Jan 23 peter 23
4281 29 Jan 23 peter 24 #include "Suite.h"
4281 29 Jan 23 peter 25
4281 29 Jan 23 peter 26 #include "yat/random/random.h"
4281 29 Jan 23 peter 27
4281 29 Jan 23 peter 28 #include <yat/statistics/AveragerPair.h>
4281 29 Jan 23 peter 29
4281 29 Jan 23 peter 30 #include "yat/utility/CholeskyDecomposer.h"
4281 29 Jan 23 peter 31 #include "yat/utility/Matrix.h"
4281 29 Jan 23 peter 32 #include "yat/utility/Vector.h"
4281 29 Jan 23 peter 33
4281 29 Jan 23 peter 34 #include <vector>
4281 29 Jan 23 peter 35
4281 29 Jan 23 peter 36 using namespace theplu::yat;
4281 29 Jan 23 peter 37
4281 29 Jan 23 peter 38 int main(int argc, char* argv[])
4281 29 Jan 23 peter 39 {
4281 29 Jan 23 peter 40   test::Suite suite(argc, argv);
4281 29 Jan 23 peter 41
4281 29 Jan 23 peter 42   size_t n = 4;
4281 29 Jan 23 peter 43   utility::Vector m(n);
4281 29 Jan 23 peter 44   m(1) = 2.0;
4281 29 Jan 23 peter 45   m(2) = 1.0;
4281 29 Jan 23 peter 46   utility::Matrix Cov(n, n);
4281 29 Jan 23 peter 47   Cov(0,0) = 1.0;
4281 29 Jan 23 peter 48   Cov(1,1) = 1.0;
4281 29 Jan 23 peter 49   Cov(2,2) = 2.0;
4281 29 Jan 23 peter 50   Cov(3,3) = 10.0;
4281 29 Jan 23 peter 51   Cov(0,1) = Cov(1,0) = 0.8;
4281 29 Jan 23 peter 52
4281 29 Jan 23 peter 53   random::MultivariateGaussian gaussian(m, utility::CholeskyDecomposer(Cov));
4281 29 Jan 23 peter 54
4281 29 Jan 23 peter 55   std::vector<statistics::Averager> mean(n);
4281 29 Jan 23 peter 56   std::vector<std::vector<statistics::AveragerPair>>
4281 29 Jan 23 peter 57     averagers(4, std::vector<statistics::AveragerPair>(4));
4281 29 Jan 23 peter 58   for (size_t i=0; i<1000000; ++i) {
4281 29 Jan 23 peter 59     utility::Vector x(n);
4281 29 Jan 23 peter 60     gaussian(x);
4281 29 Jan 23 peter 61     for (size_t j=0; j<4; ++j)
4281 29 Jan 23 peter 62       for (size_t k=0; k<4; ++k)
4281 29 Jan 23 peter 63         averagers[j][k].add(x(j), x(k));
4281 29 Jan 23 peter 64   }
4281 29 Jan 23 peter 65
4281 29 Jan 23 peter 66   suite.out() << "mean:\n";
4281 29 Jan 23 peter 67   for (size_t i=0; i<n; ++i) {
4281 29 Jan 23 peter 68     const statistics::Averager a = averagers[i][0].x_averager();
4281 29 Jan 23 peter 69     suite.out() << i << "\t" << m(i) << "\t" << a.mean() << "\n";
4281 29 Jan 23 peter 70     if (!suite.equal_fix(a.mean(), m(i), 10 * a.standard_error()))
4281 29 Jan 23 peter 71       suite.add(false);
4281 29 Jan 23 peter 72   }
4281 29 Jan 23 peter 73
4281 29 Jan 23 peter 74   suite.out() << "covariance:\n";
4281 29 Jan 23 peter 75   for (size_t i=0; i<n; ++i) {
4281 29 Jan 23 peter 76     for (size_t j=i; j<n; ++j) {
4281 29 Jan 23 peter 77       suite.out() << i << "\t" << j << "\t"
4281 29 Jan 23 peter 78                   << Cov(i,j) << "\t"
4281 29 Jan 23 peter 79                   << averagers[i][j].covariance() << "\n";
4281 29 Jan 23 peter 80       // 1/n sum X^2/sigma^2 = 1/n Xi^2(n) i.e.
4281 29 Jan 23 peter 81       // mean: 1.0; variance: 2/n
4281 29 Jan 23 peter 82       // Cov = 1/n sum X^2 = sigma^2 * 1/n Xi^2
4281 29 Jan 23 peter 83       // mean: sigma^2; variance: 2/n * sigma^4
4281 29 Jan 23 peter 84       double margin = 10 * std::sqrt(Cov(i,i)*Cov(j,j)) *
4281 29 Jan 23 peter 85         std::sqrt(2.0/averagers[i][j].n());
4281 29 Jan 23 peter 86       if (!suite.equal_fix(averagers[i][j].covariance(),
4281 29 Jan 23 peter 87                            Cov(i,j), margin))
4281 29 Jan 23 peter 88         suite.add(false);
4281 29 Jan 23 peter 89     }
4281 29 Jan 23 peter 90   }
4281 29 Jan 23 peter 91
4281 29 Jan 23 peter 92   return suite.return_value();
4281 29 Jan 23 peter 93 }

Generated on Sat Oct 14 02:21:15 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)