Changeset 4036

Timestamp:

Jan 24, 2021, 2:02:57 AM (2 years ago)

Author:

Peter

Message:

closes #963. New class PercentileConfidenceInterval?

Location:

trunk

Files:

• test/Makefile.am (modified) (1 diff)
• test/percentile_confidence_interval.cc (modified) (3 diffs)
• yat/statistics/Makefile.am (modified) (1 diff)
• yat/statistics/PercentileConfidenceInterval.cc (added)
• yat/statistics/PercentileConfidenceInterval.h (modified) (7 diffs)

trunk/test/Makefile.am

@@ -135 +135 @@
 
 # tests not passing through yet
-XFAIL_TESTS = test/percentile_confidence_interval.test
+XFAIL_TESTS =
 
 DISTRIBUTED_TESTS = \

trunk/test/percentile_confidence_interval.cc

@@ -29 +29 @@
 #include <yat/random/copy_k_of_n.h>
 
+#include <gsl/gsl_cdf.h>
+#include <gsl/gsl_randist.h>
+
 #include <set>
 
@@ -40 +43 @@
   for (size_t i=0; i<x.size(); ++i)
     x[i] = i;
-  double m = statistics::median(x.begin(), x.end());
-  suite.out() << "true median: " << m << "\n";
+  double p = 25;
+  double q1 = statistics::percentile2(x.begin(), x.end(), p, true);
+  suite.out() << "true 1st quartile: " << q1 << "\n";
+
+  size_t N = 100;
+  double alpha = 0.05;
+
+  double left_tail = 0;
+  double right_tail = 0;
+  for (size_t k=0; k<=N; ++k) {
+    double cdf = gsl_cdf_binomial_P(k, 0.01*p, N);
+    suite.out() << k << "\t" << gsl_ran_binomial_pdf(k, 0.01*p, N)
+                << "\t" << cdf << "\n";
+    if (cdf <= alpha/2) {
+      left_tail = cdf;
+    }
+    else if (cdf >= 1.0 - alpha/2) {
+      right_tail = cdf;
+      break;
+    }
+  }
+  suite.out() << "left tail: " << left_tail << "\n";
+  suite.out() << "right tail: " << 1.0-right_tail << "\n";
 
   statistics::Averager lower;
   statistics::Averager inside;
   statistics::Averager upper;
-  double alpha = 0.95;
+  statistics::PercentileConfidenceInterval interval(p, 1.0-alpha);
   for (size_t epoch=0; epoch<10000; ++epoch) {
-    std::vector<double> y(100);
-    random::copy_k_of_n(x.begin(), y.size(), x.size(), y.begin());
-    statistics::PercentileConfidenceInterval
-      ci(y.begin(), y.end(), 50.0, 100*alpha);
-    if (m < ci.lower()) {
+    std::vector<double> y(N);
+    random::copy_k_of_n(x.begin(), N, x.size(), y.begin());
+    interval(y.begin(), y.end(), true);
+    if (q1 < interval.lower()) {
       lower.add(1.0);
       upper.add(0.0);
       inside.add(0.0);
     }
-    else if (m > ci.upper()) {
+    else if (q1 > interval.upper()) {
       lower.add(0.0);
       upper.add(1.0);
@@ -71 +94 @@
               << 100*inside.mean() << "%\t"
               << 100*upper.mean() << "%\n";
-  double slack = 0.005;
-  if (!suite.add(suite.equal_fix(lower.mean(), (1.0-alpha)/2, slack)))
+  double slack = 0.002;
+  if (!suite.add(suite.equal_fix(lower.mean(), left_tail, slack)))
     suite.err() << "incorrect inside mean\n";
-  if (!suite.add(suite.equal_fix(inside.mean(), alpha, 2*slack)))
+  if (!suite.add(suite.equal_fix(inside.mean(), right_tail-left_tail, 2*slack)))
     suite.err() << "incorrect inside mean\n";
-  if (!suite.add(suite.equal_fix(upper.mean(), (1.0-alpha)/2, slack)))
+  if (!suite.add(suite.equal_fix(upper.mean(), 1.0-right_tail, slack)))
     suite.err() << "incorrect upper mean\n";
   return suite.return_value();

trunk/yat/statistics/Makefile.am

@@ -36 +36 @@
   yat/statistics/LikelihoodRatioTestBinomial.cc \
   yat/statistics/Pearson.cc \
-  yat/statistics/PearsonCorrelation.cc yat/statistics/Percentiler.cc \
+  yat/statistics/PearsonCorrelation.cc \
+  yat/statistics/PercentileConfidenceInterval.cc \
+  yat/statistics/Percentiler.cc \
   yat/statistics/ROC.cc \
   yat/statistics/SAMScore.cc yat/statistics/Score.cc \

trunk/yat/statistics/PercentileConfidenceInterval.h

@@ -28 +28 @@
 #include <boost/iterator/iterator_concepts.hpp>
 
+#include <boost/math/distributions/binomial.hpp>
+
 #include <algorithm>
 #include <vector>
@@ -36 +38 @@
 
   /**
+     Class calculates the confidence interval. It follows the method
+     suggested by Martin Bland in An Introduction to Medical
+     Statistics and <a
+     href=\"https://www-users.york.ac.uk/~mb55/intro/cicent.htm\">here</a>,
+     but uses the exact binomial distribution instead of a
+     large-sample normal approximation.
+
+     You need to load the class with data via the operator() before
+     the interval is available via the lower(void) and upper(void
+     functions.
+
+     A weighted input is currently not supported.
+
      \since New in yat 0.19
    */
@@ -42 +57 @@
   public:
     /**
+       \param k The kth percentile, for example, 50 for median
+       \param level Confidence level, for example, 0.95 for 95%
+       confidence interval
+     */
+    PercentileConfidenceInterval(double k, double level);
+
+
+    /**
        \param first First element in range
        \param last One past element in range
-       \param k The kth percentile, for example, 50 for median
-       \param alpha Level of confidence, for example, 0.95 for 95%
-       confidence interval
        \param sorted if \c true the range is assumed to be sorted;
        otherwise the range is copied, the copy sorted, before the
@@ -57 +77 @@
      */
     template<typename RandomAccessIterator>
-    PercentileConfidenceInterval(RandomAccessIterator first,
-                                 RandomAccessIterator last,
-                                 double k, double alpha,
-                                 bool sorted=false)
-      : alpha_(alpha), k_(k)
+    void operator()(RandomAccessIterator first, RandomAccessIterator last,
+                    bool sorted=false)
     {
       typedef RandomAccessIterator T;
@@ -68 +85 @@
       BOOST_CONCEPT_ASSERT((boost_concepts::ReadableIterator<T>));
       BOOST_CONCEPT_ASSERT((boost::Convertible<value_type, double>));
-
+      prepare(last - first);
       process(first, last, sorted);
     }
@@ -75 +92 @@
        \return alpha, for example 95 for 95% confidence interval
      */
-    double alpha(void) const { return alpha_; }
+    double alpha(void) const;
 
     /**
        \return lower confidence interval
      */
-    double lower(void) const { return lower_; }
+    double lower(void) const;
 
     /**
        For example, 50, for median.
      */
-    double k(void) const { return k_; }
+    double k(void) const;
 
     /**
        \return upper confidence interval
      */
-    double upper(void) const { return upper_; }
+    double upper(void) const;
 
   private:
     double alpha_;
     double k_;
+    size_t n_;
+    size_t lower_index_;
     double lower_;
+    size_t upper_index_;
     double upper_;
+
+    void prepare(size_t n);
 
     template<typename RandomAccessIterator>
@@ -106 +128 @@
         std::sort(v.begin(), v.end());
         process(v.begin(), v.end(), true);
+        return;
       }
-
-      // number of observations smaller than the percentile is an
-      // obervatioon from a Binomal distribution, Bin(n, p), where p =
-      // 0.01*perc_.
-
-      size_t n = last - first;
-      size_t lower_idx = 0;
-      size_t upper_idx = n - 1;
-      lower_ = first[lower_idx];
-      upper_ = first[upper_idx];
+      lower_ = first[lower_index_];
+      upper_ = first[upper_index_];
     }