Задайте соответствующие начальные значения параметров

Поскольку значения для подвергнутых цензуре данных не известны точно, оценка наибольшего правдоподобия запрашивает больше информации. В частности, функция плотности вероятности (PDF), кумулятивная функция распределения (cdf) и соответствующие начальные значения параметров необходима, чтобы вычислить логарифмическую правдоподобность. Можно использовать evpdf и evcdf функции, чтобы задать PDF и cdf.

Во-первых, сгенерируйте некоторые не прошедшие цензуру данные об экстремуме.

rng(0,'twister');
n = 50;
mu = 5;
sigma = 2.5;
x = evrnd(mu,sigma,n,1);

Затем подвергают цензуре любые значения, которые больше, чем предопределенное сокращение, заменяя те значения на значение сокращения. Этот тип цензурирования известен как цензурирование Типа II.

c = (x > 7);
x(c) = 7;

Проверяйте процент подвергнутых цензуре наблюдений.

sum(c)/length(c)

ans = 0.1200

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

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

[uncensCnts,Edges] = histcounts(x(~c),10);
censCnts = histcounts(x(c),Edges);
bar(Edges(1:end-1)+diff(Edges)/2,[uncensCnts' censCnts'],'stacked')
legend('Fully observed data','Censored data','Location','northwest')

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

sigma0 = sqrt(6)*std(x(~c))./pi

sigma0 = 2.3495

mu0 = mean(x(~c))-psi(1).*sigma0

mu0 = 3.5629

Найдите оценки наибольшего правдоподобия двух параметров распределения экстремума, а также аппроксимированные 95% доверительных интервалов. Задайте вектор цензурирования, PDF, cdf, и начальные значения параметров. Поскольку sigma (масштабный коэффициент) должен быть положительным, также необходимо задать более низкие границы параметра.

[paramEsts,paramCIs] = mle(x,'Censoring',c, ...
    'pdf',@evpdf,'cdf',@evcdf, ...
    'Start',[mu0 sigma0],'LowerBound',[-Inf,0])

paramEsts = 1×2

    4.5530    3.0215

paramCIs = 2×2

    3.6455    2.2937
    5.4605    3.7494

Задайте logpdf и logsf

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

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

Любое условие может вызвать проблемы когда mle вычисляет логарифмическую правдоподобность, потому что каждый приводит к значениям логарифмической правдоподобности —Inf, который алгоритм оптимизации в mle не может обработать.

Исследуйте то, что происходит с различной начальной точкой.

start = [1 1];
try
    [paramEsts,paramCIs] = mle(x,'Censoring',c, ...
        'pdf',@evpdf,'cdf',@evcdf, ...
        'Start',start,'LowerBound',[-Inf,0])
catch ME
    disp(ME.message)
end

Custom cumulative distribution function returned values greater than or equal to 1.

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

Экстремум cdf имеет форму

  p = 1 - exp(-exp((x-mu)./sigma))

Вклад подвергнутых цензуре наблюдений к логарифмической правдоподобности является журналом их значений функции выживания (SF) или log(1-cdf). Для распределения экстремума журналом SF является -exp((x-mu)./sigma). Если вы вычисляете логарифмическую правдоподобность с помощью логарифмического SF непосредственно, вместо того, чтобы вычислить log(1-(1-exp(logSF))), можно избежать округляющихся проблем с cdf. Наблюдения, cdf значения которых не различимы от 1 в двойной точности имеют логарифмические значения SF, которые являются легко представимыми как ненулевые значения. Например, cdf значение (1-1e-20) раунды к 1 в двойной точности, потому что двойная точность eps о 2e-16.

SFval = 1e-20;
cdfval = 1 - SFval

cdfval = 1

Программное обеспечение может легко представлять журнал соответствующего значения SF.

log(SFval)

ans = -46.0517

Та же ситуация верна для журнала PDF; вклад не прошедших цензуру наблюдений к логарифмической правдоподобности является журналом их значений PDF. Можно использовать журнал PDF непосредственно, вместо того, чтобы вычислить log(exp(logpdf)), избегать проблем потери значимости, где PDF не различима от нуля в двойной точности. Программное обеспечение может легко представлять журнал PDF как конечное отрицательное число. Например, значение PDF 1e-400 потери значимости в двойной точности, потому что двойная точность realmin о 2e-308.

logpdfval = -921;
pdfval = exp(logpdfval)

pdfval = 0

Используя mle функция, можно задать пользовательский дистрибутив с журналом PDF и логарифмический SF (а не PDF и cdf) путем установки logpdf и logsf аргументы name-value. В отличие от PDF и функций cdf, регистрируйте PDF и регистрируйте SF, не имеют встроенных функций. Поэтому необходимо создать анонимные функции, которые вычисляют эти значения.

evlogpdf = @(x,mu,sigma) ((x-mu)./sigma - exp((x-mu)./sigma)) - log(sigma);
evlogsf = @(x,mu,sigma) -exp((x-mu)./sigma);

Используя ту же начальную точку, альтернативный журнал спецификация SF PDF и журнала распределения экстремума делает проблему разрешимой.

start = [1 1];
[paramEsts,paramCIs] = mle(x,'Censoring',c, ...
    'logpdf',evlogpdf,'logsf',evlogsf, ...
    'Start',start,'LowerBound',[-Inf,0])

paramEsts = 1×2

    4.5530    3.0215

paramCIs = 2×2

    3.6455    2.2937
    5.4605    3.7494

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

Предоставьте градиент к оптимизационной функции fmincon

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

Можно задать OptimFun аргумент значения имени в mle как fmincon использовать fmincon функция (требует Optimization Toolbox). fmincon функция включает алгоритмы оптимизации, которые могут использовать производную информацию. Использовать в своих интересах алгоритмы в fmincon, задайте пользовательский дистрибутив с помощью функции логарифмической правдоподобности, записанной, чтобы возвратить не только логарифмическую правдоподобность, но и ее градиент также. Градиент функции логарифмической правдоподобности является вектором из своих частных производных относительно ее параметров.

Эта стратегия требует, чтобы дополнительная подготовка записала код, который вычисляет и логарифмическую правдоподобность и ее градиент. Задайте функцию с именем helper_evnegloglike в отдельном файле.

function [nll,ngrad] = helper_evnegloglike(params,x,cens,freq)
%HELPER_EVNEGLOGLIKE Negative log-likelihood for the extreme value
% distribution.
% This function supports only the example Avoid Numerical Issues When
% Fitting Custom Distributions (customdist2demo.mlx) and might change in 
% a future release.


if numel(params)~=2
    error(message('stats:probdists:WrongParameterLength',2));
end
mu = params(1);
sigma = params(2);
nunc = sum(1-cens);
z = (x - mu) ./ sigma;
expz = exp(z);
nll = sum(expz) - sum(z(~cens)) + nunc.*log(sigma);
if nargout > 1
    ngrad = [-sum(expz)./sigma + nunc./sigma, ...
        -sum(z.*expz)./sigma + sum(z(~cens))./sigma + nunc./sigma];
end

Функциональный helper_evnegloglike возвращает отрицание и значений логарифмической правдоподобности и значений градиента потому что mle минимизирует отрицательную логарифмическую правдоподобность.

Чтобы вычислить оценки наибольшего правдоподобия с помощью основанного на градиенте алгоритма оптимизации, задайте nloglf, OptimFun, и Options аргументы name-value. nloglf задает пользовательскую функцию, которая вычисляет отрицательную логарифмическую правдоподобность, OptimFun задает fmincon как оптимизационная функция и Options задает тот fmincon использует второй выход пользовательской функции для градиента.

start = [1 1];
[paramEsts,paramCIs] = mle(x,'Censoring',c,'nloglf',@helper_evnegloglike, ...
    'Start',start,'LowerBound',[-Inf,0], ...
    'OptimFun','fmincon','Options',statset('GradObj','on'))

paramEsts = 1×2

    4.5530    3.0215

paramCIs = 2×2

    3.6455    2.2937
    5.4605    3.7493