test/percentile_confidence_interval.cc

Code
Comments
Other
Rev Date Author Line
4035 23 Jan 21 peter 1 // $Id$
4035 23 Jan 21 peter 2
4035 23 Jan 21 peter 3 /*
4035 23 Jan 21 peter 4   Copyright (C) 2021 Peter Johansson
4035 23 Jan 21 peter 5
4035 23 Jan 21 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
4035 23 Jan 21 peter 7
4035 23 Jan 21 peter 8   The yat library is free software; you can redistribute it and/or
4035 23 Jan 21 peter 9   modify it under the terms of the GNU General Public License as
4035 23 Jan 21 peter 10   published by the Free Software Foundation; either version 3 of the
4035 23 Jan 21 peter 11   License, or (at your option) any later version.
4035 23 Jan 21 peter 12
4035 23 Jan 21 peter 13   The yat library is distributed in the hope that it will be useful,
4035 23 Jan 21 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
4035 23 Jan 21 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4035 23 Jan 21 peter 16   General Public License for more details.
4035 23 Jan 21 peter 17
4035 23 Jan 21 peter 18   You should have received a copy of the GNU General Public License
4035 23 Jan 21 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
4035 23 Jan 21 peter 20 */
4035 23 Jan 21 peter 21
4035 23 Jan 21 peter 22 #include <config.h>
4035 23 Jan 21 peter 23
4035 23 Jan 21 peter 24 #include "Suite.h"
4035 23 Jan 21 peter 25
4035 23 Jan 21 peter 26 #include <yat/statistics/Averager.h>
4035 23 Jan 21 peter 27 #include <yat/statistics/PercentileConfidenceInterval.h>
4035 23 Jan 21 peter 28 #include <yat/statistics/utility.h>
4035 23 Jan 21 peter 29 #include <yat/random/copy_k_of_n.h>
4035 23 Jan 21 peter 30
4036 24 Jan 21 peter 31 #include <gsl/gsl_cdf.h>
4036 24 Jan 21 peter 32 #include <gsl/gsl_randist.h>
4036 24 Jan 21 peter 33
4035 23 Jan 21 peter 34 #include <set>
4035 23 Jan 21 peter 35
4035 23 Jan 21 peter 36 using namespace theplu::yat;
4035 23 Jan 21 peter 37
4035 23 Jan 21 peter 38 int main(int argc, char* argv[])
4035 23 Jan 21 peter 39 {
4035 23 Jan 21 peter 40   test::Suite suite(argc, argv);
4035 23 Jan 21 peter 41
4035 23 Jan 21 peter 42   std::vector<double> x(10001);
4035 23 Jan 21 peter 43   for (size_t i=0; i<x.size(); ++i)
4035 23 Jan 21 peter 44     x[i] = i;
4036 24 Jan 21 peter 45   double p = 25;
4036 24 Jan 21 peter 46   double q1 = statistics::percentile2(x.begin(), x.end(), p, true);
4036 24 Jan 21 peter 47   suite.out() << "true 1st quartile: " << q1 << "\n";
4035 23 Jan 21 peter 48
4036 24 Jan 21 peter 49   size_t N = 100;
4036 24 Jan 21 peter 50   double alpha = 0.05;
4036 24 Jan 21 peter 51
4036 24 Jan 21 peter 52   double left_tail = 0;
4036 24 Jan 21 peter 53   double right_tail = 0;
4036 24 Jan 21 peter 54   for (size_t k=0; k<=N; ++k) {
4036 24 Jan 21 peter 55     double cdf = gsl_cdf_binomial_P(k, 0.01*p, N);
4036 24 Jan 21 peter 56     suite.out() << k << "\t" << gsl_ran_binomial_pdf(k, 0.01*p, N)
4036 24 Jan 21 peter 57                 << "\t" << cdf << "\n";
4036 24 Jan 21 peter 58     if (cdf <= alpha/2) {
4036 24 Jan 21 peter 59       left_tail = cdf;
4036 24 Jan 21 peter 60     }
4036 24 Jan 21 peter 61     else if (cdf >= 1.0 - alpha/2) {
4036 24 Jan 21 peter 62       right_tail = cdf;
4036 24 Jan 21 peter 63       break;
4036 24 Jan 21 peter 64     }
4036 24 Jan 21 peter 65   }
4036 24 Jan 21 peter 66   suite.out() << "left tail: " << left_tail << "\n";
4036 24 Jan 21 peter 67   suite.out() << "right tail: " << 1.0-right_tail << "\n";
4036 24 Jan 21 peter 68
4035 23 Jan 21 peter 69   statistics::Averager lower;
4035 23 Jan 21 peter 70   statistics::Averager inside;
4035 23 Jan 21 peter 71   statistics::Averager upper;
4036 24 Jan 21 peter 72   statistics::PercentileConfidenceInterval interval(p, 1.0-alpha);
4035 23 Jan 21 peter 73   for (size_t epoch=0; epoch<10000; ++epoch) {
4036 24 Jan 21 peter 74     std::vector<double> y(N);
4036 24 Jan 21 peter 75     random::copy_k_of_n(x.begin(), N, x.size(), y.begin());
4036 24 Jan 21 peter 76     interval(y.begin(), y.end(), true);
4036 24 Jan 21 peter 77     if (q1 < interval.lower()) {
4035 23 Jan 21 peter 78       lower.add(1.0);
4035 23 Jan 21 peter 79       upper.add(0.0);
4035 23 Jan 21 peter 80       inside.add(0.0);
4035 23 Jan 21 peter 81     }
4036 24 Jan 21 peter 82     else if (q1 > interval.upper()) {
4035 23 Jan 21 peter 83       lower.add(0.0);
4035 23 Jan 21 peter 84       upper.add(1.0);
4035 23 Jan 21 peter 85       inside.add(0.0);
4035 23 Jan 21 peter 86     }
4035 23 Jan 21 peter 87     else {
4035 23 Jan 21 peter 88       lower.add(0.0);
4035 23 Jan 21 peter 89       upper.add(0.0);
4035 23 Jan 21 peter 90       inside.add(1.0);
4035 23 Jan 21 peter 91     }
4035 23 Jan 21 peter 92   }
4035 23 Jan 21 peter 93   suite.out() << 100*lower.mean() << "%\t"
4035 23 Jan 21 peter 94               << 100*inside.mean() << "%\t"
4035 23 Jan 21 peter 95               << 100*upper.mean() << "%\n";
4036 24 Jan 21 peter 96   double slack = 0.002;
4036 24 Jan 21 peter 97   if (!suite.add(suite.equal_fix(lower.mean(), left_tail, slack)))
4035 23 Jan 21 peter 98     suite.err() << "incorrect inside mean\n";
4036 24 Jan 21 peter 99   if (!suite.add(suite.equal_fix(inside.mean(), right_tail-left_tail, 2*slack)))
4035 23 Jan 21 peter 100     suite.err() << "incorrect inside mean\n";
4036 24 Jan 21 peter 101   if (!suite.add(suite.equal_fix(upper.mean(), 1.0-right_tail, slack)))
4035 23 Jan 21 peter 102     suite.err() << "incorrect upper mean\n";
4035 23 Jan 21 peter 103   return suite.return_value();
4035 23 Jan 21 peter 104 }

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