Этот пример показывает, как анализировать данные о жизни с помощью цензуры. В биологических или медицинских применениях это известно как анализ выживания, и время может представлять время выживания организма или время до излечения заболевания. В технических приложениях это известно как анализ надежности, и время может представлять время до отказа единицы оборудования.
Наш пример моделирует время до отказа дросселя от автомобильной системы впрыска топлива.
Некоторые функции пожизненных данных различают их другие типы данных. Во-первых, время жизни всегда является положительными значениями, обычно представляющими время. Во-вторых, некоторые сроки жизни могут не наблюдаться точно, так что известно, что они только больше некоторого значения. В-третьих, распределения и методы анализа, которые обычно используются, довольно специфичны для данных в течение жизни
Давайте симулируем результаты проверки 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% наших данных подвергаются цензуре в 14000.
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 и используется для сравнения данных с подобранным распределением.
Вот примеры этих четырех типов графиков, используя распределение Вейбула для иллюстрации. Вейбул является общим распределением для моделирования данных о жизни.
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')
Распределение Вейбула является обобщением экспоненциального распределения. Если время жизни следует экспоненциальному распределению, то они имеют постоянную скорость опасности. Это означает, что они не стареют, в том смысле, что вероятность наблюдения отказа в интервале, учитывая выживание до начала этого интервала, не зависит от того, с чего начинается интервал. Распределение Вейбула имеет скорость опасности, которая может увеличиться или уменьшиться.
Другие распределения, используемые для моделирования данных о времени жизни, включают распределения lognormal, gamma и Birnbaum-Saunders.
Мы построим эмпирическую совокупную функцию распределения наших данных, показав долю, выходящую из строя до каждого возможного времени выживания. Пунктирные кривые дают 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')
Заметьте, что модель Weibull позволяет нам проектировать и вычислять вероятности отказа на время после окончания теста. Однако, похоже, что подобранная кривая плохо соответствует нашим данным. У нас слишком много ранние отказы до времени 2000 по сравнению с тем, что предсказала бы модель Weibull, и, как следствие, слишком мало для времени между примерно 7 000 и примерно 13 000. Это неудивительно - напомним, что мы генерировали данные с таким поведением.
Предопределенные функции, предоставляемые с Statistics and Machine Learning 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, падает до более низких значений, затем снова повышается. Это типично для уровня опасности для компонента, который более восприимчив к отказу в начале своей жизни (младенческая смертность), и снова в более позднем возрасте (старение).
Также заметьте, что скорость опасности не может быть оценена выше самого большого наблюдения без цензуры для непараметрической модели, и график падает до нуля.
Для моделируемых данных, которые мы использовали в этом примере, мы обнаружили, что распределение Weibull не является подходящей подгонкой. Мы смогли хорошо подогнать данные к непараметрической подгонке, но эта модель была полезна только в области значений данных.
Другой альтернативой является использование другого параметрического распределения. Statistics and Machine Learning Toolbox включает функции для других распространенных распределений в течение жизни, таких как lognormal, gamma и Birnbaum-Saunders, а также многих других распределений, которые обычно не используются в моделях в течение жизни. Можно также задать и подгонять пользовательские параметрические модели к пожизненным данным, как описано в примере Fitting Custom Univariate Distributions, Part 2.
Другой альтернативой было бы использование смеси двух параметрических распределений - одного, представляющего ранний отказ, и другого, представляющего остальную часть распределения. Подгонка смесей распределений описана в примере Fitting Custom Univariate Distributions.