Шаг 1. Загрузка и подготовка образцов данных.

Загрузите образцы данных.

load('lightbulb.mat')

Первый столбец данных имеет срок службы (в часах) двух типов лампочек. Второй столбец содержит информацию о типе лампочки. 0 - флуоресцентные лампы, а 1 - лампа накаливания. Третий столбец содержит цензурную информацию. 1 указывает цензурные данные, а 0 - точное время сбоя. Это смоделированные данные.

Создайте переменную для каждого типа лампочки, а также включите информацию о цензуре.

fluo = [lightbulb(lightbulb(:,2)==0,1),...
			lightbulb(lightbulb(:,2)==0,3)];
insc = [lightbulb(lightbulb(:,2)==1,1),...
			lightbulb(lightbulb(:,2)==1,3)];

Шаг 2. Постройте график оценочных функций выживших.

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

figure()
[f,x,flow,fup] = ecdf(fluo(:,1),'censoring',fluo(:,2),...
				'function','survivor');
ax1 = stairs(x,f);
hold on
stairs(x,flow,':')
stairs(x,fup,':')
[f,x,flow,fup] = ecdf(insc(:,1),'censoring',insc(:,2),...
				'function','survivor');
ax2 = stairs(x,f,'color','r');
stairs(x,flow,':r')
stairs(x,fup,':r')
legend([ax1,ax2],{'Fluorescent','Incandescent'})
xlabel('Lifetime (hours)')
ylabel('Survival probability')

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

Шаг 3. Подогнать распределение Burr типа XII.

Соответствие распределения Burr данным о времени жизни люминесцентных ламп и ламп накаливания.

pd = fitdist(fluo(:,1),'burr','Censoring',fluo(:,2))

pd = 
  BurrDistribution

  Burr distribution
    alpha = 29143.4   [0.903922, 9.39617e+08]
        c = 3.44582   [2.13013, 5.57417]
        k = 33.7039   [8.10737e-14, 1.40114e+16]

pd2 = fitdist(insc(:,1),'burr','Censoring',insc(:,2))

pd2 = 
  BurrDistribution

  Burr distribution
    alpha = 2650.76   [430.773, 16311.4]
        c = 3.41898   [2.16794, 5.39197]
        k =  4.5891   [0.0307809, 684.185]

Наложить функции выжившего Бёрра типа XII.

ax3 = plot(0:500:15000,1-cdf('burr',0:500:15000,29143.5,...
			3.44582,33.704),'m');
ax4 = plot(0:500:5000,1-cdf('burr',0:500:5000,2650.76,...
			3.41898,4.5891),'g');
legend([ax1;ax2;ax3;ax4],'Festimate','Iestimate','FBurr','IBurr')

Распределение Burr обеспечивает хорошую пригодность для срока службы лампочек в этом примере.

Шаг 4. Подберите модель пропорциональных рисков Кокса.

Установите пропорциональную регрессию рисков Кокса, где тип колбы является пояснительной переменной.

[b,logl,H,stats] = coxphfit(lightbulb(:,2),lightbulb(:,1),...
'Censoring',lightbulb(:,3));
stats

stats = struct with fields:
                    covb: 1.0757
                    beta: 4.7262
                      se: 1.0372
                       z: 4.5568
                       p: 5.1936e-06
                   csres: [100x1 double]
                  devres: [100x1 double]
                 martres: [100x1 double]
                  schres: [100x1 double]
                 sschres: [100x1 double]
                  scores: [100x1 double]
                 sscores: [100x1 double]
    LikelihoodRatioTestP: 0

Значение $$$$p, pуказывает, что тип лампочки является статистически значимым. Оценка коэффициента опасности exp $$$$$$($$b) = 112,8646. Это означает, что опасность для ламп накаливания в 112,86 раза превышает опасность для люминесцентных ламп .