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

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

Statistics and Machine Learning Toolbox включает функциональный evfit, который соответствует распределению экстремума к данным, включая данные с цензурированием. Однако в целях этого примера, мы проигнорируем evfit, и покажите, как использовать mle и пользовательские дистрибутивы, чтобы подобрать модель к подвергнутым цензуре данным, с помощью распределения экстремума.

Поскольку значения для подвергнутых цензуре данных не известны точно, оценка наибольшего правдоподобия становится более трудной. В частности, и PDF и CDF необходимы, чтобы вычислить логарифмическую правдоподобность. Поэтому необходимо обеспечить mle с функциями для обоих из тех для того, чтобы соответствовать подвергнутым цензуре данным. Statistics and Machine Learning Toolbox включает функции evpdf и evcdf, таким образом для этого примера, работа записи кода была уже сделана.

Мы будем подбирать модель к симулированным данным. Первый шаг должен сгенерировать некоторые не прошедшие цензуру данные об экстремуме.

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

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

c = (x > 7);
x(c) = 7;
sum(c)/length(c)

ans = 0.1200

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

[uncensCnts,binCtrs] = hist(x(~c));
censCnts = hist(x(c),binCtrs);
bar(binCtrs,[uncensCnts' censCnts'],'stacked');

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

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

sigma0 = 2.3495

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

mu0 = 3.5629

В дополнение к передаче данных, x, и указателей на PDF и CDF функционирует в mle, мы также используем параметр 'цензурирования', чтобы передать в векторе цензурирования, c. Поскольку масштабный коэффициент, сигма, должен быть положительным, мы задаем более низкие границы параметра. mle возвращает оценки наибольшего правдоподобия двух параметров распределения экстремума, mu и сигмы, а также аппроксимированные 95% доверительных интервалов.

[paramEsts,paramCIs] = mle(x, 'censoring',c, 'pdf',@evpdf, 'cdf',@evcdf, ...
                           'start',[mu0 sigma0], 'lower',[-Inf,0])

paramEsts = 1×2

    4.5530    3.0215

paramCIs = 2×2

    3.6455    2.2937
    5.4605    3.7494

Некоторые числовые проблемы в подборе кривой пользовательским дистрибутивам

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

Во-первых, одно из значений PDF может стать столь маленьким, что оно потеряло значимость, чтобы обнулить в арифметике двойной точности. Во-вторых, одно из значений CDF может стать так близко к 1, что оно окружило в двойной точности. (Также возможно, что значение CDF может стать столь маленьким, что теряет значимость, но это, оказывается, не проблема.)

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

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

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

The CDF function returned values greater than or equal to 1.

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

Заметьте, что экстремум CDF имеет форму

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

Вклад подвергнутых цензуре наблюдений к логарифмической правдоподобности является журналом их значений функции выживания (SF), т.е. журналом (1-CDF). Для распределения экстремума журнал SF только-exp ((x-mu)./сигма). Если мы могли бы вычислить логарифмическую правдоподобность с помощью логарифмического SF непосредственно, (вместо, в действительности, вычислив журнал (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, т.е. журнал (1-CDF), легко представлен.

log(SFval)

ans = -46.0517

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

logPDFval = -921;
PDFval = exp(logPDFval)

PDFval = 0

mle обеспечивает синтаксис для определения пользовательского дистрибутива с помощью журнала PDF и логарифмический SF (а не PDF и CDF), через 'logpdf' и 'logsf' параметры. В отличие от PDF и функций CDF, нет никаких существующих функций, таким образом, мы создадим анонимные функции, которые вычисляют эти значения:

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

Используя ту же начальную точку, альтернатива logPDF/logSF спецификация распределения экстремума делает проблему разрешимой:

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

paramEsts = 1×2

    4.5530    3.0215

paramCIs = 2×2

    3.6455    2.2937
    5.4605    3.7494

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

Предоставление градиента

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

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

Если ваша установка MATLAB® включает Optimization Toolbox™, mle позволяет вам использовать функциональный fmincon, который включает алгоритмы оптимизации, которые могут использовать производную информацию. Воспользоваться лучшим преимуществом алгоритмов в fmincon, можно задать пользовательский дистрибутив с помощью функции логарифмической правдоподобности, записанной, чтобы возвратить не только саму логарифмическую правдоподобность, но и ее градиент также. Градиент функции логарифмической правдоподобности является просто вектором своих частных производных относительно ее параметров.

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

type evnegloglike.m

function [nll,ngrad] = evnegloglike(params,x,cens,freq)
%EVNEGLOGLIKE Negative log-likelihood for the extreme value distribution.

%   Copyright 1984-2004 The MathWorks, Inc.

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

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

Чтобы вычислить оценки наибольшего правдоподобия с помощью основанного на градиенте алгоритма оптимизации, мы используем 'nloglf' параметр, указывая, что мы обеспечиваем указатель на функцию, которая вычисляет отрицательную логарифмическую правдоподобность и 'optimfun' параметр, задавая fmincon как оптимизационная функция. mle автоматически обнаружит тот evnegloglike может возвратить и отрицательную логарифмическую правдоподобность и ее градиент.

start = [1 1];
[paramEsts,paramCIs] = mle(x, 'censoring',c, 'nloglf',@evnegloglike, ...
                           'start',start, 'lower',[-Inf,0], 'optimfun','fmincon')

paramEsts = 1×2

    4.5530    3.0215

paramCIs = 2×2

    3.6455    2.2937
    5.4605    3.7493