В этом примере показано, как соответствовать обобщенному распределению экстремума с помощью оценки наибольшего правдоподобия. Распределение экстремума используется, чтобы смоделировать самое большое или наименьшее значение от группы или блока данных.
Три типа распределений экстремума распространены, каждый как ограничивающий случай для различных типов базовых распределений. Например, экстремум типа 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;
Во-первых, мы построим масштабированную гистограмму данных, наложенных с 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');
