Подбор экспоненциального распределения с помощью методом наименьших квадратов

Термин «наименьшие квадраты» чаще всего используется в контексте подбора кривой регрессионой линии или поверхности для моделирования переменной отклика как функции одной или нескольких переменных предиктора. Метод, описанный здесь, является очень другим приложением методом наименьших квадратов: одномерных подборов кривой распределения с одной переменной.

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

rng(37,'twister');
n = 100;
x = exprnd(2,n,1);

Затем вычислите эмпирическую совокупную функцию распределения (ECDF) данных. Это просто функция шага с переходом совокупной вероятности, p, 1/n в каждой точке данных, x.

x = sort(x);
p = ((1:n)-0.5)' ./ n;
stairs(x,p,'k-');
xlabel('x');
ylabel('Cumulative probability (p)');

Мы подберем экспоненциальное распределение к этим данным. Один из способов сделать это - найти экспоненциальное распределение, совокупная функция распределения (CDF) которого лучше всего аппроксимирует (в смысле, который будет объяснен ниже) ECDF данных. Экспоненциальный CDF является p = Pr {X < = x} = 1 - exp (-x/mu). Преобразование этого в -журнал (1-p) * mu = x дает линейное соотношение между -журнал (1-p) и x. Если данные происходят из экспоненциального, мы должны увидеть, по крайней мере, приблизительно, линейное соотношение, если мы подключаем вычисленные значения x и p из ECDF к этому уравнению. Если мы используем наименьшие квадраты для аппроксимации прямой линии через источник x против -журнал (1-p), то эта подобранная линия представляет экспоненциальное распределение, которое «ближе всего» к данным. Уклон линии является оценкой параметра mu.

Эквивалентно, мы можем думать о y = -журнал (1-p) как об «идеализированной выборке» из стандартного (среднего 1) экспоненциального распределения. Эти идеализированные значения равномерно разнесены по шкале вероятностей. Q-Q график x и y должен быть примерно линейным, если данные получены из экспоненциального распределения, и мы подгоняем линию наименьших квадратов через источник x против y.

y = -log(1 - p);
muHat = y \ x

muHat =

    1.8627

Постройте график данных и установленной линии.

plot(x,y,'+', y*muHat,y,'r--');
xlabel('x');
ylabel('y = -log(1-p)');

Заметьте, что линейная подгонка, которую мы сделали, минимизирует сумму квадратичных невязок в горизонтальном, или «x», направлении. Это потому, что значения для y = -журнал (1-p) детерминированны, и это значения x случайны. Также можно регрессировать y против x или использовать другие типы линейных подгонок, например, взвешенную регрессию, ортогональную регрессию или даже устойчивую регрессию. Мы не будем здесь изучать эти возможности.

Для сравнения подгонки данные по максимальной вероятности.

muMLE = expfit(x)

muMLE =

    1.7894

Теперь постройте график двух оцененных распределений по нетрансформированной совокупной шкале вероятностей.

stairs(x,p,'k-');
hold on
xgrid = linspace(0,1.1*max(x),100)';
plot(xgrid,expcdf(xgrid,muHat),'r--', xgrid,expcdf(xgrid,muMLE),'b--');
hold off
xlabel('x'); ylabel('Cumulative Probability (p)');
legend({'Data','LS Fit','ML Fit'},'location','southeast');

Эти два метода дают очень похожие подобранные распределения, хотя на аппроксимацию LS больше влияли наблюдения в хвосте распределения.

Подбор распределения Вейбула

Для немного более сложного примера моделируйте некоторые выборочные данные из распределения Вейбула и вычислите ECDF x.

n = 100;
x = wblrnd(2,1,n,1);
x = sort(x);
p = ((1:n)-0.5)' ./ n;

Чтобы подогнать распределение Вейбула к этим данным, заметьте, что CDF для Вейбула является p = Pr {X < = x} = 1 - exp (- (x/a) ^ b). Преобразование этого в log (a) + журнал (-журнал (1-p)) * (1/b) = log (x) снова дает линейную связь, на этот раз между журнал (-журнал (1-p)) и log (x). Мы можем использовать наименьшие квадраты, чтобы соответствовать прямой линии в преобразованной шкале, используя p и x из ECDF, и наклон и точка пересечения этой линии приводят к оценкам a и b.

logx = log(x);
logy = log(-log(1 - p));
poly = polyfit(logy,logx,1);
paramHat = [exp(poly(2)) 1/poly(1)]

paramHat =

    2.1420    1.0843

Постройте график данных и установленной линии на преобразованной шкале.

plot(logx,logy,'+', log(paramHat(1)) + logy/paramHat(2),logy,'r--');
xlabel('log(x)');
ylabel('log(-log(1-p))');

Для сравнения подгоните данные по максимальной вероятности и постройте два предполагаемых распределения по нетрансформированной шкале.

paramMLE = wblfit(x)
stairs(x,p,'k');
hold on
xgrid = linspace(0,1.1*max(x),100)';
plot(xgrid,wblcdf(xgrid,paramHat(1),paramHat(2)),'r--', ...
     xgrid,wblcdf(xgrid,paramMLE(1),paramMLE(2)),'b--');
hold off
xlabel('x'); ylabel('Cumulative Probability (p)');
legend({'Data','LS Fit','ML Fit'},'location','southeast');

paramMLE =

    2.1685    1.0372

Пример порогового параметра

Иногда необходимо подгонять положительные распределения, такие как Weibull или lognormal, с параметром порога. Например, случайная переменная Вейбула принимает значения над (0, Inf), а пороговый параметр, c, смещает этот диапазон на (c, Inf). Если пороговый параметр известен, то сложности нет. Но если пороговый параметр не известен, он должен быть вместо этого оценен. Эти модели трудно подогнать к максимальной вероятности - вероятность может иметь несколько режимов или даже стать бесконечной для значений параметров, которые не разумны для данных, и поэтому максимальная вероятность часто не является хорошим методом. Но с небольшим сложением к процедуре наименьших квадратов мы можем получить стабильные оценки.

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

n = 100;
x = wblrnd(4,2,n,1) + 4;
hist(x,20); xlim([0 16]);

Как мы можем подогнать трехпараметровое распределение Вейбула к этим данным? Если бы мы знали, что такое пороговое значение, например, 1, мы могли бы вычесть это значение из данных и затем использовать процедуру наименьших квадратов, чтобы оценить форму Вейбула и параметры шкалы.

x = sort(x);
p = ((1:n)-0.5)' ./ n;
logy = log(-log(1-p));
logxm1 = log(x-1);
poly1 = polyfit(log(-log(1-p)),log(x-1),1);
paramHat1 = [exp(poly1(2)) 1/poly1(1)]
plot(logxm1,logy,'b+', log(paramHat1(1)) + logy/paramHat1(2),logy,'r--');
xlabel('log(x-1)');
ylabel('log(-log(1-p))');

paramHat1 =

    7.4305    4.5574

Это не очень хорошая подгонка - журнал (x-1) и журнал (-log (1-p)) не имеют линейной связи. Конечно, это потому, что мы не знаем правильного порогового значения. Если мы попробуем вычесть различные пороговые значения, мы получим различные графики и различные оценки параметров.

logxm2 = log(x-2);
poly2 = polyfit(log(-log(1-p)),log(x-2),1);
paramHat2 = [exp(poly2(2)) 1/poly2(1)]

paramHat2 =

    6.4046    3.7690

logxm4 = log(x-4);
poly4 = polyfit(log(-log(1-p)),log(x-4),1);
paramHat4 = [exp(poly4(2)) 1/poly4(1)]

paramHat4 =

    4.3530    1.9130

plot(logxm1,logy,'b+', logxm2,logy,'r+', logxm4,logy,'g+', ...
     log(paramHat1(1)) + logy/paramHat1(2),logy,'b--', ...
     log(paramHat2(1)) + logy/paramHat2(2),logy,'r--', ...
     log(paramHat4(1)) + logy/paramHat4(2),logy,'g--');
xlabel('log(x - c)');
ylabel('log(-log(1 - p))');
legend({'Threshold = 1' 'Threshold = 2' 'Threshold = 4'}, 'location','northwest');

Отношение между логарифмической (x-4) и логарифмической (-журнал (1-p)) выглядит приблизительно линейно. Поскольку мы ожидаем увидеть приблизительно линейный график, если вычесть истинный пороговый параметр, это является доказательством того, что 4 может быть разумным значением для порога. С другой стороны, графики для 2 и 3 более систематически отличаются от линейных, что является доказательством того, что эти значения не соответствуют данным.

Этот аргумент может быть формализован. Для каждого предварительного значения порогового параметра соответствующая предварительная подгонка Вейбула может быть охарактеризована как значения параметров, которые максимизируют значение R ^ 2 линейной регрессии на преобразованных переменных log (x-c) и журнал (-журнал (1-p)). Чтобы оценить пороговый параметр, мы можем перенести этот один шаг дальше и максимизировать значение R ^ 2 над всеми возможными пороговыми значениями.

r2 = @(x,y) 1 - norm(y - polyval(polyfit(x,y,1),x)).^2 / norm(y - mean(y)).^2;
threshObj = @(c) -r2(log(-log(1-p)),log(x-c));
cHat = fminbnd(threshObj,.75*min(x), .9999*min(x));
poly = polyfit(log(-log(1-p)),log(x-cHat),1);
paramHat = [exp(poly(2)) 1/poly(1) cHat]
logx = log(x-cHat);
logy = log(-log(1-p));
plot(logx,logy,'b+', log(paramHat(1)) + logy/paramHat(2),logy,'r--');
xlabel('log(x - cHat)');
ylabel('log(-log(1 - p))');

paramHat =

    4.7448    2.3839    3.6029

Шкалы

Экспоненциальное распределение является семейством шкал, и в журнал масштабе распределение Вейбула является семейством шкал местоположения, поэтому этот метод наименьших квадратов был простым в этих двух случаях. Общая процедура для соответствия распределению шкалы местоположения

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

Но другие распределения, которые не являются семействами шкалы местоположения, как гамма, немного хитрее. Нет преобразования CDF, которое даст линейную связь. Однако мы можем использовать подобную идею, только на этот раз работая над нетрансформированной совокупной шкалой вероятностей. График P-P является подходящим способом визуализации этой процедуры аппроксимации.

Если эмпирические вероятности из ECDF построены относительно подобранных вероятностей из параметрической модели, плотное рассеяние вдоль линии 1:1 от нуля до единицы указывает, что значения параметров определяют распределение, которое хорошо объясняет наблюдаемые данные, потому что подобранная CDF аппроксимирует эмпирическую скважину CDF. Идея состоит в том, чтобы найти значения параметров, которые делают график вероятности как можно ближе к линии 1:1. Это может даже не быть возможным, если распределение не является хорошей моделью для данных. Если график P-P показывает систематическое отклонение от линии 1:1, то модель может оказаться сомнительной. Однако важно помнить, что поскольку точки на этих графиках не являются независимыми, интерпретация не совсем такая же, как регрессионный остаточный график.

Для примера мы моделируем некоторые данные и подбираем гамма- распределение.

n = 100;
x = gamrnd(2,1,n,1);

Вычислите ECDF х.

x = sort(x);
pEmp = ((1:n)-0.5)' ./ n;

Мы можем сделать график вероятности, используя любое начальное предположение для параметров гамма-распределения, например, a = 1 и b = 1. Это предположение не очень хорошо - вероятности из параметрического CDF не близки к вероятностям из ECDF. Если бы мы попробовали другие a и b, мы бы получили другое рассеяние на графике P-P, с другим расхождением с линией 1:1. Поскольку мы знаем истинные a и b в этом примере, мы попробуем эти значения.

a0 = 1; b0 = 1;
p0Fit = gamcdf(x,a0,b0);
a1 = 2; b1 = 1;
p1Fit = gamcdf(x,a1,b1);
plot([0 1],[0 1],'k--', pEmp,p0Fit,'b+', pEmp,p1Fit,'r+');
xlabel('Empirical Probabilities');
ylabel('(Provisionally) Fitted Gamma Probabilities');
legend({'1:1 Line','a=1, b=1', 'a=2, b=1'}, 'location','southeast');

Второе множество значений для a и b делает для намного лучшего графика, и, таким образом, более совместимы с данными, если вы измеряете «совместимость» по тому, как прямо вы можете сделать график P-P.

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

wgt = 1 ./ sqrt(pEmp.*(1-pEmp));
gammaObj = @(params) sum(wgt.*(gamcdf(x,exp(params(1)),exp(params(2)))-pEmp).^2);
paramHat = fminsearch(gammaObj,[log(a1),log(b1)]);
paramHat = exp(paramHat)

paramHat =

    2.2759    0.9059

pFit = gamcdf(x,paramHat(1),paramHat(2));
plot([0 1],[0 1],'k--', pEmp,pFit,'b+');
xlabel('Empirical Probabilities');
ylabel('Fitted Gamma Probabilities');

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

Модель Мисспецификации

График P-P также может быть полезен для сравнения подгонек из различных семейств распределения. Что произойдет, если мы попытаемся подогнать логнормальное распределение к этим данным?

wgt = 1 ./ sqrt(pEmp.*(1-pEmp));
LNobj = @(params) sum(wgt.*(logncdf(x,params(1),exp(params(2)))-pEmp).^2);
mu0 = mean(log(x)); sigma0 = std(log(x));
paramHatLN = fminsearch(LNobj,[mu0,log(sigma0)]);
paramHatLN(2) = exp(paramHatLN(2))

paramHatLN =

    0.5331    0.7038

pFitLN = logncdf(x,paramHatLN(1),paramHatLN(2));
hold on
plot(pEmp,pFitLN,'rx');
hold off
ylabel('Fitted Probabilities');
legend({'1:1 Line', 'Fitted Gamma', 'Fitted Lognormal'},'location','southeast');

Заметьте, как логнормальная подгонка систематически отличается от гамма-подгонки в хвостах. Медленнее растёт в левом хвосте, медленнее гибнет в правом хвосте. Гамма, по-видимому, немного лучше соответствует данным.

Пример Lognormal Threshold Parameter

Lognormal distribution является простым в соответствии с максимальной вероятностью, потому что, как только преобразование логарифмического сигнала применяется к данным, максимальная вероятность идентична подбору кривой нормального. Но иногда необходимо оценить параметр порога в логнормальной модели. Вероятность для такой модели неограниченна, и поэтому максимальная вероятность не работает. Однако метод наименьших квадратов предоставляет способ сделать оценки. Поскольку двухпараметровое логнормальное распределение может быть логарифмически преобразовано в семейство шкал местоположения, мы могли бы следовать тем же шагам, что и в более раннем примере, который показал подбор кривой распределения Вейбула с пороговым параметром. Здесь, однако, мы сделаем оценку по совокупной шкале вероятностей, как в предыдущем примере, показывающем подгонку с гамма- распределением.

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

n = 200;
x = lognrnd(0,.5,n,1) + 10;
hist(x,20); xlim([8 15]);

Вычислите ECDF x и найдите параметры для оптимального трехпараметрического логнормального распределения.

x = sort(x);
pEmp = ((1:n)-0.5)' ./ n;
wgt = 1 ./ sqrt(pEmp.*(1-pEmp));
LN3obj = @(params) sum(wgt.*(logncdf(x-params(3),params(1),exp(params(2)))-pEmp).^2);
c0 = .99*min(x);
mu0 = mean(log(x-c0)); sigma0 = std(log(x-c0));
paramHat = fminsearch(LN3obj,[mu0,log(sigma0),c0]);
paramHat(2) = exp(paramHat(2))

paramHat =

   -0.0698    0.5930   10.1045

pFit = logncdf(x-paramHat(3),paramHat(1),paramHat(2));
plot(pEmp,pFit,'b+', [0 1],[0 1],'k--');
xlabel('Empirical Probabilities');
ylabel('Fitted 3-param Lognormal Probabilities');

Меры точности

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

Однако симуляция Монте-Карло предоставляет другой способ оценки стандартных ошибок. Если мы используем подобранную модель, чтобы сгенерировать большое количество наборов данных, мы можем аппроксимировать стандартную ошибку оценщиков со стандартным отклонением Монте-Карло. Для простоты мы определили функцию аппроксимации в отдельном файле, logn3fit.m.

estsSim = zeros(1000,3);
for i = 1:size(estsSim,1)
    xSim = lognrnd(paramHat(1),paramHat(2),n,1) + paramHat(3);
    estsSim(i,:) = logn3fit(xSim);
end
std(estsSim)

ans =

    0.1542    0.0908    0.1303

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

subplot(3,1,1), hist(estsSim(:,1),20);
title('Log-Location Parameter Bootstrap Estimates');
subplot(3,1,2), hist(estsSim(:,2),20);
title('Log-Scale Parameter Bootstrap Estimates');
subplot(3,1,3), hist(estsSim(:,3),20);
title('Threshold Parameter Bootstrap Estimates');

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

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

[paramHat; mean(estsSim)]

ans =

   -0.0698    0.5930   10.1045
   -0.0690    0.5926   10.0905

Наконец, мы могли бы также использовать функцию bootstrp для вычисления стандартных оценок ошибок bootstrap. Они не делают никаких параметрических предположений о данных.

estsBoot = bootstrp(1000,@logn3fit,x);
std(estsBoot)

ans =

    0.1490    0.0785    0.1180

Стандартные ошибки bootstrap не за горами от расчетов Монте-Карло. Это неудивительно, поскольку подобранная модель является такой же, из которой были сгенерированы данные примера.

Сводные данные

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

Графики P-P, используемые здесь для иллюстрации метода аппроксимации, полезны сами по себе, как визуальное указание на отсутствие подгонки при подборе одномерного распределения.