Итеративная квантильная оценка в Matlab

Я пытаюсь реализовать интерактивный алгоритм для оценки квантилей в данных, созданных в результате моделирования Монте-Карло. Я хочу сделать его итеративным, потому что у меня есть много итераций и переменных, поэтому сохранение всех точек данных и использование функции Matlab quantile займет большую часть памяти, которая мне действительно нужна для моделирования.

Я нашел несколько подходов на основе -W/full" rel="nofollow">процесс Роббина-Монро, заданный

Формула

Реализация с управляющей последовательностью ct = c / t, где c — константа, довольно проста. В цитируемой работе показано, что c = 2 * sqrt(2 * pi) дает неплохие результаты, по крайней мере, для медианы. Но они также предлагают адаптивный подход, основанный на оценке гистограммы. К сожалению, я пока не придумал, как реализовать эту адаптацию.

Я протестировал реализацию с константой c для трех тестовых образцов с 10 000 точек данных. Значение 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

Есть ли более надежное решение для итеративной квантильной оценки или у кого-нибудь есть реализация адаптивного подхода с небольшим потреблением памяти?


matlab iteration estimation quantile
person JotWe    schedule 27.10.2016    source источник
comment
Несколько потенциальных проблем: indices представляет собой массив 1 и 0, не уверен, что должен делать probabilities(indices). Кроме того, я думаю, что вам нужно что-то вроде quantile(index) = quantile(index-1) + control * indicator;. Наконец, я думаю, что вы неправильно реализовали c/t. Я бы подумал, что t - это время, если только экземпляры между вашими точками данных не равны 1сек.   -  person mpaskov    schedule 27.10.2016
comment
Спасибо за комментарий. На мой взгляд, индекс t просто обозначает счетчик итераций, поэтому время не требуется. Переменная quantile представляет собой вектор того же размера, что и вероятности, в данном случае 1x3, содержащий итеративные квантильные оценки для probabilities = [0.1, 0.4, 0.7]. Последняя строка цикла for обновляет эти оценки. Конструкция индексов/индикаторов — это моя реализация индикаторной функции I, которая выбирает, когда использовать probabilities или probabilities - 1.   -  person JotWe    schedule 28.10.2016

Ответы (1)


arrow_upward
0
arrow_downward

Попробовав несколько адаптивных итерационных подходов, которые я нашел в литературе, но без особого успеха (не уверен, правильно ли я все сделал), я пришел к решению, которое дает хорошие результаты для моих тестовых образцов, а также для реального метода Монте-Карло. -Моделирование.

Я буферизирую подмножество результатов моделирования, вычисляю квантили выборки и в конце усредняю ​​по всем квантилям выборки подмножества. Кажется, это работает довольно хорошо и без настройки многих параметров. Единственным параметром является размер буфера, который в моем случае равен 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);
person JotWe    schedule 28.10.2016
comment
Я только что нашел этот пост на Cross Validated. Кажется, это связано. - person JotWe; 28.10.2016
comment
Имеет смысл, что quantileEstimation3 будет быстро сходиться или насыщаться, поскольку фактор изменения control быстро падает, поэтому новые изменения адаптируются все медленнее. - person mpaskov; 28.10.2016
comment
Да, но даже приспособления управляющего фактора, которые я нашел и испробовал, не улучшили для меня результатов. - person JotWe; 28.10.2016