Обобщенное распределение экстремальных значений

Распределение Обобщенных Экстремальных Значений (GEV) объединяет распределения экстремальных значений типа I, типа II и типа III в одно семейство, чтобы позволить непрерывная область значений возможных форм. Он параметризован с параметрами местоположения и шкалы, mu и sigma, и параметром формы, k. Когда k < 0, GEV эквивалентно экстремальному значению типа III. Когда k > 0, GEV эквивалентно типу II. В пределе при приближении k к 0 GEV становится типом I.

x = linspace(-3,6,1000);
plot(x,gevpdf(x,-.5,1,0),'-', x,gevpdf(x,0,1,0),'-', x,gevpdf(x,.5,1,0),'-');
xlabel('(x-mu) / sigma');
ylabel('Probability Density');
legend({'k < 0, Type III' 'k = 0, Type I' 'k > 0, Type II'});

Заметьте, что для k < 0 или k > 0 плотность имеет нулевую вероятность выше или ниже, соответственно, верхнюю или нижнюю границу - (1/k). В пределе при приближении k к 0 GEV неограниченен. Это может быть суммировано как ограничение, что 1 + k * (y-mu )/sigma должны быть положительными.

Симуляция максимальных данных блоков

GEV может быть определено конструктивно как ограничение распределения блока maxima (или minima). То есть, если вы генерируете большое количество независимых случайных значений из одного распределения вероятностей и принимаете их максимальное значение, распределение этого максимума составляет приблизительно GEV.

Исходное распределение определяет параметр k формы полученного распределения GEV. Распределения, чьи хвосты падают как полином, такой как t Студента, приводят к положительному параметру формы. Распределения, чьи хвосты уменьшаются экспоненциально, такие как normal, соответствуют параметру нулевой формы. Распределения с конечными хвостами, такими как бета, соответствуют отрицательному параметру формы.

Реальные приложения для GEV могут включать моделирование наибольшего возврата для запаса в течение каждого месяца. Здесь мы будем симулировать данные, взяв максимум 25 значений из распределение Студента с двумя степенями свободы. Моделируемые данные будут включать 75 максимальных значений случайных блоков.

rng(0,'twister');
y = max(trnd(2,25,75),[],1);

Подбор распределения по максимальному правдоподобию

Функция gevfit возвращает обе оценки максимальных параметров правдоподобия и (по умолчанию) 95% доверительных интервалов.

[paramEsts,paramCIs] = gevfit(y);

kMLE = paramEsts(1)        % Shape parameter
sigmaMLE = paramEsts(2)    % Scale parameter
muMLE = paramEsts(3)       % Location parameter

kMLE =

    0.4901


sigmaMLE =

    1.4856


muMLE =

    2.9710

kCI = paramCIs(:,1)
sigmaCI = paramCIs(:,2)
muCI = paramCIs(:,3)

kCI =

    0.2020
    0.7782


sigmaCI =

    1.1431
    1.9307


muCI =

    2.5599
    3.3821

Заметьте, что 95% доверительный интервал для k не включает значение нуль. Распределение экстремальных значений типа I, по-видимому, не является хорошей моделью для этих данных. Это имеет смысл, потому что базовое распределение для симуляции имело намного более тяжелые хвосты, чем нормальное, и распределение экстремальных значений типа II теоретически является правильным, когда размер блока становится большим.

В качестве альтернативы доверительным интервалам мы также можем вычислить приближение к асимптотической ковариационной матрице оценок параметров, и из них извлечь стандартные ошибки параметров.

[nll,acov] = gevlike(paramEsts,y);
paramSEs = sqrt(diag(acov))

paramSEs =

    0.1470
    0.1986
    0.2097

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

Чтобы визуально оценить, насколько хороша подгонка, мы рассмотрим графики установленной функции плотности вероятностей (PDF) и совокупной функции распределения (CDF).

Поддержка GEV зависит от значений параметров. В этом случае оценка для k положительна, поэтому установленное распределение имеет нулевую вероятность ниже нижней границы.

lowerBnd = muMLE-sigmaMLE./kMLE;

Во-первых, мы построим масштабированную гистограмму данных, наложенную на PDF для подобранной модели GEV. Эта гистограмма масштабируется так, чтобы высоты штриха умножали их ширину на 1, чтобы сделать ее сопоставимой с PDF.

ymax = 1.1*max(y);
bins = floor(lowerBnd):ceil(ymax);
h = bar(bins,histc(y,bins)/length(y),'histc');
h.FaceColor = [.9 .9 .9];
ygrid = linspace(lowerBnd,ymax,100);
line(ygrid,gevpdf(ygrid,kMLE,sigmaMLE,muMLE));
xlabel('Block Maximum');
ylabel('Probability Density');
xlim([lowerBnd ymax]);

Мы также можем сравнить подгонку к данным с точки зрения совокупной вероятности, путем наложения эмпирического CDF и подобранного CDF.

[F,yi] = ecdf(y);
plot(ygrid,gevcdf(ygrid,kMLE,sigmaMLE,muMLE),'-');
hold on;
stairs(yi,F,'r');
hold off;
xlabel('Block Maximum');
ylabel('Cumulative Probability');
legend('Fitted Generalized Extreme Value CDF','Empirical CDF','location','southeast');
xlim([lowerBnd ymax]);

Оценка квантилей модели

Хотя оценки параметров могут быть важны сами по себе, количество подобранной модели GEV часто является количеством интереса к анализу данных максимумов блоков.

Для примера уровень Rm возврата определяется как максимальное значение блока, которое, как ожидается, будет превышено только один раз в m блоках. Это всего лишь 1-1/m 'й квантиль. Мы можем включить оценки максимального параметра правдоподобия в обратный CDF, чтобы оценить Rm для m = 10.

R10MLE = gevinv(1-1./10,kMLE,sigmaMLE,muMLE)

R10MLE =

    9.0724

Мы могли бы вычислить доверие пределы для R10, используя асимптотические приближения, но они могут быть недопустимыми. Вместо этого мы будем использовать метод, основанный на вероятностях, для вычисления доверительных пределов. Этот метод часто дает более точные результаты, чем один, на основе предполагаемой ковариационной матрицы оценок параметров.

Учитывая любое множество значений для параметров mu, sigma и k, мы можем вычислить логарифмическую правдоподобность - например, MLE являются значениями параметров, которые максимизируют логарифмическую правдоподобность GEV. Когда значения параметров отходят от MLE, их логарифмическая правдоподобность обычно становится значительно меньше максимальной. Если мы рассмотрим набор значений параметров, которые генерируют логарифмическую правдоподобность, большую заданного критического значения, это сложная область в пространстве параметров. Однако для подходящего критического значения это область доверия для параметров модели. Область содержит значения параметров, которые «совместимы с данными». Критическое значение, которое определяет область, основано на хи-квадратном приближении, и мы будем использовать 95% в качестве нашего доверительного уровня. (Обратите внимание, что мы на самом деле будем работать с отрицательным логарифмической логарифмической правдоподобностью.)

nllCritVal = gevlike([kMLE,sigmaMLE,muMLE],y) + .5*chi2inv(.95,1)

nllCritVal =

  170.3044

Для любого набора значений параметров mu, sigma и k, мы можем вычислить R10. Поэтому мы можем найти наименьшее R10 значение, достигнутое в критической области пространства параметров, где отрицательная логарифмическая правдоподобность больше критического значения. Это наименьшее значение является нижним основанным на вероятностях пределом доверия для R10.

Это трудно визуализировать во всех трёх размерностях параметра, но в качестве мысленного эксперимента мы можем зафиксировать параметр формы, k, мы можем увидеть, как процедура работала бы по двум оставшимся параметрам, sigma и mu.

sigmaGrid = linspace(.8, 2.25, 110);
muGrid = linspace(2.4, 3.6);
nllGrid = zeros(length(sigmaGrid),length(muGrid));
R10Grid = zeros(length(sigmaGrid),length(muGrid));
for i = 1:size(nllGrid,1)
    for j = 1:size(nllGrid,2)
        nllGrid(i,j) = gevlike([kMLE,sigmaGrid(i),muGrid(j)],y);
        R10Grid(i,j) = gevinv(1-1./10,kMLE,sigmaGrid(i),muGrid(j));
    end
end
nllGrid(nllGrid>gevlike([kMLE,sigmaMLE,muMLE],y)+6) = NaN;
contour(muGrid,sigmaGrid,R10Grid,6.14:.64:12.14,'LineColor','r');
hold on
contour(muGrid,sigmaGrid,R10Grid,[7.42 11.26],'LineWidth',2,'LineColor','r');
contour(muGrid,sigmaGrid,nllGrid,[168.7 169.1 169.6 170.3:1:173.3],'LineColor','b');
contour(muGrid,sigmaGrid,nllGrid,[nllCritVal nllCritVal],'LineWidth',2,'LineColor','b');
hold off
axis([2.4 3.6 .8 2.25]);
xlabel('mu');
ylabel('sigma');

Синие контуры представляют поверхность логарифмической правдоподобности, а жирный синий контур является границей критической области. Красные контуры представляют поверхность для R10 - большие значения верхней части справа, ниже внизу слева. Контуры прямых линий, потому что для фиксированного k, Rm является линейной функцией sigma и mu. Жирные красные контуры являются самыми низкими и самыми высокими значениями R10, которые попадают в критическую область. В полном трехмерном пространстве параметров логарифмической правдоподобности контуры будут эллипсоидальными, а R10 контуры будут поверхностями.

Нахождение нижнего доверительного предела для R10 является задачей оптимизации с нелинейными ограничениями неравенства, и поэтому мы будем использовать функцию fmincon из Optimization Toolbox™. Нам нужно найти наименьшее R10 значение, и поэтому цель, которая будет минимизирована, сама R10 себе, равна обратной CDF, оцененной для p = 1-1/m. Мы создадим функцию обертки, которая вычисляет Rm специально для m = 10.

CIobjfun = @(params) gevinv(1-1./10,params(1),params(2),params(3));

Чтобы выполнить ограниченную оптимизацию, нам также понадобится функция, которая задает ограничение, то есть что отрицательная логарифмическая правдоподобность меньше критического значения. Функция ограничения должна возвращать положительные значения, когда ограничение нарушено. Создадим анонимную функцию, используя моделируемые данные и критическое значение логарифмической правдоподобности. Это также возвращает пустое значение, потому что мы не используем здесь никаких ограничений равенства.

CIconfun = @(params) deal(gevlike(params,y) - nllCritVal, []);

Наконец, мы звоним fmincon, используя алгоритм активного набора, чтобы выполнить ограниченную оптимизацию.

opts = optimset('Algorithm','active-set', 'Display','notify', 'MaxFunEvals',500, ...
                'RelLineSrchBnd',.1, 'RelLineSrchBndDuration',Inf);
[params,R10Lower,flag,output] = ...
    fmincon(CIobjfun,paramEsts,[],[],[],[],[],[],CIconfun,opts);

Feasible point with lower objective function value found.

Чтобы найти верхний предел доверия правдоподобия для R10, мы просто обращаем знак на целевой функции, чтобы найти самое большое R10 значение в критической области и вызываем fmincon второй раз.

CIobjfun = @(params) -gevinv(1-1./10,params(1),params(2),params(3));
[params,R10Upper,flag,output] = ...
    fmincon(CIobjfun,paramEsts,[],[],[],[],[],[],CIconfun,opts);
R10Upper = -R10Upper;

R10CI = [R10Lower, R10Upper]

R10CI =

    7.0841   13.4452

plot(ygrid,gevcdf(ygrid,kMLE,sigmaMLE,muMLE),'-');
hold on;
stairs(yi,F,'r');
plot(R10CI([1 1 1 1 2 2 2 2]), [.88 .92 NaN .9 .9 NaN .88 .92],'k-')
hold off;
xlabel('Block Maximum');
ylabel('Cumulative Probability');
legend('Fitted Generalized Extreme Value CDF','Empirical CDF', ...
       'R_{10} 95% CI','location','southeast');
xlim([lowerBnd ymax]);

Профиль правдоподобия для квантиля

Иногда только интервал не дает достаточной информации об оценаемом количестве, и вместо этого необходима вероятность профиля. Чтобы найти профиль логарифмической правдоподобности для R10, мы исправим возможное значение для R10, а затем максимизируем логарифмическую вероятность GEV с параметрами, ограниченными, чтобы они соответствовали этому текущему значению R10. Это нелинейное ограничение равенства. Если мы сделаем это в области значений R10 значений, мы получим профиль вероятности.

Как и в случае доверительного интервала, основанного на вероятностях, мы можем думать о том, какой была бы эта процедура, если бы мы фиксировали k и работали над двумя оставшимися параметрами, sigma и mu. Каждая красная контурная линия на контурном графике, показанном ранее, представляет фиксированное значение R10; оптимизация вероятности профиля состоит из шага вдоль одной R10 контурной линии, чтобы найти самый высокий логарифмический (синий) контур.

В этом примере мы вычислим вероятность профиля для R10 над значениями, которые были включены в доверительный интервал вероятности.

R10grid = linspace(R10CI(1)-.05*diff(R10CI), R10CI(2)+.05*diff(R10CI), 51);

Целевой функцией для оптимизации правдоподобия профиля является просто логарифмическая правдоподобность, используя моделируемые данные.

PLobjfun = @(params) gevlike(params,y);

Как использовать fmincon, нам понадобится функция, которая возвращает ненулевые значения, когда ограничение нарушено, то есть когда параметры не соответствуют текущему значению R10. Для каждого значения R10 мы создадим анонимную функцию для конкретного значения рассматриваемых R10. Это также возвращает пустое значение, потому что мы не используем здесь никаких ограничений неравенства.

Наконец-то мы позвоним fmincon при каждом значении R10, чтобы найти соответствующий ограниченный максимум логарифмической правдоподобности. Мы начнем с максимальной оценки вероятности R10 и поработаем в обоих направлениях.

Lprof = nan(size(R10grid));
params = paramEsts;
[dum,peak] = min(abs(R10grid-R10MLE));
for i = peak:1:length(R10grid)
    PLconfun = ...
        @(params) deal([], gevinv(1-1./10,params(1),params(2),params(3)) - R10grid(i));
    [params,Lprof(i),flag,output] = ...
        fmincon(PLobjfun,params,[],[],[],[],[],[],PLconfun,opts);
end
params = paramEsts;
for i = peak-1:-1:1
    PLconfun = ...
        @(params) deal([], gevinv(1-1./10,params(1),params(2),params(3)) - R10grid(i));
    [params,Lprof(i),flag,output] = ...
        fmincon(PLobjfun,params,[],[],[],[],[],[],PLconfun,opts);
end

plot(R10grid,-Lprof,'-', R10MLE,-gevlike(paramEsts,y),'ro', ...
     [R10grid(1), R10grid(end)],[-nllCritVal,-nllCritVal],'k--');
xlabel('R_{10}');
ylabel('Log-Likelihood');
legend('Profile likelihood','MLE','95% Conf. Limit');