|
23 | 23 |
|
24 | 24 | #pragma once |
25 | 25 |
|
| 26 | +#include <boost/accumulators/accumulators.hpp> |
| 27 | +#include <boost/accumulators/statistics/mean.hpp> |
| 28 | +#include <boost/accumulators/statistics/stats.hpp> |
| 29 | +#include <boost/accumulators/statistics/variance.hpp> |
| 30 | + |
26 | 31 | namespace QuantExt { |
27 | 32 |
|
28 | 33 | template <class I> |
29 | | -void distributionCount(I begin, I end, const Size steps, std::vector<Real>& bounds, std::vector<Size>& counts) { |
30 | | - Real xmin = *std::min_element(begin, end); |
31 | | - Real xmax = *std::max_element(begin, end); |
| 34 | +void distributionCount(I begin, I end, const Size steps, std::vector<Real>& bounds, std::vector<Size>& counts, |
| 35 | + const Real coveredStdDevs = Null<Real>()) { |
| 36 | + |
32 | 37 | std::vector<Real> v(begin, end); |
33 | 38 | std::sort(v.begin(), v.end()); |
| 39 | + |
| 40 | + Real xmin, xmax; |
| 41 | + if (coveredStdDevs == Null<Real>()) { |
| 42 | + |
| 43 | + // cover [xmin, xmax] in output |
| 44 | + |
| 45 | + xmin = *std::min_element(begin, end); |
| 46 | + xmax = *std::max_element(begin, end); |
| 47 | + |
| 48 | + } else { |
| 49 | + |
| 50 | + // cover given number of std devs around mean in output |
| 51 | + |
| 52 | + boost::accumulators::accumulator_set< |
| 53 | + double, boost::accumulators::stats<boost::accumulators::tag::mean, boost::accumulators::tag::variance>> |
| 54 | + acc; |
| 55 | + std::for_each(v.begin(), v.end(), [&acc](double x) { acc(x); }); |
| 56 | + double mean = boost::accumulators::mean(acc); |
| 57 | + double stdDev = std::sqrt(boost::accumulators::variance(acc)); |
| 58 | + xmin = mean - coveredStdDevs * stdDev; |
| 59 | + xmax = mean + coveredStdDevs * stdDev; |
| 60 | + } |
| 61 | + |
34 | 62 | Real h = (xmax - xmin) / static_cast<Real>(steps); |
35 | | - Size idx0 = 0; |
| 63 | + Size idx0 = coveredStdDevs == Null<Real>() ? 0 : std::upper_bound(v.begin(), v.end(), xmin) - v.begin(); |
36 | 64 | counts.resize(steps); |
37 | 65 | bounds.resize(steps); |
38 | 66 | for (Size i = 0; i < steps; ++i) { |
39 | 67 | Real v1 = xmin + static_cast<Real>(i + 1) * h; |
40 | | - // the code for "i == steps - 1" ensures all observations are accounted for |
41 | | - Size idx1 = i == steps - 1 ? v.size() : std::upper_bound(v.begin(), v.end(), v1) - v.begin(); |
| 68 | + // if all data should be covered, ensure the counts sum up to v.size() |
| 69 | + Size idx1 = coveredStdDevs == Null<Real>() && i == steps - 1 |
| 70 | + ? v.size() |
| 71 | + : std::upper_bound(v.begin(), v.end(), v1) - v.begin(); |
42 | 72 | counts[i] = idx1 - idx0; |
43 | 73 | bounds[i] = v1; |
44 | 74 | idx0 = idx1; |
|
0 commit comments