Этот пример показывает, как соответствовать обобщенному распределению экстремума с помощью оценки наибольшего правдоподобия. Распределение экстремума используется, чтобы смоделировать самое большое или наименьшее значение от группы или блока данных.
Три типа дистрибутивов экстремума распространены, каждый как ограничивающий случай для различных типов базовых дистрибутивов. Например, экстремум типа I является предельным распределением максимума (или минимум) блока нормально распределенных данных, когда размер блока становится большим. В этом примере мы проиллюстрируем, как соответствовать таким данным с помощью одного распределения, которое включает все три типа дистрибутивов экстремума как особый случай, и привлеките основанные на вероятности доверительные интервалы по делу о квантилях подходящего распределения.
Распределение Обобщенного экстремума (GEV) объединяет тип I, тип II, и дистрибутивы экстремума типа III в одно семейство, чтобы позволить непрерывную область значений возможных форм. Это параметризовано с местоположением и масштабными коэффициентами, mu и сигмой и параметром формы, 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 приближается 0, GEV неограничен. Это может быть получено в итоге как ограничение, которое 1+k* (y-mu) / сигма должно быть положительно.
GEV может быть задан конструктивно как ограничивающее распределение максимумов блока (или минимумы). Таким образом, если вы генерируете большое количество независимых случайных значений от одного распределения вероятностей и принимаете их максимальное значение, распределением того максимума является приблизительно GEV.
Исходное распределение определяет параметр формы, k, получившегося распределения GEV. Дистрибутивы, хвосты которых уменьшаются как полином, такой как t Студента, вывод к положительному параметру формы. Дистрибутивы, хвосты которых уменьшаются экспоненциально, такой как нормальное, соответствуют нулевому параметру формы. Дистрибутивы с конечными хвостами, такими как бета, соответствуют отрицательному параметру формы.
Действительные приложения для GEV могут включать моделирование самого большого возврата для запаса в течение каждого месяца. Здесь, мы моделируем данные путем взятия максимума 25 значений от t распределения Студента с двумя степенями свободы. Моделируемые данные будут включать 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;
Во-первых, мы построим масштабированную гистограмму данных, overlayed с 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 часто является количеством интереса к анализу данных о максимумах блока.
Например, Комната уровня возврата задана, когда максимальное значение блока ожидало быть превышенным только однажды в блоках m. Это только (1-1/m) 'th квантиль. Мы можем включить оценки параметра наибольшего правдоподобия в обратный CDF, чтобы оценить Комнату для m=10.
R10MLE = gevinv(1-1./10,kMLE,sigmaMLE,muMLE)
R10MLE =
9.0724
Мы могли вычислить пределы достоверности для R10 с помощью асимптотических приближений, но те не могут быть допустимыми. Вместо этого мы будем использовать основанный на вероятности метод, чтобы вычислить пределы достоверности. Этот метод часто приводит к более точным результатам, чем один на основе предполагаемой ковариационной матрицы оценок параметра.
Учитывая любое множество значений для параметров mu, сигмы и k, мы можем вычислить логарифмическую вероятность - например, MLEs являются значениями параметров, которые максимизируют логарифмическую вероятность GEV. Когда значения параметров переезжают от MLEs, их логарифмическая вероятность обычно становится значительно меньше, чем максимум. Если мы смотрим на набор значений параметров, которые производят логарифмическую вероятность, больше, чем заданное критическое значение, это - сложная область в пространстве параметров. Однако для подходящего критического значения, это - область уверенности для параметров модели. Область содержит значения параметров, которые "совместимы с данными". Критическое значение, которое определяет область, основано на приближении хи-квадрата, и мы будем использовать 95% в качестве нашего доверительного уровня. (Обратите внимание на то, что мы будем на самом деле работать с отрицанием логарифмической вероятности.)
nllCritVal = gevlike([kMLE,sigmaMLE,muMLE],y) + .5*chi2inv(.95,1)
nllCritVal = 170.3044
Для любого набора значений параметров mu, сигмы и k, мы можем вычислить R10. Поэтому мы можем найти наименьшее значение R10 достигнутым в критической области пространства параметров, где отрицательная логарифмическая вероятность больше, чем критическое значение. То наименьшее значение является более низким основанным на вероятности пределом достоверности для R10.
Это трудно визуализировать во всех трех размерностях параметра, но как мысленный эксперимент, мы можем зафиксировать параметр формы, k, мы видим, как процедура переработала бы два остающихся параметра, сигму и 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, Комната является линейной функцией сигмы и mu. Полужирные красные контуры являются самыми низкими и самыми высокими значениями R10, которые находятся в пределах критической области. В полном трехмерном пространстве параметров контуры логарифмической вероятности были бы эллипсоидальными, и контуры R10 будут поверхностями.
Нахождение более низкого предела достоверности для R10 является задачей оптимизации с нелинейными ограничениями неравенства, и таким образом, мы будем использовать функциональный fmincon от Optimization Toolbox™. Мы должны найти наименьшее значение R10, и поэтому цель, которая будет минимизирована, является самим R10, равный обратному CDF, оцененному для p=1-1/m. Мы создадим функцию обертки, которая вычисляет Комнату специально для 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);
Чтобы найти верхний предел достоверности вероятности для 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 и переработали два остающихся параметра, сигму и 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');
