Changeset 1418
- Timestamp:
- Aug 20, 2008, 1:15:32 AM (14 years ago)
- Location:
- trunk
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/test/Makefile.am
r1407 r1418 48 48 49 49 TESTS = $(check_PROGRAMS) 50 51 XFAIL_TESTS = statistics_test 50 52 51 53 if HAVE_DOXYGEN -
trunk/test/statistics_test.cc
r1404 r1418 33 33 #include "yat/utility/Vector.h" 34 34 35 #include < vector>35 #include <cmath> 36 36 #include <cstdlib> 37 37 #include <iostream> 38 #include <cmath> 38 #include <map> 39 #include <vector> 40 41 using namespace theplu::yat; 42 void test_percentiler(test::Suite&); 43 44 template<typename RandomAccessIterator> 45 void test_percentiler(test::Suite&, RandomAccessIterator, 46 RandomAccessIterator, 47 double p, double correct); 48 49 template<typename RandomAccessIterator1, typename RandomAccessIterator2> 50 void cmp_percentiler(test::Suite&, 51 RandomAccessIterator1, 52 RandomAccessIterator1, 53 RandomAccessIterator2, 54 RandomAccessIterator2); 39 55 40 56 int main(int argc, char* argv[]) 41 57 { 42 using namespace theplu::yat;43 58 test::Suite suite(argc, argv); 44 59 … … 58 73 statistics::median(data.begin(), data.end()); 59 74 // testing percentile2 60 data.resize(5); 61 std::copy(gsl_vec.begin(), gsl_vec.begin()+5, data.begin()); 62 std::map<double, double> correct; 63 correct[0]=0; 64 correct[10]=0; 65 correct[30]=1; 66 correct[40]=1.5; 67 correct[50]=2; 68 correct[70]=3; 69 correct[90]=4; 70 correct[100]=4; 71 for (std::map<double, double>::const_iterator i=correct.begin(); 72 i!=correct.end(); ++i) { 73 if (!suite.add(suite.equal(statistics::percentile2(data.begin(), data.end(), 74 i->first), 75 i->second))) { 76 suite.err() << "error: unweighted percentile2 " << i->first 77 << "th percentile = " 78 << statistics::percentile2(data.begin(),data.end(),i->first) 79 << " expected to be " << i->second 80 << std::endl; 81 } 82 } 83 84 utility::MatrixWeighted mw(1,data.size(),1); 85 for (size_t i=0; i<mw.columns(); ++i) 86 mw(0,i).data()=data[i]; 87 if (!suite.equal(statistics::median(mw.begin_row(0), mw.end_row(0), true),2)) 88 suite.err() << "error: weighted median sorted failed" << std::endl; 89 if (!suite.equal(statistics::median(mw.begin_row(0), mw.end_row(0), false),2)) 90 suite.err() << "error: weighted median failed" << std::endl; 91 for (std::map<double, double>::const_iterator i=correct.begin(); 92 i!=correct.end(); ++i) { 93 if (!suite.add(suite.equal(statistics::percentile2(mw.begin_row(0), 94 mw.end_row(0), 95 i->first), 96 i->second))) { 97 suite.err() << "error: weighted percentile2 " << i->first 98 << "th percentile = " 99 << statistics::percentile2(mw.begin_row(0), 100 mw.end_row(0), i->first) 101 << " expected to be " << i->second 102 << std::endl; 103 } 104 } 75 test_percentiler(suite); 105 76 106 77 double skewness_gsl=statistics::skewness(gsl_vec); … … 119 90 return suite.return_value(); 120 91 } 92 93 void test_percentiler(test::Suite& suite) 94 { 95 suite.err() << "testing unweighted percentile2" << std::endl; 96 std::vector<double> x; 97 x.reserve(6); 98 for (unsigned int i=0; i<5; i++){ 99 x.push_back(static_cast<double>(i)); 100 } 101 test_percentiler(suite, x.begin(), x.end(), 50, 2); 102 x.push_back(5); 103 test_percentiler(suite, x.begin(), x.end(), 50, 2.5); 104 test_percentiler(suite, x.begin(), x.end(), 25, 1); 105 test_percentiler(suite, x.begin(), x.end(), 0, 0); 106 test_percentiler(suite, x.begin(), x.end(), 10, 0); 107 108 suite.err() << "testing duplication of data\n"; 109 std::vector<double> x2(x); 110 for (size_t i=0; i<x.size(); ++i) 111 x2.push_back(x[i]); 112 cmp_percentiler(suite, x.begin(), x.end(), x2.begin(), x2.end()); 113 114 115 // testing weighted 116 117 suite.err() << "testing weighted percentile2" << std::endl; 118 std::vector<utility::DataWeight> xw(x.size()); 119 for (size_t i=0; i<xw.size(); ++i) { 120 xw[i].data() = x[i]; 121 xw[i].weight() = 1.0; 122 } 123 const std::vector<utility::DataWeight> xw_orig(xw); 124 suite.err() << "testing weighted with unity weights" << std::endl; 125 cmp_percentiler(suite, x.begin(), x.end(), xw.begin(), xw.end()); 126 127 suite.err() << "testing that w=0 equals removed data point\n"; 128 xw=xw_orig; 129 std::vector<utility::DataWeight> xw2(xw_orig); 130 xw[3].weight() = 0.0; 131 xw2.erase(xw2.begin()+3); 132 cmp_percentiler(suite, xw.begin(), xw.end(), xw2.begin(), xw2.end()); 133 134 suite.err() << "testing rescaling of weights\n"; 135 xw2 = xw; 136 for (size_t i=0; i<xw2.size(); ++i) 137 xw2[i].weight()*=2; 138 cmp_percentiler(suite, xw.begin(), xw.end(), xw2.begin(), xw2.end()); 139 140 } 141 142 template<typename RandomAccessIterator> 143 void test_percentiler(test::Suite& suite, 144 RandomAccessIterator first, 145 RandomAccessIterator last, 146 double p, double correct) 147 { 148 using statistics::percentile2; 149 double x = percentile2(first, last, p); 150 if (!suite.add(suite.equal(x, correct, 10))) { 151 suite.err() << "Error in percentile2\n"; 152 suite.err() << " calculated value: " << x << "\n"; 153 suite.err() << " expected value: " << correct << "\n"; 154 } 155 } 156 157 template<typename RandomAccessIterator1, typename RandomAccessIterator2> 158 void cmp_percentiler(test::Suite& suite, 159 RandomAccessIterator1 first1, 160 RandomAccessIterator1 last1, 161 RandomAccessIterator2 first2, 162 RandomAccessIterator2 last2) 163 { 164 for (double p=0; p<100; p+=10) { 165 double correct=statistics::percentile2(first1, last1, p); 166 test_percentiler(suite, first2, last2, p, correct); 167 } 168 169 } -
trunk/yat/statistics/Percentiler.h
r1404 r1418 106 106 // range is one value only is a special case 107 107 if (first+1 == last) 108 return utility::iterator_traits<RandomAccessIterator>().data(first);108 *first; 109 109 if (sorted) { 110 110 size_t n = last - first; 111 111 // have to take care of this special case 112 if (perc_>= 100.0 - 50.0/n)112 if (perc_>= 100.0) 113 113 return *(--last); 114 if (perc_<= 50.0/n)114 if (perc_<= 0.0) 115 115 return *first; 116 double j = n * (perc_-50.0/n) / 100.0; 117 int i = static_cast<int>(j); 118 return (1-j+floor(j))*first[i] + (j-floor(j))*first[i+1]; 116 double j = n * perc_ / 100.0; 117 size_t i = static_cast<size_t>(j); 118 if (i==j) 119 return (first[i]+first[i-1])/2; 120 return first[i]; 119 121 } 120 122 … … 142 144 std::back_inserter(accum_w)); 143 145 144 double sum_w=0;145 146 double w_bound=perc_/100.0*accum_w.back(); 146 RandomAccessIterator upper(first);147 while (sum_w + 0.5*trait.weight(upper) < w_bound) {148 sum_w += trait.weight(upper);147 std::vector<double>::const_iterator upper(accum_w.begin()); 148 double allowed_error=1e-10; 149 while (upper != accum_w.end() && *upper < w_bound + allowed_error) 149 150 ++upper; 150 } 151 if (upper==last) 152 return trait.data(--last); 153 if (sum_w==0.0) { 154 while (first!=last) 155 if (trait.weight(first)) 156 return trait.data(first); 157 // all weights zero, hmm 158 return 0; 159 } 160 // lower should have non-zero weight 161 // and (lower, upper) should have zero weights 162 RandomAccessIterator lower=upper-1; 163 while (trait.weight(lower)==0) 151 std::vector<double>::const_iterator lower(upper); 152 while (lower!=accum_w.begin() && *lower > w_bound - allowed_error) 164 153 --lower; 165 166 double p1=100.0/accum_w.back()*(sum_w + 0.5*trait.weight(upper)); 167 double p0=100.0/accum_w.back()*(sum_w - 0.5*trait.weight(lower)); 168 double alpha = (perc_-p0)/(p1-p0); 169 utility::yat_assert<std::runtime_error>(alpha<=1.0, 170 "Percentile: alpha out of bounds"); 171 return trait.data(lower)+alpha*(trait.data(upper)-trait.data(lower)); 154 return (trait.data(first+(upper-accum_w.begin()))+ 155 trait.data(first+(lower-accum_w.begin()+1)))/2; 172 156 } 173 157
Note: See TracChangeset
for help on using the changeset viewer.