Вычислите оценку Каплана-Мейера эмпирической кумулятивной функции распределения (cdf) для симулированных данных о выживании.

Сгенерируйте данные о выживании из распределения Weibull параметрами 3 и 1.

rng('default')  % For reproducibility
failuretime = random('wbl',3,1,15,1);

Вычислите оценку Каплана-Мейера эмпирического cdf для данных о выживании.

[f,x] = ecdf(failuretime);
[f,x]

ans = 16×2

         0    0.0895
    0.0667    0.0895
    0.1333    0.1072
    0.2000    0.1303
    0.2667    0.1313
    0.3333    0.2718
    0.4000    0.2968
    0.4667    0.6147
    0.5333    0.6684
    0.6000    1.3749
      ⋮

Постройте предполагаемый эмпирический cdf.

ecdf(failuretime)

Сгенерируйте подвергнутые цензуре правом данные о выживании и сравните эмпирическую кумулятивную функцию распределения (cdf) с известным cdf.

Сгенерируйте времена отказа от экспоненциального распределения со средним временем отказа 15.

rng('default')  % For reproducibility
y = exprnd(15,75,1);

Сгенерируйте времена уволенного от экспоненциального распределения со средним временем отказа 30.

d = exprnd(30,75,1);

Сгенерируйте наблюдаемые времена отказа, то есть, минимум сгенерированных времен отказа и времен уволенного.

t = min(y,d);

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

censored = (y>d);

Вычислите эмпирический cdf и доверительные границы.

[f,x,flo,fup] = ecdf(t,'Censoring',censored);

Постройте эмпирический cdf и доверительные границы.

ecdf(t,'Censoring',censored,'Bounds','on')
hold on

Наложите график известного населения cdf.

xx = 0:.1:max(t);
yy = 1-exp(-xx/15);
plot(xx,yy,'g-','LineWidth',2)
axis([0 max(t) 0 1])
legend('Empirical cdf','Lower confidence bound', ...
    'Upper confidence bound','Known population cdf', ...
    'Location','southeast')
hold off

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

Сгенерируйте пожизненные данные из распределения Weibull параметрами 100 и 2.

rng('default')  % For reproducibility
R = wblrnd(100,2,100,1);

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

ecdf(R,'Function','survivor','Alpha',0.01,'Bounds','on')
hold on

Наложите график функции оставшегося в живых Weibull.

x = 1:1:250;
wblsurv = 1-cdf('weibull',x,100,2);
plot(x,wblsurv,'g-','LineWidth',2)
legend('Empirical survivor function','Lower confidence bound', ...
    'Upper confidence bound','Weibull survivor function', ...
    'Location','northeast')

Функция оставшегося в живых Weibull на основе фактического распределения в доверительных границах.

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

Сгенерируйте времена отказа от распределения Бирнбаума-Сондерса.

rng('default')  % For reproducibility
failuretime = random('BirnbaumSaunders',0.3,1,[100,1]);

Примите, что исследование запускается во время 0.1 и концы во время 0.9. Предположение подразумевает, что времена отказа, меньше чем 0,1 оставляют подвергнутыми цензуре, и времена отказа, больше, чем 0,9, правильные подвергнутый цензуре.

Создайте вектор, в котором каждый элемент указывает на состояние цензуры соответствующего наблюдения в failuretime. Используйте –1, 1, и 0, чтобы указать на лево-подвергнутые цензуре, подвергнутые цензуре правом, и полностью наблюдаемые наблюдения, соответственно.

L = 0.1;
U = 0.9;
left_censored = (failuretime<L);
right_censored = (failuretime>U);
c = right_censored - left_censored;

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

ecdf(failuretime,'Function','cumulative hazard', ...
    'Censoring',c,'Bounds','on')

Вычислите и постройте эмпирический cdf подвергнутых цензуре интервалом данных.

Загрузите cities набор данных. Данные включают оценки для девяти различных индикаторов качества жизни в 329 городах США: климат, корпус, здоровье, преступление, транспортировка, образование, искусства, воссоздание и экономика. Для каждого индикатора более высокая оценка лучше.

load cities

Выберите первый индикатор (климат) как выборочные данные.

Y = ratings(:,1);

Примите что индикаторы в Y значения, округленные до ближайшего целого числа. Затем можно обработать значения в Y как подвергнутые цензуре интервалом наблюдения. Наблюдение y в Y указывает, что фактическая оценка между y–0.5 и y+0.5.

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

intervalY = [Y-0.5, Y+0.5];

Вычислите эмпирические cdf значения.

[f,x] = ecdf(intervalY);

Постройте эмпирические cdf значения.

figure
ecdf(intervalY)

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

idx_roi = 21:30;
xlim([x(idx_roi(1),1) x(idx_roi(end),2)])

Отобразите соответствующий x и f значения.

table(idx_roi',x(idx_roi,:),f(idx_roi,:), ...
    'VariableNames',{'Index','x','Empirical cdf F(x)'})

ans=10×3 table
    Index          x           Empirical cdf F(x)
    _____    ______________    __________________

     21      377.5    378.5         0.069909     
     22      382.5    383.5         0.075988     
     23      384.5    385.5         0.079027     
     24      390.5    391.5         0.082067     
     25      395.5    396.5         0.085106     
     26      397.5    398.5         0.091185     
     27      400.5    401.5         0.094225     
     28      401.5    402.5         0.097264     
     29      403.5    404.5          0.10334     
     30      409.5    410.5          0.10638     

Закрашенные прямоугольники указывают на изменение эмпирических cdf значений F (x) в соответствующих интервалах. Например, второй закрашенный прямоугольник слева в масштабируемом графике соответствует интервалу (382.5,383.5]. F (382.5) 0.075988, F (383.5) 0.079027, и изменение от 0,075988 до 0,079027 происходит в интервале (382.5,383.5]. Точная синхронизация изменения сомнительна.

Можно построить оценки интервала по-разному. Если вы принимаете, что изменение вероятности происходит в начале каждого интервала, можно построить F (x) значения с помощью первого столбца x.

figure
stairs(x(:,1),f)
title('Probability changes at the start')
xlabel('x')
ylabel('F(x)')
xlim([x(idx_roi(1),1) x(idx_roi(end),2)])

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

figure
stairs(x(:,2),f)
title('Probability changes at the end')
xlabel('x')
ylabel('F(x)')
xlim([x(idx_roi(1),1) x(idx_roi(end),2)])

Постройте предыдущие два графика, чтобы визуализировать интервалы.

figure
stairs(x(:,1),f)
hold on
stairs(x(:,2),f)
title('Probability changes in the interval')
xlabel('x')
ylabel('F(x)')
xlim([x(idx_roi(1),1) x(idx_roi(end),2)])
hold off

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

Загрузите выборочные данные. Визуализируйте терпеливые данные о весе с помощью гистограммы.

load patients
histogram(Weight(strcmp(Gender,'Female')))
hold on
histogram(Weight(strcmp(Gender,'Male')))
legend('Female','Male')

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

Вычислите эмпирический cdf для данных.

[f,x] = ecdf(Weight);

Создайте кусочную линейную аппроксимацию к эмпирическому cdf путем принятия значение каждые пять точек.

f = f(1:5:end);
x = x(1:5:end);

Постройте эмпирический cdf и приближение.

figure
ecdf(Weight)
hold on
plot(x,f,'ko-','MarkerFace','r') 
legend('Empirical cdf','Piecewise linear approximation', ...
    'Location','best')

Создайте кусочный линейный объект вероятностного распределения использование кусочного приближения эмпирического cdf.

pd = makedist('PiecewiseLinear','x',x,'Fx',f)

pd = 
  PiecewiseLinearDistribution

F(111) = 0
F(118) = 0.05
F(124) = 0.13
F(130) = 0.25
F(135) = 0.37
F(142) = 0.5
F(163) = 0.55
F(171) = 0.61
F(178) = 0.7
F(183) = 0.82
F(189) = 0.94
F(202) = 1

Сгенерируйте 100 случайных чисел от распределения.

rng('default') % For reproducibility
rw = random(pd,[100,1]);

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

figure
histogram(Weight)
hold on
histogram(rw)
legend('Original data','Generated data')

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