В этом примере показано, как соответствовать данным о хвосте к Обобщенному распределению Парето по оценке наибольшего правдоподобия.
Подбор кривой параметрическому распределению к данным иногда приводит к модели, которая соглашается хорошо с данными в областях высокой плотности, но плохо в областях низкой плотности. Для одномодовых распределений, таких как t нормального или Студента, эти низкие области плотности известны как "хвосты" распределения. Одна причина, почему сила модели соответствует плохо в хвостах, состоит в том, что по определению, существует меньше данных в хвостах, на которых можно основывать выбор модели, и таким образом, модели часто выбираются на основе их способности соответствовать данным около режима. Другая причина может состоять в том, что распределение действительных данных часто более сложно, чем обычные параметрические модели.
Однако во многих приложениях, подгонка данных в хвосте является основным беспокойством. Обобщенное распределение Парето (GP) было разработано как распределение, которое может смоделировать хвосты большого разнообразия распределений, на основе теоретических аргументов. Один подход к распределению, соответствующему, который включает GP, должен использовать непараметрическую подгонку (эмпирическая кумулятивная функция распределения, например) в областях, где существует много наблюдений, и соответствовать GP к хвосту (хвостам) данных.
Обобщенный Парето (GP) является скошенным правом распределением, параметризованным параметром формы, k, и масштабным коэффициентом, сигмой. k также известен как "параметр" индекса хвоста, и может быть положителен, нуль, или отрицателен.
x = linspace(0,10,1000); plot(x,gppdf(x,-.4,1),'-', x,gppdf(x,0,1),'-', x,gppdf(x,2,1),'-'); xlabel('x / sigma'); ylabel('Probability density'); legend({'k < 0' 'k = 0' 'k > 0'});
Заметьте, что для k <0, GP имеет нулевую вероятность выше верхнего предела - (1/К). Для k> = 0, GP не имеет никакого верхнего предела. Кроме того, GP часто используется в сочетании с одной третью, пороговый параметр, который переключает нижний предел далеко от нуля. Нам не будет нужна та общность здесь.
Распределение GP является обобщением обоих экспоненциальное распределение (k = 0) и распределение Парето (k> 0). GP включает те два распределения в более многочисленное семейство так, чтобы непрерывная область значений форм была возможна.
Распределение GP может быть задано конструктивно в терминах exceedances. Начиная с вероятностного распределения, правый хвост которого понижается, чтобы обнулить, такой как нормальное, мы можем произвести случайные значения независимо от того распределения. Если мы фиксируем пороговое значение, выводим все значения, которые являются ниже порога и вычитают порог прочь значений, которые не выводятся, результат известен exceedances. Распределением exceedances является приблизительно GP. Точно так же мы можем установить порог в левом хвосте распределения и проигнорировать все значения выше того порога. Порог должен далеко достаточно отсутствовать в хвосте исходного распределения для приближения, чтобы быть разумным.
Исходное распределение определяет параметр формы, k, получившегося распределения GP. Распределения, хвосты которых уменьшаются как полином, такой как t Студента, вывод к положительному параметру формы. Распределения, хвосты которых уменьшаются экспоненциально, такой как нормальное, соответствуют нулевому параметру формы. Распределения с конечными хвостами, такими как бета, соответствуют отрицательному параметру формы.
Реальные приложения для распределения GP включают экстремальные значения моделирования фондового рынка, возвращается, и моделирование экстремальных лавинных рассылок. В данном примере мы будем использовать симулированные данные, сгенерированные от t распределения Студента с 5 степенями свободы. Мы возьмем самые большие 5% из 2 000 наблюдений от t распределения, и затем вычтем от 95%-го квантиля, чтобы получить exceedances.
rng(3,'twister');
x = trnd(5,2000,1);
q = quantile(x,.95);
y = x(x>q) - q;
n = numel(y)
n = 100
Распределение GP задано для 0 <сигма и-Inf <k <Inf. Однако интерпретация результатов оценки наибольшего правдоподобия проблематична когда k <-1/2. К счастью, те случаи соответствуют подходящим хвостам от распределений как бета или треугольный, и так не представят проблему здесь.
paramEsts = gpfit(y); kHat = paramEsts(1) % Tail index parameter sigmaHat = paramEsts(2) % Scale parameter
kHat =
0.0987
sigmaHat =
0.7156
Как может ожидаться, поскольку симулированные данные были сгенерированы с помощью t распределения, оценка k положительна.
Чтобы визуально оценить, насколько хороший подгонка, мы построим масштабированную гистограмму данных о хвосте, наложенных с функцией плотности GP, который мы оценили. Гистограмма масштабируется так, чтобы времена высот панели их сумма ширины к 1.
bins = 0:.25:7; h = bar(bins,histc(y,bins)/(length(y)*.25),'histc'); h.FaceColor = [.9 .9 .9]; ygrid = linspace(0,1.1*max(y),100); line(ygrid,gppdf(ygrid,kHat,sigmaHat)); xlim([0,6]); xlabel('Exceedance'); ylabel('Probability Density');
Мы использовали довольно маленькую ширину интервала, таким образом, существует много шума в гистограмме. Несмотря на это, подходящая плотность следует за формой данных, и таким образом, модель GP, кажется, хороший выбор.
Мы можем также сравнить эмпирический CDF с подходящим CDF.
[F,yi] = ecdf(y); plot(yi,gpcdf(yi,kHat,sigmaHat),'-'); hold on; stairs(yi,F,'r'); hold off; legend('Fitted Generalized Pareto CDF','Empirical CDF','location','southeast');
Чтобы определить количество точности оценок, мы будем использовать стандартные погрешности, вычисленные из асимптотической ковариационной матрицы средств оценки наибольшего правдоподобия. Функциональный gplike вычисляет, как его второй выход, числовое приближение к той ковариационной матрице. В качестве альтернативы мы, возможно, вызвали gpfit с двумя выходными аргументами, и это возвратило бы доверительные интервалы для параметров.
[nll,acov] = gplike(paramEsts, y); stdErr = sqrt(diag(acov))
stdErr =
0.1158
0.1093
Эти стандартные погрешности указывают, что относительная точность оценки для k вполне немного ниже, чем это для сигмы - ее стандартная погрешность находится на порядке самой оценки. Параметры формы часто затрудняют оценку. Важно иметь в виду, что расчет этих стандартных погрешностей принял, что модель GP правильна, и что у нас есть достаточно данных для асимптотического приближения к ковариационной матрице, чтобы содержать.
Интерпретация стандартных погрешностей обычно включает предположение, что, если та же подгонка могла бы много раз повторяться на данных, которые прибыли из того же источника, оценки наибольшего правдоподобия параметров будут приблизительно следовать за нормальным распределением. Например, доверительные интервалы часто базируются это предположение.
Однако то нормальное приближение может или не может быть хорошим. Чтобы оценить, насколько хороший это находится в этом примере, мы можем использовать симуляцию начальной загрузки. Мы сгенерируем 1000, реплицируют наборы данных путем передискретизации из данных, соответствуют распределению GP к каждому и сохраняют все реплицировать оценки.
replEsts = bootstrp(1000,@gpfit,y);
Как грубое проверяют распределение выборки средств оценки параметра, мы можем посмотреть на гистограммы начальной загрузки, реплицирует.
subplot(2,1,1); hist(replEsts(:,1)); title('Bootstrap estimates of k'); subplot(2,1,2); hist(replEsts(:,2)); title('Bootstrap estimates of sigma');
Гистограмма оценок начальной загрузки для k, кажется, только немного асимметрична, в то время как это для оценок сигмы определенно кажется скошенным направо. Общее средство от той скошенности должно оценить параметр и его стандартную погрешность на логарифмической шкале, где нормальное приближение может быть более разумным. График Q-Q является лучшим способом оценить нормальность, чем гистограмма, потому что ненормальность обнаруживается как точки, которые не приблизительно следуют за прямой линией. Давайте проверять, что видеть, преобразовывает ли журнал для сигмы, является соответствующим.
subplot(1,2,1); qqplot(replEsts(:,1)); title('Bootstrap estimates of k'); subplot(1,2,2); qqplot(log(replEsts(:,2))); title('Bootstrap estimates of log(sigma)');
Оценки начальной загрузки для k и журнала (сигма) появляются приемлемо близко к нормальности. График Q-Q для оценок сигмы, по нерегистрируемой шкале, подтвердил бы скошенность, которую мы уже видели в гистограмме. Таким образом было бы более разумно создать доверительный интервал для сигмы первым вычислением того для журнала (сигма) под предположением о нормальности, и затем возведение в степень, чтобы преобразовать тот интервал назад к исходной шкале для сигмы.
На самом деле это точно что функциональный gpfit делает негласно.
[paramEsts,paramCI] = gpfit(y);
kHat kCI = paramCI(:,1)
kHat =
0.0987
kCI =
-0.1283
0.3258
sigmaHat sigmaCI = paramCI(:,2)
sigmaHat =
0.7156
sigmaCI =
0.5305
0.9654
Заметьте, что, в то время как 95%-й доверительный интервал для k симметричен об оценке наибольшего правдоподобия, доверительный интервал для сигмы не. Поэтому это было создано путем преобразования симметричного CI для журнала (сигма).