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

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

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

Импортируем необходимые пакеты:

import Pkg
Pkg.add(["Distributions", "Random", "Plots", "Statistics", "StatsBase", "StatsPlots", "Optim", "KernelDensity", "MATLAB"])
using Distributions, Random, Plots, Statistics, StatsBase, StatsPlots, Optim, KernelDensity, MATLAB

Особые свойства данных о времени жизни

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

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

Random.seed!(2)
lifetime = [rand(Weibull(3, 15000), 90); rand(Weibull(3, 1500), 10)]

В данном примере предположим, что мы тестируем дроссельные заслонки в стрессовых условиях, так что каждый час испытаний эквивалентен 100 часам реальной эксплуатации. По практическим соображениям испытания на надёжность часто останавливают по истечении фиксированного времени. В этом примере мы используем 140 часов, что эквивалентно общему сроку реальной службы в 14000 часов. Одни образцы выходят из строя во время испытаний, а другие выдерживают все 140 часов. В реальном испытании для последних было бы зафиксировано время 14000 часов, что мы и воспроизводим в смоделированных данных. Также общепринятой практикой является сортировка времени отказов.

T = 14000
obstime = sort(min.(T, lifetime))

Мы знаем, что любые дроссели, прошедшие испытание, в конечном счёте выйдут из строя, но тест недостаточно продолжителен, чтобы наблюдать их фактическое время до отказа. Их времена работы известны лишь как больше 14000 часов. Эти значения называют цензурированными.

failed = obstime[obstime .< T]
nfailed = length(failed)
survived = obstime[obstime .== T]
nsurvived = length(survived)
censored = obstime .>= T
p1 = plot()
for i in 1:length(obstime)
    plot!([0, obstime[i]], [i, i], color=:blue, line=:solid, label="")
end

if nsurvived > 0
    for i in 1:nsurvived
        plot!([T, 3e4], [nfailed + i, nfailed + i], color=:blue, line=:dot, label="")
    end
end
plot!([T, T], [0, nfailed + nsurvived], color=:black, line=:solid, label="")
annotate!(T + 500, 30, text("<-- Неизвестное время жизни", 10, :left))
xlabel!("Время жизни")
ylabel!("Количество наблюдений")
xlims!(0, 3e4)
display(p1)

На данной диаграмме видно, что около 40% наших данных являются цензурированными при значении 14000 часов.

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

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

Функция плотности вероятности указывает относительную вероятность отказа в разные моменты времени.
Функция выживания показывает вероятность выживания как функцию времени и является просто 1 минус кумулятивная функция распределения.
Интенсивность отказов даёт мгновенную вероятность отказа при условии выживания до данного момента. Она определяется как делённая на функцию выживания. В этом примере интенсивность отказов оказывается возрастающей, что означает — изделия становятся более склонны к отказу с течением времени (старение).
Вероятностный график — это перенормированная кумулятивная функция распределения и используется для сравнения данных с подогнанным распределением.

x = range(1, 30000, length=100)

p2 = plot(layout=(2,2), size=(800,600))
plot!(p2[1], x, pdf.(Weibull(2, 14000), x), label=nothing, linewidth=2)
plot!(p2[1], x, pdf.(Weibull(2, 18000), x), label=nothing, linewidth=2)
plot!(p2[1], x, pdf.(Weibull(1.1, 14000), x), label=nothing, linewidth=2)
title!(p2[1], "Плотность вероятности")

plot!(p2[2], x, ccdf.(Weibull(2, 14000), x), label=nothing, linewidth=2)
plot!(p2[2], x, ccdf.(Weibull(2, 18000), x), label=nothing, linewidth=2)
plot!(p2[2], x, ccdf.(Weibull(1.1, 14000), x), label=nothing, linewidth=2)
title!(p2[2], "Выживание")

wblhaz(x, a, b) = pdf(Weibull(b, a), x) / ccdf(Weibull(b, a), x)
plot!(p2[3], x, wblhaz.(x, 14000, 2), label=nothing, linewidth=2)
plot!(p2[3], x, wblhaz.(x, 18000, 2), label=nothing, linewidth=2)
plot!(p2[3], x, wblhaz.(x, 14000, 1.1), label=nothing, linewidth=2)
title!(p2[3], "Интенсивность отказов")

weibull_data = rand(Weibull(2, 14000), 40)
sorted_data = sort(weibull_data)
n = length(sorted_data)
empirical_probs = (1:n) ./ (n + 1)
theoretical_quantiles = quantile.(Weibull(2, 14000), empirical_probs)
scatter!(p2[4], theoretical_quantiles, sorted_data, marker=:xcross, markercolor=:blue, markersize=2, label=nothing)
plot!(p2[4], theoretical_quantiles, theoretical_quantiles, linecolor=:black, linestyle=:dash, linewidth=1.5, label=nothing)
title!(p2[4], "Вероятностный график")
#xlabel!(p2[4], "Данные")
#ylabel!(p2[4], "Вероятность")
plot!(p2[4], xticks=[], yticks=[]) 
display(p2)

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

Подгонка распределения Вейбулла

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

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

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

p3 = plot()
function ecdf_with_ci(times, censored; alpha=0.05)
    sorted_indices = sortperm(times)
    sorted_times = times[sorted_indices]
    sorted_censored = censored[sorted_indices]
    n = length(times)
    F = zeros(n+1)
    F_lo = zeros(n+1)
    F_up = zeros(n+1)
    x = zeros(n+1)
    F[1] = 0.0
    F_lo[1] = 0.0
    F_up[1] = 0.0
    x[1] = 0.0
    at_risk = n 
    for i in 1:n
        if !sorted_censored[i] 
            survival_prob = (at_risk - 1) / at_risk
            F[i+1] = 1 - (1 - F[i]) * survival_prob
            se = sqrt(F[i+1] * (1 - F[i+1]) / at_risk) 
            z = quantile(Normal(), 1 - alpha/2)
            F_lo[i+1] = max(0, F[i+1] - z * se)
            F_up[i+1] = min(1, F[i+1] + z * se)
        else 
            F[i+1] = F[i]
            F_lo[i+1] = F_lo[i]
            F_up[i+1] = F_up[i]
        end
        at_risk -= 1
        x[i+1] = sorted_times[i]
    end
    return x, F, F_lo, F_up
end

x_ecdf, empF, empFlo, empFup = ecdf_with_ci(obstime, censored)

function create_stairs(x, y)
    x_stairs = Float64[]
    y_stairs = Float64[]
    
    for i in 1:length(x)
        if i > 1
            push!(x_stairs, x[i])
            push!(y_stairs, y[i-1])
        end
        push!(x_stairs, x[i])
        push!(y_stairs, y[i])
    end
    
    return x_stairs, y_stairs
end

x_stairs, F_stairs = create_stairs(x_ecdf, empF)
x_stairs_lo, F_stairs_lo = create_stairs(x_ecdf, empFlo)
x_stairs_up, F_stairs_up = create_stairs(x_ecdf, empFup)
plot!(x_stairs, F_stairs, color=:blue, linewidth=1.5, label=nothing)
plot!(x_stairs_lo, F_stairs_lo, color=:red, linestyle=:dot, linewidth=1, label=nothing)
plot!(x_stairs_up, F_stairs_up, color=:orange, linestyle=:dot, linewidth=1, label=nothing)

xlabel!("Время")
ylabel!("Соотношение отказов")
title!("Эмпирическая кумулятивная функция распределения")
display(p3)

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

Распределение Вейбулла часто является хорошей моделью надёжности оборудования. Функция weibull_fit подгоняет распределение Вейбулла к данным, включая данные с цензурированием. После оценки параметров мы оценим кумулятивную функцию распределения для подогнанной модели Вейбулла, используя эти оценки. Поскольку значения кумулятивной функции распределения основаны на оценённых параметрах, мы также вычислим доверительные границы для них.

function weibull_fit(data, censored)
    function weibull_loglik(θ, data, censored)
        scale, shape = θ
        if scale <= 0 || shape <= 0
            return Inf
        end
        loglik = 0.0
        for i in 1:length(data)
            if !censored[i]
                loglik += logpdf(Weibull(shape, scale), data[i])
            else
                loglik += logccdf(Weibull(shape, scale), data[i])
            end
        end
        return -loglik
    end
    
    uncensored_data = data[.!censored]
    if isempty(uncensored_data)
        initial_scale = mean(data)
    else
        initial_scale = mean(uncensored_data)
    end
    initial_shape = 2.0
    initial_guess = [initial_scale, initial_shape]
    lower = [1.0, 0.1] 
    upper = [1e6, 100.0] 
    
    result = optimize(θ -> weibull_loglik(θ, data, censored), lower, upper, initial_guess, Fminbox(LBFGS()))
    return result.minimizer
end

function weibull_cdf_ci(x, scale, shape, data, censored; alpha=0.05)
    wblF = cdf.(Weibull(shape, scale), x)
    n_uncensored = sum(.!censored)
    z = quantile(Normal(), 1 - alpha/2)
    wblFlo = similar(wblF)
    wblFup = similar(wblF)
    
    for i in 1:length(x)
        se = z * sqrt(wblF[i] * (1 - wblF[i]) / n_uncensored)
        wblFlo[i] = max(0.0, wblF[i] - se)
        wblFup[i] = min(1.0, wblF[i] + se)
    end
    
    return wblF, wblFlo, wblFup
end

paramEsts = weibull_fit(obstime, censored)
xx = range(1, 2*T, length=500)
wblF, wblFlo, wblFup = weibull_cdf_ci(xx, paramEsts[1], paramEsts[2], obstime, censored)

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

p4 = plot()
plot!(x_stairs, F_stairs, color=:blue, linewidth=1.5, label=nothing)
plot!(xx, wblF, color=:red, linewidth=2, label=nothing)
plot!(xx, wblFlo, color=:red, linewidth=1, linestyle=:dash, label=nothing)
plot!(xx, wblFup, color=:red, linewidth=1, linestyle=:dash, label=nothing)
xlabel!("Время")
ylabel!("Вероятность отказа")
title!("Модели: Эмпирическая и Вейбулла")
xlims!(0, 2*T)
ylims!(0, 1)
display(p4)

Обратите внимание: модель Вейбулла позволяет делать прогнозы и вычислять вероятность отказа за пределами конца испытания. Однако видно, что подогнанная кривая не очень хорошо совпадает с нашими данными. У нас слишком много ранних отказов до 2000 часов по сравнению с тем, что предсказывала модель Вейбулла, и в результате меньше отказов, чем ожидалось, для промежутка примерно от 7000 до 13000 часов.

Это неудивительно — напомним, что мы сгенерировали данные именно с таким поведением.

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

Мы можем построить сглаженную непараметрическую кривую через эмпирическую с помощью функции kde. Мы удалим доверительные и добавим две кривые: одну с параметром сглаживания по умолчанию, и одну с параметром сглаживания, равным 1/3 от значения по умолчанию. Меньший параметр сглаживания позволяет кривой точнее следовать данным.

if !@isdefined(empF) || !@isdefined(x_ecdf)
    empirical_cdf = ecdf(obstime)
    x_ecdf = range(0, maximum(obstime), length=100)
    empF = empirical_cdf.(x_ecdf)
end

kde_result = kde(obstime)
u = KernelDensity.default_bandwidth(obstime)
kde_standard = kde(obstime, bandwidth = u)
kde_small = kde(obstime, bandwidth = u/3)
npF_pdf = pdf(kde_standard, xx)
npF3_pdf = pdf(kde_small, xx)
npF_cdf = cumsum(npF_pdf) .* step(xx)
npF_cdf = npF_cdf ./ maximum(npF_cdf) 
npF3_cdf = cumsum(npF3_pdf) .* step(xx)
npF3_cdf = npF3_cdf ./ maximum(npF3_cdf)

p5 = plot()
x_stairs, F_stairs = create_stairs(x_ecdf, empF)
plot!(x_stairs, F_stairs,  color=:blue,  linewidth=1.5, label="Эмпирическая")
plot!(xx, wblF, color=:red, linewidth=2, label="Вейбулла")
plot!(xx, npF_cdf, color=:green, linewidth=2, label="Непараметрическая, стандартная")
plot!(xx, npF3_cdf, color=:magenta, linewidth=2, label="Непараметрическая, 1/3 от стандартной")
xlims!(0, 1.3*T)
title!("Модели")
xlabel!("Время")
ylabel!("Вероятность")
display(p5)

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

Давайте вычислим интенсивность отказов для этой непараметрической модели и построим её график в диапазоне имеющихся данных.

hzr = npF3_pdf ./ (1 .- npF3_cdf)
hzr[isinf.(hazrate)] .= 0
hzr[isnan.(hazrate)] .= 0
p6 = plot(xx, hzr, linewidth = 2, legend = false)
title!("Интенсивность отказов непараметрической модели")
xlims!(0, T)
ylims!(0, 0.00017)
display(p6)

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

Сравнение с Matlab

Для наглядного сравнения, давайте построим сглаженную непараметрическую кривую, используя вместо функции kde Julia, функцию ksdensity Matlab.

@mput obstime censored xx
mat"""
[npF,ignore,u] = ksdensity(obstime, xx, 'cens', censored, 'function', 'cdf');
npF3 = ksdensity(obstime, xx, 'cens', censored, 'function', 'cdf', 'width', u/3);
hz = ksdensity(obstime, xx, 'cens', censored, 'width', u/3);
"""
@mget npF npF3 hz

p7 = plot()
x_stairs, F_stairs = create_stairs(x_ecdf, empF)
plot!(x_stairs, F_stairs,  color=:blue,  linewidth=1.5, label="Эмпирическая")
plot!(xx, wblF, color=:red, linewidth=2, label="Вейбулла")
plot!(xx, npF, color=:green, linewidth=2, label="Непараметрическая, стандартная")
plot!(xx, npF3, color=:magenta, linewidth=2, label="Непараметрическая, 1/3 от стандартной")
xlims!(0, 1.3*T)
title!("Модели Matlab")
xlabel!("Время")
ylabel!("Вероятность")
display(p7)

Вычислим интенсивность отказов для непараметрической модели и построим её график на основе вычислений в Matlab.

hazrate = hz ./ (1 .- npF3)
p8 = plot(xx, hazrate, linewidth = 2, legend = false)
title!("Интенсивность отказов непараметрической модели Matlab")
xlims!(0, T)
display(p8)

Выводы:

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

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

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

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

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

Тип

Краткое описание

Связанные материалы

Теги

Категории

Языки

Форматы

Уровень

Статус

Опубликовано

Обновлено

Источник