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

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

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

rng(2,'twister');
lifetime = [wblrnd(15000,3,90,1); wblrnd(1500,3,10,1)];

В этом примере предположим, что мы тестируем дроссели в стрессовых условиях, так что каждый час тестирования эквивалентен 100 часам фактического использования в полевых условиях. По прагматическим причинам часто бывает, что тесты надежности останавливаются через определенное время. Для этого примера мы будем использовать 140 часов, эквивалентных в общей сложности 14 000 часам реального обслуживания. Некоторые элементы отказывают во время теста, в то время как другие выживают в течение 140 часов. В реальном тесте время для последнего было бы записано как 14000, и мы имитируем это в смоделированных данных. Также обычной практикой является сортировка времени отказа.

T = 14000;
obstime = sort(min(T, lifetime));

Мы знаем, что любые дроссели, которые переживут испытание, в конечном итоге потерпят неудачу, но испытание не достаточно долго, чтобы наблюдать их фактическое время до отказа. Известно, что их продолжительность жизни превышает 14 000 часов. Говорят, что эти ценности подвергаются цензуре. Этот сюжет показывает, что около 40% наших данных подвергаются цензуре на уровне 14 000.

failed = obstime(obstime<T); nfailed = length(failed);
survived = obstime(obstime==T); nsurvived = length(survived);
censored = (obstime >= T);
plot([zeros(size(obstime)),obstime]', repmat(1:length(obstime),2,1), ...
     'Color','b','LineStyle','-')
line([T;3e4], repmat(nfailed+(1:nsurvived), 2, 1), 'Color','b','LineStyle',':');
line([T;T], [0;nfailed+nsurvived],'Color','k','LineStyle','-')
text(T,30,'<--Unknown survival time past here')
xlabel('Survival time'); ylabel('Observation number')

Способы просмотра распределений

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

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

x = linspace(1,30000);
subplot(2,2,1);
plot(x,wblpdf(x,14000,2),x,wblpdf(x,18000,2),x,wblpdf(x,14000,1.1))
title('Prob. Density Fcn')
subplot(2,2,2);
plot(x,1-wblcdf(x,14000,2),x,1-wblcdf(x,18000,2),x,1-wblcdf(x,14000,1.1))
title('Survivor Fcn')
subplot(2,2,3);
wblhaz = @(x,a,b) (wblpdf(x,a,b) ./ (1-wblcdf(x,a,b)));
plot(x,wblhaz(x,14000,2),x,wblhaz(x,18000,2),x,wblhaz(x,14000,1.1))
title('Hazard Rate Fcn')
subplot(2,2,4);
probplot('weibull',wblrnd(14000,2,40,1))
title('Probability Plot')

Подгонка распределения Вейбулла

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

Другие распределения, используемые для моделирования данных о времени жизни, включают в себя распределения логнормального, гамма и Бирнбаума-Сондерса.

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

subplot(1,1,1);
[empF,x,empFlo,empFup] = ecdf(obstime,'censoring',censored);
stairs(x,empF);
hold on;
stairs(x,empFlo,':'); stairs(x,empFup,':');
hold off
xlabel('Time'); ylabel('Proportion failed'); title('Empirical CDF')

Этот график показывает, например, что доля отказа к времени 4000 составляет около 12%, и 95% уверенности, связанной с вероятностью отказа к этому времени, составляет от 6% до 18%. Обратите внимание, что, поскольку наш тест длился всего 14 000 часов, эмпирический CDF позволяет нам вычислять вероятности отказов только до этого предела. Почти половина данных была подвергнута цензуре на уровне 14 000, и поэтому эмпирический CDF повышается только до около 0,53, вместо 1,0.

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

paramEsts = wblfit(obstime,'censoring',censored);
[nlogl,paramCov] = wbllike(paramEsts,obstime,censored);
xx = linspace(1,2*T,500);
[wblF,wblFlo,wblFup] = wblcdf(xx,paramEsts(1),paramEsts(2),paramCov);

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

stairs(x,empF);
hold on
handles = plot(xx,wblF,'r-',xx,wblFlo,'r:',xx,wblFup,'r:');
hold off
xlabel('Time'); ylabel('Fitted failure probability'); title('Weibull Model vs. Empirical')

Обратите внимание, что модель Вейбулла позволяет проецировать и вычислять вероятности отказов на время после окончания теста. Тем не менее, похоже, что подогнанная кривая плохо соответствует нашим данным. У нас слишком много ранних сбоев до времени 2000 по сравнению с тем, что предсказывала бы модель Вейбулла, и, как следствие, слишком мало времени между примерно 7000 и примерно 13000. Это неудивительно - вспомним, что мы создавали данные именно таким образом.

Добавление гладкой непараметрической оценки

Предварительно определенные функции, предоставляемые Toolbox™ статистики и машинного обучения, не включают в себя распределения, которые имеют избыток ранних сбоев, подобных этому. Вместо этого может потребоваться нарисовать гладкую непараметрическую кривую через эмпирический CDF, используя функцию ksdensity. Мы удалим доверительные полосы для CDF Вейбулла и добавим две кривые, одна с параметром сглаживания по умолчанию, а другая с параметром сглаживания 1/3 значение по умолчанию. Меньший параметр сглаживания позволяет кривой более точно следовать данным.

delete(handles(2:end))
[npF,ignore,u] = ksdensity(obstime,xx,'cens',censored,'function','cdf');
line(xx,npF,'Color','g');
npF3 = ksdensity(obstime,xx,'cens',censored,'function','cdf','width',u/3);
line(xx,npF3,'Color','m');
xlim([0 1.3*T])
title('Weibull and Nonparametric Models vs. Empirical')
legend('Empirical','Fitted Weibull','Nonparametric, default','Nonparametric, 1/3 default', ...
       'location','northwest');

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

Давайте рассчитаем степень опасности для этой непараметрической посадки и построим ее график по диапазону данных.

hazrate = ksdensity(obstime,xx,'cens',censored,'width',u/3) ./ (1-npF3);
plot(xx,hazrate)
title('Hazard Rate for Nonparametric Model')
xlim([0 T])

Эта кривая имеет немного форму «ванны», с высоким уровнем опасности около 2000, падает до более низких значений, затем снова повышается. Это типично для уровня опасности для компонента, который более подвержен отказу в раннем возрасте (младенческая смертность) и снова в более позднем возрасте (старение).

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

Альтернативные модели

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

Одним из вариантов было бы использование другого параметрического распределения. Инструментарий статистики и машинного обучения включает функции для других распространенных распределений времени жизни, таких как логнормальное, гамма и Бирнбаум-Сондерс, а также для многих других распределений, которые обычно не используются в моделях времени жизни. Можно также определить и подогнать пользовательские параметрические модели к данным о сроке службы, как описано в примере Фитинг пользовательских одномерных распределений, деталь 2.

Другой альтернативой было бы использование смеси двух параметрических распределений: одно представляет ранний сбой, а другое представляет остальную часть распределения. Фитинг-смеси распределений описаны в примере Фитинг пользовательских одномерных распределений (Fitting Custom Univariate Distributions).