Matlab 中的迭代分位数估计

Question

我正在尝试实现一种交互式算法来估计从蒙特卡洛模拟生成的数据中的分位数。我想让它迭代，因为我有很多迭代和变量，所以存储所有数据点和使用 Matlab 的 quantile 函数会占用我模拟实际需要的大部分内存。

我找到了一些基于 Robbin-Monro process 的方法，由

给出

控制序列 c_t = c / t 的实现非常简单，其中 c 是常量。在引用的论文中，他们表明 c = 2 * sqrt(2 * pi) 给出了相当好的结果，至少对于中位数而言。但他们也提出了一种基于直方图估计的自适应方法。不幸的是，我还没有弄清楚如何实现这种适配。

我用 10.000 个数据点测试了三个测试样本的 implementation with a constant c。值 c = 2 * sqrt(2 * pi) 对我来说效果不佳，但 c = 100 对于测试样本来说看起来相当不错。然而，这个选择不是很稳健，在实际的蒙特卡洛模拟中失败了，给出的结果大相径庭。

probabilities = [0.1, 0.4, 0.7];
controlFactor = 100;
quantile = zeros(size(probabilities));
indicator = zeros(size(probabilities));
for index = 1:length(data)
    control = controlFactor / index;
    indices = (data(index) >= quantile);
    indicator(indices) = probabilities(indices);
    indices = (data(index) < quantile);
    indicator(indices) = probabilities(indices) - 1;
    quantile = quantile + control * indicator;
end

是否有更稳健的迭代分位数估计解决方案，或者是否有人实现了内存消耗小的自适应方法？

Answer 1

在尝试了一些我在文献中发现的自适应迭代方法但没有取得巨大成功之后（不确定，如果我做对了），我想出了一个解决方案，它为我的测试样本和实际的蒙特卡洛模拟。

我缓冲了模拟结果的一个子集，计算样本分位数，最后对所有子集样本分位数进行平均。这似乎工作得很好并且无需调整许多参数。唯一的参数是缓冲区大小，在我的例子中是 100。

结果收敛得非常快，增加样本量不会显着改善结果。似乎有一个小但恒定的偏差，大概是子集样本分位数的平均误差。这就是我的解决方案的缺点。通过选择缓冲区大小，可以确定可实现的精度。增加缓冲区大小可以减少这种偏差。最后，这似乎是一个内存和准确性的权衡。

% Generate data
rng('default');
data = sqrt(0.5) * randn(10000, 1) + 5 * rand(10000, 1) + 10;

% Set parameters
probabilities = 0.2;

% Compute reference sample quantiles
quantileEstimation1 = quantile(data, probabilities);

% Estimate quantiles with computing the mean over a number of subset
% sample quantiles
subsetSize = 100;
quantileSum = 0;
for index = 1:length(data) / subsetSize;

    quantileSum = quantileSum + quantile(data(((index - 1) * subsetSize + 1):(index * subsetSize)), probabilities);

end
quantileEstimation2 = quantileSum / (length(data) / subsetSize);

% Estimate quantiles with iterative computation
quantileEstimation3 = zeros(size(probabilities));
indicator = zeros(size(probabilities));
controlFactor = 2 * sqrt(2 * pi);
for index = 1:length(data)

    control = controlFactor / index;
    indices = (data(index) >= quantileEstimation3);
    indicator(indices) = probabilities(indices);
    indices = (data(index) < quantileEstimation3);
    indicator(indices) = probabilities(indices) - 1;
    quantileEstimation3 = quantileEstimation3 + control * indicator;

end

fprintf('Reference result: %f\nSubset result: %f\nIterative result: %f\n\n', quantileEstimation1, quantileEstimation2, quantileEstimation3);

Matlab 中的迭代分位数估计

Iterative quantile estimation in Matlab

iteration

matlab

estimation

quantile