Документация

Обобщенное распределение Парето

Обобщенный Парето (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 включает те два дистрибутива в более многочисленное семейство так, чтобы непрерывная область значений форм была возможна.

Симуляция данные Exceedance

Распределение 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 положительна.

Проверка подгонки визуально

Чтобы визуально оценить, насколько хороший подгонка, мы построим масштабированную гистограмму данных о хвосте, overlayed с функцией плотности 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 для журнала (сигма).