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

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

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

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). Преобразование этого в -log (1-p) * mu = x дает линейное отношение между -log (1-p) и x. Если данные происходят из экспоненциального, мы должны видеть, по крайней мере приблизительно, линейное отношение, если мы подключаем вычисленные значения x и p из ECDF в это уравнение. Если мы используем наименьшие квадраты, чтобы вписать прямую линию через начало координат в x против -log (1-p), то эта подгоняемая линия представляет экспоненциальное распределение, которое является «ближайшим» к данным. Наклон линии является оценкой параметра mu.

Эквивалентно, мы можем думать о y = -log (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 = -log (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) + log (-log (1-p)) * (1/b) = log (x) снова дает линейную связь, на этот раз между log (-log (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

Это не очень хорошая посадка - log (x-1) и log (-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');

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

Этот аргумент может быть формализован. Для каждого временного значения порогового параметра соответствующее временное вписывание Вейбулла можно охарактеризовать как значения параметра, которые максимизируют значение R ^ 2 линейной регрессии на преобразованных переменных log (x-c) и log (-log (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.

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

Мы можем сделать график вероятности, используя любое начальное предположение для параметров гамма-распределения, а = 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');

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

Пример параметра логнормального порога

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

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

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');

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

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

[paramHat; mean(estsSim)]

ans =

   -0.0698    0.5930   10.1045
   -0.0690    0.5926   10.0905

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

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

ans =

    0.1490    0.0785    0.1180

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

Резюме

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

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