Документация

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

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

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, и таким образом, логарифмическая вероятность бесконечна. Мы могли попытаться установить параметр управления 'FunValCheck' mle на '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