yat/statistics/PercentileConfidenceInterval.h

Code
Comments
Other
Rev Date Author Line
4035 23 Jan 21 peter 1 #ifndef _theplu_yat_statistics_percentile_confidence_interval_
4035 23 Jan 21 peter 2 #define _theplu_yat_statistics_percentile_confidence_interval_
4035 23 Jan 21 peter 3
4035 23 Jan 21 peter 4 // $Id$
4035 23 Jan 21 peter 5
4035 23 Jan 21 peter 6 /*
4035 23 Jan 21 peter 7   Copyright (C) 2021 Peter Johansson
4035 23 Jan 21 peter 8
4035 23 Jan 21 peter 9   This file is part of the yat library, http://dev.thep.lu.se/yat
4035 23 Jan 21 peter 10
4035 23 Jan 21 peter 11   The yat library is free software; you can redistribute it and/or
4035 23 Jan 21 peter 12   modify it under the terms of the GNU General Public License as
4035 23 Jan 21 peter 13   published by the Free Software Foundation; either version 3 of the
4035 23 Jan 21 peter 14   License, or (at your option) any later version.
4035 23 Jan 21 peter 15
4035 23 Jan 21 peter 16   The yat library is distributed in the hope that it will be useful,
4035 23 Jan 21 peter 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
4035 23 Jan 21 peter 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
4035 23 Jan 21 peter 19   General Public License for more details.
4035 23 Jan 21 peter 20
4035 23 Jan 21 peter 21   You should have received a copy of the GNU General Public License
4035 23 Jan 21 peter 22   along with yat. If not, see <http://www.gnu.org/licenses/>.
4035 23 Jan 21 peter 23 */
4035 23 Jan 21 peter 24
4035 23 Jan 21 peter 25 #include "yat/utility/yat_assert.h"
4035 23 Jan 21 peter 26
4035 23 Jan 21 peter 27 #include <boost/concept_check.hpp>
4035 23 Jan 21 peter 28 #include <boost/iterator/iterator_concepts.hpp>
4035 23 Jan 21 peter 29
4036 24 Jan 21 peter 30 #include <boost/math/distributions/binomial.hpp>
4036 24 Jan 21 peter 31
4035 23 Jan 21 peter 32 #include <algorithm>
4035 23 Jan 21 peter 33 #include <vector>
4035 23 Jan 21 peter 34
4035 23 Jan 21 peter 35 namespace theplu {
4035 23 Jan 21 peter 36 namespace yat {
4035 23 Jan 21 peter 37 namespace statistics {
4035 23 Jan 21 peter 38
4035 23 Jan 21 peter 39   /**
4036 24 Jan 21 peter 40      Class calculates the confidence interval. It follows the method
4036 24 Jan 21 peter 41      suggested by Martin Bland in An Introduction to Medical
4036 24 Jan 21 peter 42      Statistics and <a
4036 24 Jan 21 peter 43      href=\"https://www-users.york.ac.uk/~mb55/intro/cicent.htm\">here</a>,
4036 24 Jan 21 peter 44      but uses the exact binomial distribution instead of a
4036 24 Jan 21 peter 45      large-sample normal approximation.
4036 24 Jan 21 peter 46
4036 24 Jan 21 peter 47      You need to load the class with data via the operator() before
4036 24 Jan 21 peter 48      the interval is available via the lower(void) and upper(void
4036 24 Jan 21 peter 49      functions.
4036 24 Jan 21 peter 50
4036 24 Jan 21 peter 51      A weighted input is currently not supported.
4036 24 Jan 21 peter 52
4035 23 Jan 21 peter 53      \since New in yat 0.19
4035 23 Jan 21 peter 54    */
4035 23 Jan 21 peter 55   class PercentileConfidenceInterval
4035 23 Jan 21 peter 56   {
4035 23 Jan 21 peter 57   public:
4035 23 Jan 21 peter 58     /**
4036 24 Jan 21 peter 59        \param k The kth percentile, for example, 50 for median
4036 24 Jan 21 peter 60        \param level Confidence level, for example, 0.95 for 95%
4036 24 Jan 21 peter 61        confidence interval
4036 24 Jan 21 peter 62      */
4036 24 Jan 21 peter 63     PercentileConfidenceInterval(double k, double level);
4036 24 Jan 21 peter 64
4036 24 Jan 21 peter 65
4036 24 Jan 21 peter 66     /**
4035 23 Jan 21 peter 67        \param first First element in range
4035 23 Jan 21 peter 68        \param last One past element in range
4035 23 Jan 21 peter 69        \param sorted if \c true the range is assumed to be sorted;
4035 23 Jan 21 peter 70        otherwise the range is copied, the copy sorted, before the
4035 23 Jan 21 peter 71        confidence interval calculated.
4035 23 Jan 21 peter 72
4035 23 Jan 21 peter 73        Type Requirements:
4035 23 Jan 21 peter 74        - \c RandomAccessIterator must be a \random_access_traversal_iterator
4035 23 Jan 21 peter 75        - \c RandomAccessIterator must be a \readable_iterator
4035 23 Jan 21 peter 76        - \c value_type of RandomAccessIterator is convertible to \c double
4035 23 Jan 21 peter 77      */
4035 23 Jan 21 peter 78     template<typename RandomAccessIterator>
4036 24 Jan 21 peter 79     void operator()(RandomAccessIterator first, RandomAccessIterator last,
4036 24 Jan 21 peter 80                     bool sorted=false)
4035 23 Jan 21 peter 81     {
4035 23 Jan 21 peter 82       typedef RandomAccessIterator T;
4035 23 Jan 21 peter 83       typedef typename std::iterator_traits<T>::value_type value_type;
4035 23 Jan 21 peter 84       BOOST_CONCEPT_ASSERT((boost_concepts::RandomAccessTraversal<T>));
4035 23 Jan 21 peter 85       BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<T>));
4035 23 Jan 21 peter 86       BOOST_CONCEPT_ASSERT((boost::Convertible<value_type, double>));
4036 24 Jan 21 peter 87       prepare(last - first);
4035 23 Jan 21 peter 88       process(first, last, sorted);
4035 23 Jan 21 peter 89     }
4035 23 Jan 21 peter 90
4035 23 Jan 21 peter 91     /**
4035 23 Jan 21 peter 92        \return alpha, for example 95 for 95% confidence interval
4035 23 Jan 21 peter 93      */
4036 24 Jan 21 peter 94     double alpha(void) const;
4035 23 Jan 21 peter 95
4035 23 Jan 21 peter 96     /**
4035 23 Jan 21 peter 97        \return lower confidence interval
4035 23 Jan 21 peter 98      */
4036 24 Jan 21 peter 99     double lower(void) const;
4035 23 Jan 21 peter 100
4035 23 Jan 21 peter 101     /**
4035 23 Jan 21 peter 102        For example, 50, for median.
4035 23 Jan 21 peter 103      */
4036 24 Jan 21 peter 104     double k(void) const;
4035 23 Jan 21 peter 105
4035 23 Jan 21 peter 106     /**
4035 23 Jan 21 peter 107        \return upper confidence interval
4035 23 Jan 21 peter 108      */
4036 24 Jan 21 peter 109     double upper(void) const;
4035 23 Jan 21 peter 110
4035 23 Jan 21 peter 111   private:
4035 23 Jan 21 peter 112     double alpha_;
4035 23 Jan 21 peter 113     double k_;
4036 24 Jan 21 peter 114     size_t n_;
4036 24 Jan 21 peter 115     size_t lower_index_;
4035 23 Jan 21 peter 116     double lower_;
4036 24 Jan 21 peter 117     size_t upper_index_;
4035 23 Jan 21 peter 118     double upper_;
4035 23 Jan 21 peter 119
4036 24 Jan 21 peter 120     void prepare(size_t n);
4036 24 Jan 21 peter 121
4035 23 Jan 21 peter 122     template<typename RandomAccessIterator>
4035 23 Jan 21 peter 123     void process(RandomAccessIterator first, RandomAccessIterator last,
4035 23 Jan 21 peter 124                  bool sorted)
4035 23 Jan 21 peter 125     {
4035 23 Jan 21 peter 126       if (!sorted) {
4035 23 Jan 21 peter 127         std::vector<double> v(first, last);
4035 23 Jan 21 peter 128         std::sort(v.begin(), v.end());
4035 23 Jan 21 peter 129         process(v.begin(), v.end(), true);
4036 24 Jan 21 peter 130         return;
4035 23 Jan 21 peter 131       }
4036 24 Jan 21 peter 132       lower_ = first[lower_index_];
4036 24 Jan 21 peter 133       upper_ = first[upper_index_];
4035 23 Jan 21 peter 134     }
4035 23 Jan 21 peter 135
4035 23 Jan 21 peter 136   };
4035 23 Jan 21 peter 137
4035 23 Jan 21 peter 138 }}} // of namespace statistics, yat, and theplu
4035 23 Jan 21 peter 139
4035 23 Jan 21 peter 140 #endif