Этот пример показывает, как соответствовать одномерным распределениям с помощью оценок методом наименьших квадратов кумулятивных функций распределения. Это - обычно применимый метод, который может быть полезным в случаях, когда наибольшее правдоподобие перестало работать, например, некоторые модели, которые включают пороговый параметр.
Наиболее распространенный метод для подбора кривой одномерному распределению к данным является наибольшим правдоподобием. Но наибольшее правдоподобие не работает во всех случаях, и другие методы оценки, такие как Метод Моментов, иногда необходимы. Когда применимо, наибольшее правдоподобие является, вероятно, лучшим выбором методов, потому что это часто более эффективно. Но метод, описанный здесь, обеспечивает другой инструмент, который может использоваться при необходимости.
Термин "наименьшие квадраты" обычно используется в контексте подбора кривой линии регрессии или поверхности, чтобы смоделировать переменную отклика как функцию одной или нескольких переменных прогноза. Метод, описанный здесь, является совсем другим приложением наименьших квадратов: подбор кривой одномерного распределения, только с одной переменной.
Чтобы начаться, сначала моделируйте некоторые выборочные данные. Мы будем использовать экспоненциальное распределение, чтобы сгенерировать данные. В целях этого примера, как на практике, мы примем, что данные, как известно, не прибыли из конкретной модели.
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 влияли больше наблюдения в хвосте распределения.
Для немного более комплексного примера моделируйте некоторые выборочные данные от распределения Weibull и вычислите ECDF x.
n = 100; x = wblrnd(2,1,n,1); x = sort(x); p = ((1:n)-0.5)' ./ n;
Чтобы соответствовать распределению Weibull к этим данным, заметьте, что CDF для Weibull является p = PR {X <= x} = 1 - exp (-(x/a) ^b). При преобразовании этого, чтобы регистрировать (a) + журнал (-журнал (1-p)) * (1/b) = журнал (x) снова дает линейное соотношение, на этот раз между журналом (-журнал (1-p)) и журналом (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 или логарифмически нормальный с пороговым параметром. Например, случайная переменная Weibull принимает значения по (0, Inf), и пороговый параметр, c, сдвиги, которые располагаются к (c, Inf). Если пороговый параметр известен, то нет никакой трудности. Но если пороговый параметр не известен, он должен вместо этого быть оценен. Этим моделям трудно соответствовать наибольшему правдоподобию - вероятность может иметь несколько режимов, или даже стать бесконечной для значений параметров, которые не разумны для данных, и таким образом, наибольшее правдоподобие часто является не хорошим методом. Но с маленьким добавлением к процедуре наименьших квадратов, мы можем получить устойчивые оценки.
Чтобы проиллюстрировать, мы моделируем некоторые данные из распределения Weibull с тремя параметрами с пороговым значением. Как выше, мы примем в целях примера, что данные, как известно, не прибыли из конкретной модели, и что порог не известен.
n = 100; x = wblrnd(4,2,n,1) + 4; hist(x,20); xlim([0 16]);
Как мы можем соответствовать распределению Weibull с тремя параметрами к этим данным? Если бы мы знали, каково пороговое значение было, 1, например, мы могли вычесть то значение из данных и затем использовать процедуру наименьших квадратов, чтобы оценить форму Weibull и масштабные коэффициенты.
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) и журнал (-журнал (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 отличаются более систематически от линейного, который является доказательством, что те значения не сопоставимы с данными.
Этот аргумент может быть формализован. Для каждого предварительного значения порогового параметра соответствующая временная подгонка Weibull может быть охарактеризована как значения параметров, которые максимизируют значение R^2 линейной регрессии на преобразованном журнале переменных (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
Экспоненциальное распределение является семейством шкал, и на логарифмической шкале, распределение Weibull является семейством шкал местоположения, таким образом, этот метод наименьших квадратов был прямым в тех двух случаях. Общая процедура, чтобы соответствовать распределению шкалы местоположения
Вычислите ECDF наблюдаемых данных.
Преобразуйте CDF распределения, чтобы получить линейное соотношение между некоторой функцией данных и некоторой функцией интегральной вероятности. Эти две функции не включают параметры распределения, но наклон и прерывание строки делают.
Включите значения x, и p от ECDF в это преобразовал CDF, и соответствуйте прямой линии с помощью наименьших квадратов.
Решите для параметров распределения с точки зрения наклона и прерывания строки.
Мы также видели, что подбор кривой распределению, которое является семейством шкал местоположения с дополнительным, пороговый параметр является незначительно более трудным.
Но другие дистрибутивы, которые не являются семействами шкал местоположения, как гамма, немного более хитры. Нет никакого преобразования 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;
Мы можем сделать график вероятности с помощью любого исходного предположения для гамма параметров распределения, 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');
Заметьте, как логарифмически нормальная подгонка систематически отличается от гаммы, помещаются в хвосты. Это растет более медленно в левом хвосте и умирает более медленно в правом хвосте. Гамма, кажется, немного лучшая подгонка к данным.
Логарифмически нормальное распределение просто соответствовать наибольшим правдоподобием, потому что, если логарифмическое преобразование применяется к данным, наибольшее правдоподобие идентично подбору кривой нормальному. Но иногда необходимо оценить пороговый параметр в логарифмически нормальной модели. Вероятность для такой модели неограниченна, и таким образом, наибольшее правдоподобие не работает. Однако метод наименьших квадратов обеспечивает способ сделать оценки. Начиная с 2D параметра логарифмически нормальное распределение может быть преобразовано в журнал семейству шкал местоположения, мы могли выполнить те же шаги как в более раннем примере, который показал подбор кривой распределению Weibull с пороговым параметром. Здесь, однако, мы сделаем оценку в шкале интегральной вероятности, как в предыдущем примере, показывающем подгонку с гамма распределением.
Чтобы проиллюстрировать, мы моделируем некоторые данные из логарифмически нормального распределения с тремя параметрами с порогом.
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, чтобы вычислить оценки стандартной погрешности начальной загрузки. Они не делают параметрических предположений о данных.
estsBoot = bootstrp(1000,@logn3fit,x); std(estsBoot)
ans =
0.1490 0.0785 0.1180
Стандартные погрешности начальной загрузки не далеки от вычислений Монте-Карло. Это не удивительно, поскольку подобранная модель - тот же, от которого были сгенерированы данные в качестве примера.
Подходящий метод, описанный здесь, является альтернативой наибольшему правдоподобию, которое может использоваться, чтобы соответствовать одномерным распределениям, когда наибольшему правдоподобию не удается обеспечить полезные оценки параметра. Одно важное приложение находится в подходящих дистрибутивах, включающих пороговый параметр, такой как логарифмически нормальное с тремя параметрами. Стандартные погрешности более трудно вычислить, чем для оценок наибольшего правдоподобия, потому что аналитические приближения не существуют, но симуляция обеспечивает выполнимую альтернативу.
Графики P-P, используемые здесь, чтобы проиллюстрировать подходящий метод, полезны самостоятельно как визуальная индикация относительно отсутствия подгонки при подборе кривой одномерному распределению.