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

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

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

Поскольку значения для данных, подвергнутых цензуре, точно не известны, оценка максимального правдоподобия становится более сложной задачей. В частности, как PDF, так и CDF необходимы для вычисления логарифмической вероятности. Поэтому необходимо предоставить mle с функциями для обоих, чтобы соответствовать цензурированным данным. Инструментарий статистики и машинного обучения включает в себя следующие функции: 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 и sigma, которые соответствуют наблюдаемому среднему и стандартному отклонению неподцензурных данных.

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. Поскольку параметр масштаба, sigma, должен быть положительным, мы указываем более низкие границы параметров. mle возвращает оценки максимального правдоподобия двух параметров распределения экстремальных значений, mu и sigma, а также приблизительные 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» должен быть «выключен», что отключит проверку на наличие неграничных значений правдоподобия, а затем надеется на лучшее. Но правильный способ решения этой численной задачи лежит в ее основе, и в данном случае это не трудно сделать.

Обратите внимание, что крайнее значение CDF имеет вид

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

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

SFval = 1e-20;
CDFval = 1 - SFval

CDFval = 1

Однако журнал регистрации соответствующего значения SF, т.е. log (1-CDF), легко представляется.

log(SFval)

ans = -46.0517

Аналогичное наблюдение можно сделать по поводу использования log PDF, а не самого PDF - вклад незацензурных наблюдений в log-правдоподобие является журналом их значений 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». В отличие от функций 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 ® включает Toolbox™ оптимизации ,mle позволяет использовать функцию fmincon, который включает в себя алгоритмы оптимизации, которые могут использовать производную информацию. Чтобы использовать наилучшие преимущества алгоритмов в fmincon, можно указать пользовательское распределение с помощью функции log-правдоподобия, записанной для возврата не только самого log-правдоподобия, но и его градиента. Градиент логарифмической функции правдоподобия является просто вектором его частных производных относительно его параметров.

Эта стратегия требует дополнительной подготовки для написания кода, который вычисляет как логарифмическое правдоподобие, так и его градиент. Для этого примера мы создали код, чтобы сделать это для распределения экстремальных значений в виде отдельного файла 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