Анализ данных о выживании или надежности

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

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

Специальные свойства пожизненных данных

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

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

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

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

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')

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

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

  • Функция плотности вероятности (PDF) указывает на относительную вероятность отказа в разное время.

  • Функция оставшегося в живых дает вероятность выживания как функция времени и просто один минус кумулятивная функция распределения (1-CDF).

  • Показатель риска дает мгновенную вероятность отказа, данного выживание данному времени. Это - PDF, разделенный на функцию оставшегося в живых. В этом примере показатели риска оказываются увеличением, означая, что элементы более восприимчивы к отказу, когда время передает (старение).

  • График вероятности является перемасштабированным CDF и используется, чтобы сравнить данные с подходящим распределением.

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

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')

Подбор кривой распределению Weibull

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

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

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

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

Распределение Weibull часто является хорошей моделью для отказа оборудования. Функциональный wblfit соответствует распределению Weibull к данным, включая данные с цензурированием. После вычислительных оценок параметра мы оценим CDF для подбиравшей модели Weibull, с помощью тех оценок. Поскольку значения 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, чтобы судить как хорошо модели распределения Weibull данные о надежности дросселя.

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')

Заметьте, что модель Weibull позволяет нам проекту, и вычислите вероятности отказа в течение многих времен вне конца теста. Однако кажется, что кривая по экспериментальным точкам не совпадает с нашими данными хорошо. У нас есть слишком много ранних отказов перед временем 2,000 по сравнению с тем, что модель Weibull предсказала бы, и в результате лишь немногие в течение многих времен между приблизительно 7 000 и приблизительно 13 000. Это не удивительно - вспоминают, что мы сгенерировали данные с только этим видом поведения.

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

Предопределенные функции, которым предоставляют Statistics and Machine Learning Toolbox™, не включают распределений, которые имеют избыток ранних отказов как это. Вместо этого мы можем хотеть чертить сглаженную, непараметрическую кривую через эмпирический CDF, с помощью функционального ksdensity. Мы удалим полосы уверенности для Weibull 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])

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

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

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

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

Одна альтернатива должна была бы использовать различное параметрическое распределение. Statistics and Machine Learning Toolbox включает функции для других общих пожизненных распределений такой как логарифмически нормальное, гамма, и Бирнбаум-Сондерс, а также много других распределений, которые обычно не используются в пожизненных моделях. Можно также задать и подбирать пользовательские параметрические модели к пожизненным данным, как описано в Подходящих Пользовательских Одномерных распределениях, примере Части 2.

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