Шаг 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 Type 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]

Наложение функций выжившего Burr типа 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 раза превышает опасность для лампочек флуоресценции.