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

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

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 и 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, мы также используем параметр 'censoring', чтобы пройти в векторе цензуры, с. Потому что параметр шкалы, 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' управляет параметром 'off', который отключал бы проверку на несрочные значения правдоподобия, а затем надеялся на лучшее. Но правильный способ решить эту числовую задачу в корне, и в этом случае сделать это не сложно.

Заметьте, что экстремальное значение CDF имеет вид

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

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