这个 Quantile.h 文件是一个分位数计算库，主要用于实时数据流的分位数估计。 基于 Andrey Akinshin 滑动窗口扩展P²分位数估计算法移植到 C++ 并进行了增强。参考原文：aakinshin.net/posts/movin…
详细解释其作用：
核心功能
- 分位数计算：实现高效的流式数据分位数估计算法，能够在不存储全部数据的情况下，实时计算各种分位数值（如中位数、95%分位数等）。
主要组件
- 抽象基类 QuantileBase 定义了分位数计算的标准接口 支持添加数据点、计算CDF、获取分位数等操作
- 具体实现类 MovingExtendedP2AndreyAkinshin 基于滑动窗口扩展P²算法，具有以下特点：
算法优势
- 内存高效：仅维护固定数量的标记点，不存储全部数据
- 实时计算：支持流式数据，逐个数据点处理
- 滑动窗口：支持窗口化处理，适应数据分布变化
- 精度可控：通过标记点数量平衡精度和内存使用 核心特性：
- 多分位数支持：同时计算多个分位数（0%, 1%, 5%, 10%, ..., 99%, 100%）
- 权重支持：支持带权重的数据点
- 边界处理：完善的边界值处理机制
- 内存监控：提供内存使用量统计
应用场景 这个库特别适合以下场景：
- 时间序列分析：实时监控系统指标的分位数
- 性能监控：计算延迟、吞吐量等指标的分位数
- 数据流处理：大数据流中的统计计算
- 资源监控：CPU使用率、内存使用等的分位数监控
技术特点
- P²算法核心：使用抛物线插值和线性插值进行平滑调整
- 滑动窗口机制：结合当前窗口和前一窗口的估计结果
- 自适应调整：根据数据分布动态调整标记点位置
- 边界插值：对边界概率进行特殊处理，保证计算稳定性
#pragma once
#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <map>
#include <random>
#include <stdexcept>
#include <string>
#include <vector>
// 分位数计算抽象基类
class QuantileBase {
public:
virtual void add(double x) = 0; // 添加单个数据点
virtual void add(double x, double w) = 0; // 添加带权重的数据点
virtual double cdf(double x) = 0; // 计算累积分布函数值
virtual double quantile(double q) = 0; // 计算分位数值
[[nodiscard]] virtual long totalWeight() const = 0; // 获取总权重（数据点数）
[[nodiscard]] virtual long memoryUsage() const = 0; // 获取内存使用量
virtual void reset() = 0; // 重置估计器状态
virtual ~QuantileBase() = default; // 虚析构函数
};
/**
* 基于Andrey Akinshin提出的滑动窗口扩展P²分位数估计算法
* 论文参考：https://aakinshin.net/posts/moving-ex-p2-quantile-estimator/
*
* 算法原理：
* 1. 使用P²算法维护一组标记点(markers)来估计分位数
* 2. 支持滑动窗口，结合当前窗口和前一窗口的估计结果
* 3. 通过线性插值和平滑过渡处理窗口边界
*/
class MovingExtendedP2AndreyAkinshin : public QuantileBase {
public:
/**
* 构造函数
* @param window 滑动窗口大小，0表示不使用滑动窗口
* @param probs 目标概率值列表，默认为常用分位数概率
*/
explicit MovingExtendedP2AndreyAkinshin(int window = 0,
const std::vector<double>& probs = {0, 0.01, 0.05, 0.1, 0.15, 0.2, 0.25,
0.3, 0.4, 0.5, 0.75, 0.9, 0.95, 0.99,
1})
: windowSize(window), // 滑动窗口大小
probabilities(probs), // 目标概率列表
m(static_cast<int>(probs.size())), // 目标概率个数
markerCount(2 * m + 3), // 标记点总数：每个概率对应2个标记点 + 3个边界点
n(0), // 当前数据点计数
currentMin(std::numeric_limits<double>::max()), // 全局最小值
currentMax(std::numeric_limits<double>::lowest()), // 全局最大值
windowMin(std::numeric_limits<double>::max()), // 当前窗口最小值
windowMax(std::numeric_limits<double>::lowest()), // 当前窗口最大值
previousWindowEstimations(m, 0.0), // 前一窗口的分位数估计值
n_pos(markerCount, 0), // 标记点的实际位置索引
ns(markerCount, 0.0), // 标记点的期望位置
q(markerCount, 0.0), // 标记点对应的分位数值
totalCount(0), // 内部估计器数据点数
initialized(false) { // 初始化标志
// 参数验证
if (m <= 0) {
throw std::invalid_argument("Probability list cannot be empty");
}
if (markerCount < 5) {
throw std::invalid_argument("Too few probabilities");
}
}
// 添加单个数据点（权重为1）
void add(double x) override { add(x, 1.0); }
/**
* 添加带权重的数据点
* 实现说明：将权重转换为重复添加次数（简化实现）
* 注意：生产环境应考虑更精确的权重处理方式
*/
void add(double x, double w) override {
if (w <= 0) return;
// 简化实现：将权重转换为重复添加
int repeatCount = static_cast<int>(std::round(w));
for (int i = 0; i < repeatCount; i++) {
addSingleDataPoint(x);
}
}
/**
* 计算指定概率的分位数值
* 算法流程：
* 1. 参数验证和边界处理
* 2. 非滑动窗口模式直接返回内部估计值
* 3. 滑动窗口模式结合当前窗口和前一窗口的估计值
* 4. 使用线性插值处理不在预定义概率列表中的分位数
*/
double quantile(double p) override {
if (p < 0 || p > 1) {
throw std::invalid_argument("Quantile must be between 0 and 1");
}
if (n == 0) {
throw std::runtime_error("No data points added");
}
// 边界分位数特殊处理：直接返回最小/最大值
if (p == 0.0) return getMinValue();
if (p == 1.0) return getMaxValue();
// 非滑动窗口模式回退
if (windowSize == 0) return getInternalQuantile(p);
// 数据量不足窗口大小时直接使用内部估计
if (n < windowSize) return getInternalQuantile(p);
// 查找预定义概率列表中的匹配项
for (int i = 0; i < m; i++) {
if (std::abs(probabilities[i] - p) < 1e-10) {
// 滑动窗口融合：w1*前一窗口 + w2*当前窗口
double estimation1 = previousWindowEstimations[i]; // 前一窗口估计
double estimation2 = getInternalQuantile(p); // 当前窗口估计
double w2 = (n % windowSize + 1) * 1.0 / windowSize; // 当前窗口权重
double w1 = 1.0 - w2; // 前一窗口权重
return w1 * estimation1 + w2 * estimation2; // 加权融合
}
}
// 概率小于最小预定义概率：边界插值处理
if (p <= probabilities[0]) {
double estimation1 = previousWindowEstimations[0];
double estimation2 = getInternalQuantile(p);
double w2 = (n % windowSize + 1) * 1.0 / windowSize;
double w1 = 1.0 - w2;
double blended = w1 * estimation1 + w2 * estimation2;
// 额外边界插值：在0分位数和混合估计值之间线性插值
double fraction = p / probabilities[0];
return blended * fraction + getInternalQuantile(0.0) * (1 - fraction);
}
// 概率大于最大预定义概率：边界插值处理
if (p >= probabilities[m - 1]) {
double estimation1 = previousWindowEstimations[m - 1];
double estimation2 = getInternalQuantile(p);
double w2 = (n % windowSize + 1) * 1.0 / windowSize;
double w1 = 1.0 - w2;
double blended = w1 * estimation1 + w2 * estimation2;
// 额外边界插值：在混合估计值和1分位数之间线性插值
double fraction = (p - probabilities[m - 1]) / (1 - probabilities[m - 1]);
return blended + (getInternalQuantile(1.0) - blended) * fraction;
}
// 中间概率插值：在两个预定义概率之间进行线性插值
for (int i = 0; i < m - 1; i++) {
if (p >= probabilities[i] && p <= probabilities[i + 1]) {
double estimation1_low = previousWindowEstimations[i];
double estimation1_high = previousWindowEstimations[i + 1];
double estimation2 = getInternalQuantile(p);
double w2 = (n % windowSize + 1) * 1.0 / windowSize;
double w1 = 1.0 - w2;
// 对前一窗口的估计值进行线性插值
double fraction1 = (p - probabilities[i]) / (probabilities[i + 1] - probabilities[i]);
double blended_estimation1 = estimation1_low + fraction1 * (estimation1_high - estimation1_low);
// 融合两个窗口的估计值
return w1 * blended_estimation1 + w2 * estimation2;
}
}
throw std::runtime_error("Failed to compute quantile");
}
/**
* 计算累积分布函数值（经验CDF）
* 实现说明：基于当前标记点进行线性插值估计
*/
double cdf(double x) override {
if (totalCount == 0) return 0.0;
// 数据量较少时使用精确排序计算
if (totalCount <= markerCount) {
std::vector<double> sortedData(q.begin(), q.begin() + totalCount);
std::sort(sortedData.begin(), sortedData.end());
auto it = std::lower_bound(sortedData.begin(), sortedData.end(), x);
size_t count = std::distance(sortedData.begin(), it);
return static_cast<double>(count) / static_cast<double>(totalCount);
}
// 边界值处理
if (x < q[0]) return 0.0;
if (x >= q[markerCount - 1]) return 1.0;
// 在标记点之间线性插值计算CDF
for (int i = 1; i < markerCount; i++) {
if (x < q[i]) {
double fraction = (x - q[i - 1]) / (q[i] - q[i - 1]);
double cdf_lower = static_cast<double>(n_pos[i - 1]) / static_cast<double>(totalCount);
double cdf_upper = static_cast<double>(n_pos[i]) / static_cast<double>(totalCount);
return cdf_lower + fraction * (cdf_upper - cdf_lower);
}
}
return 1.0;
}
long totalWeight() const override { return n; }
/**
* 计算内存使用量（字节）
* 注意：此为近似值，实际内存使用可能因内存对齐等有所不同
*/
long memoryUsage() const override {
size_t memory = 0;
memory += probabilities.capacity() * sizeof(double);
memory += previousWindowEstimations.capacity() * sizeof(double);
memory += n_pos.capacity() * sizeof(int);
memory += ns.capacity() * sizeof(double);
memory += q.capacity() * sizeof(double);
return static_cast<long>(memory);
}
// 重置估计器状态
void reset() override {
n = 0;
currentMin = std::numeric_limits<double>::max();
currentMax = std::numeric_limits<double>::lowest();
windowMin = std::numeric_limits<double>::max();
windowMax = std::numeric_limits<double>::lowest();
std::fill(previousWindowEstimations.begin(), previousWindowEstimations.end(), 0.0);
resetInternalEstimator();
}
// 获取配置的概率列表
const std::vector<double>& getProbabilities() const { return probabilities; }
private:
// 配置参数
int windowSize; // 滑动窗口大小
std::vector<double> probabilities; // 目标概率值列表
int m; // 目标概率个数
int markerCount; // 标记点总数 (2m + 3)
// 数据统计
int n; // 总数据点数
double currentMin, currentMax; // 全局最小最大值
double windowMin, windowMax; // 当前窗口最小最大值
// 窗口状态管理
std::vector<double> previousWindowEstimations; // 前一窗口的分位数估计值
// P²算法内部状态（核心组件）
std::vector<int> n_pos; // 标记点的实际位置索引
std::vector<double> ns; // 标记点的期望位置（理论位置）
std::vector<double> q; // 标记点对应的分位数值
int totalCount; // 内部估计器数据点数
bool initialized; // 初始化完成标志
/**
* P²算法核心原理说明：
* 1. 维护一组固定数量的标记点(markers)，每个标记点包含：
* - 位置索引(n_pos)：标记点在数据流中的排序位置
* - 分位数值(q)：该位置对应的数据值
* - 期望位置(ns)：基于目标概率的理论位置
*
* 2. 算法通过调整标记点的位置和值来逼近真实分位数
* 3. 使用抛物线插值(Parabolic)和线性插值(Linear)进行平滑调整
*/
private:
// 添加单个数据点（核心处理逻辑）
void addSingleDataPoint(double value) {
n++;
updateMinMax(value);
// 非滑动窗口模式处理
if (windowSize == 0) {
addSingleValue(value);
return;
}
// 窗口切换检查：达到窗口大小时保存统计信息并重置内部估计器
if (n % windowSize == 0) {
saveWindowStatistics();
resetInternalEstimator();
}
addSingleValue(value);
}
// 更新最小最大值统计
void updateMinMax(double value) {
if (value < currentMin) currentMin = value;
if (value > currentMax) currentMax = value;
if (value < windowMin) windowMin = value;
if (value > windowMax) windowMax = value;
}
// 保存当前窗口的统计信息
void saveWindowStatistics() {
for (int i = 0; i < m; i++) {
previousWindowEstimations[i] = getInternalQuantile(probabilities[i]);
}
// 更新窗口边界信息
windowMin = currentMin;
windowMax = currentMax;
}
// 获取最小值（考虑滑动窗口）
double getMinValue() const {
if (windowSize == 0) return currentMin;
if (n < windowSize) return currentMin;
return std::min(windowMin, currentMin);
}
// 获取最大值（考虑滑动窗口）
double getMaxValue() const {
if (windowSize == 0) return currentMax;
if (n < windowSize) return currentMax;
return std::max(windowMax, currentMax);
}
/**
* P²算法核心：添加单个值到内部估计器
* 处理流程：
* 1. 初始化阶段：直接存储数据点
* 2. 确定新数据点的插入位置
* 3. 更新标记点位置信息
* 4. 调整标记点位置和分位数值
*/
void addSingleValue(double value) {
// 初始化阶段：直接存储前markerCount个数据点
if (totalCount < markerCount) {
q[totalCount] = value;
totalCount++;
if (totalCount == markerCount) {
initializeIfNeeded(); // 初始化标记点
}
return;
}
// 确定新数据点的插入区间
int k = -1;
if (value < q[0]) {
q[0] = value; // 更新最小值
k = 0;
} else {
for (int i = 1; i < markerCount; i++) {
if (value < q[i]) {
k = i - 1;
break;
}
}
if (k == -1) {
q[markerCount - 1] = value; // 更新最大值
k = markerCount - 2;
}
}
// 更新标记点位置索引（k之后的标记点位置+1）
for (int i = k + 1; i < markerCount; i++) n_pos[i]++;
UpdateNs(totalCount);
// 调整标记点：从中间向两边处理
int leftI = 1, rightI = markerCount - 2;
while (leftI <= rightI) {
int i;
// 选择距离中位数最近的标记点优先调整
if (std::abs(ns[leftI] / static_cast<double>(totalCount) - 0.5) <=
std::abs(ns[rightI] / static_cast<double>(totalCount) - 0.5))
i = leftI++;
else
i = rightI--;
if (i >= 1 && i < markerCount - 1) {
Adjust(i); // 调整第i个标记点
}
}
totalCount++;
}
/**
* 更新标记点的期望位置
* 标记点布局：
* 索引：0, 1, 2, 3, ..., 2m, 2m+1, 2m+2
* 类型：min, mid1, p0, mid2, ..., pm-1, midm, max
*/
void UpdateNs(int maxIndex) {
if (maxIndex <= 0) return;
// 主标记点位置（对应目标概率）
ns[0] = 0;
for (int i = 0; i < m; i++) {
if (i * 2 + 2 < ns.size()) {
ns[i * 2 + 2] = maxIndex * probabilities[i];
}
}
ns[markerCount - 1] = maxIndex;
// 中间标记点位置（概率之间的中点）
ns[1] = maxIndex * probabilities[0] / 2;
for (int i = 1; i < m; i++) {
if (2 * i + 1 < ns.size()) {
ns[2 * i + 1] = maxIndex * (probabilities[i - 1] + probabilities[i]) / 2;
}
}
ns[markerCount - 2] = maxIndex * (1 + probabilities[m - 1]) / 2;
}
/**
* 调整标记点位置和分位数值（P²算法核心调整逻辑）
* 调整条件：当标记点的期望位置与实际位置偏差超过1时进行调整
*/
void Adjust(int i) {
if (i <= 0 || i >= markerCount - 1) return;
if (i >= n_pos.size() || i >= ns.size()) return;
double d = ns[i] - n_pos[i]; // 位置偏差
// 检查调整条件：有偏差且相邻标记点之间有足够空间
if ((d >= 1 && n_pos[i + 1] - n_pos[i] > 1) || (d <= -1 && n_pos[i - 1] - n_pos[i] < -1)) {
int dInt = (d > 0) ? 1 : -1;
double qs = Parabolic(i, dInt); // 首选抛物线插值
// 检查插值结果是否在合理范围内
if (q[i - 1] < qs && qs < q[i + 1])
q[i] = qs;
else
q[i] = Linear(i, dInt); // 次选线性插值
n_pos[i] += dInt; // 更新位置索引
}
}
/**
* 抛物线插值：使用二次函数进行平滑插值
* 公式：q[i] += d/(n_pos[i+1]-n_pos[i-1]) *
* [ (n_pos[i]-n_pos[i-1]+d)*(q[i+1]-q[i])/(n_pos[i+1]-n_pos[i]) +
* (n_pos[i+1]-n_pos[i]-d)*(q[i]-q[i-1])/(n_pos[i]-n_pos[i-1]) ]
*/
double Parabolic(int i, double d) const {
assert(i > 0 && i < markerCount - 1);
return q[i] + d / (n_pos[i + 1] - n_pos[i - 1]) *
((n_pos[i] - n_pos[i - 1] + d) * (q[i + 1] - q[i]) / (n_pos[i + 1] - n_pos[i]) +
(n_pos[i + 1] - n_pos[i] - d) * (q[i] - q[i - 1]) / (n_pos[i] - n_pos[i - 1]));
}
// 线性插值：简单的线性插值方法
double Linear(int i, int d) const {
if (i + d < 0 || i + d >= markerCount) return q[i];
return q[i] + d * (q[i + d] - q[i]) / (n_pos[i + d] - n_pos[i]);
}
/**
* 初始化标记点：对前markerCount个数据点进行排序和初始化
* 这是P²算法的关键初始化步骤
*/
void initializeIfNeeded() {
if (!initialized && totalCount >= markerCount) {
// 排序初始数据点
std::sort(q.begin(), q.end());
// 计算期望位置和实际位置
UpdateNs(markerCount - 1);
for (int i = 0; i < markerCount; i++) {
if (i < n_pos.size()) {
n_pos[i] = static_cast<int>(std::round(ns[i]));
}
}
// 根据位置索引设置初始分位数值
std::vector<double> temp(q);
for (int i = 0; i < markerCount; i++) {
if (i < temp.size() && n_pos[i] < temp.size()) {
q[i] = temp[n_pos[i]];
}
}
UpdateNs(markerCount - 1);
initialized = true;
}
}
/**
* 获取内部估计器的分位数值
* 处理不同数据量情况下的分位数计算
*/
double getInternalQuantile(double p) {
if (totalCount == 0) return NAN;
// 边界分位数直接返回
if (p == 0.0) return currentMin;
if (p == 1.0) return currentMax;
// 数据量不足时使用精确计算
if (totalCount <= markerCount) {
std::vector<double> sortedData(q.begin(), q.begin() + totalCount);
std::sort(sortedData.begin(), sortedData.end());
int index = static_cast<int>(std::round(p * static_cast<double>(totalCount - 1)));
index = std::max(0, std::min(totalCount - 1, index));
return sortedData[index];
}
// 查找预定义概率的对应标记点
for (int i = 0; i < m; i++) {
if (std::abs(probabilities[i] - p) < 1e-10) {
int targetIndex = 2 * i + 2; // 主标记点索引
if (targetIndex < q.size()) return q[targetIndex];
}
}
// 边界概率插值
if (p <= probabilities[0]) {
double fraction = p / probabilities[0];
if (2 < q.size()) return q[0] + fraction * (q[2] - q[0]);
}
if (p >= probabilities[m - 1]) {
double fraction = (p - probabilities[m - 1]) / (1 - probabilities[m - 1]);
int lastIndex = 2 * m;
if (lastIndex < q.size() && lastIndex + 2 < q.size())
return q[lastIndex] + fraction * (q[lastIndex + 2] - q[lastIndex]);
}
// 中间概率插值
for (int i = 0; i < m - 1; i++) {
if (p >= probabilities[i] && p <= probabilities[i + 1]) {
double fraction = (p - probabilities[i]) / (probabilities[i + 1] - probabilities[i]);
int index1 = 2 * i + 2;
int index2 = 2 * i + 4;
if (index1 < q.size() && index2 < q.size()) return q[index1] + fraction * (q[index2] - q[index1]);
}
}
throw std::runtime_error("Failed to compute internal quantile");
}
// 重置内部估计器状态
void resetInternalEstimator() {
std::fill(n_pos.begin(), n_pos.end(), 0);
std::fill(ns.begin(), ns.end(), 0.0);
std::fill(q.begin(), q.end(), 0.0);
totalCount = 0;
initialized = false;
windowMin = std::numeric_limits<double>::max();
windowMax = std::numeric_limits<double>::lowest();
}
};
