Fordeling til sviktintensitet - AssetLife verktøykasse

Formålet med denne notebooken er å finne fordelingen på sviktintensiteten til et datasett. For å gjøre det kan man plotte estimatoren for kumulativ sviktsannsynlighet og se på formen til grafen.

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from reliability.Nonparametric import NelsonAalen
from reliability.Fitters import Fit_Weibull_2P, Fit_Exponential_1P, Fit_Gamma_2P

# les og formatter datasettet
df = pd.read_excel("ttt_data.xlsx")

failures = df.loc[df["Havarert"] == 1, "Alder"].values
right_censored = df.loc[df["Havarert"] == 0,"Alder"].values

# Forberede datasett på lineærapproksimasjon
NA = NelsonAalen(failures=failures, right_censored=right_censored, label='Nelson-Aalen', print_results=False, show_plot=False, plot_type='CHF')

# Renser dataen slik at vi bare har én H-verdi per t
df = pd.DataFrame({'t': NA.xvals, 'H': NA.CHF})
df = df[(df['t'] > 0) & (df['H'] > 0) & np.isfinite(df['H'])]
df = df.groupby('t', as_index=False).last()

t_f = df['t'].to_numpy()
H_f = df['H'].to_numpy()

print("Antall unike punkter:", len(t_f))
print(df.head(10))

Antall unike punkter: 24
    t         H
0  21 0.0232558
1  25 0.0232558
2  36 0.0765986
3  38  0.106902
4  40  0.106902
5  41  0.106902
6  42  0.143939
7  45  0.185605
8  46  0.185605
9  47   0.23106

Nelson-Aalen¶

Her plottes Nelson-Aalen-estimatoren for kumulativ sviktintensitet $$ \hat{Z}(t)¶

\sum_{t_i \le t} \left( \frac{d_i}{n_i} \right) $$

Hvis den kan tilpasses med en rett linje vil det være en indikasjon på at sviktintensiteten kan tilpasses med en eksponensialfordeling som har konstant sviktsannsynlighet lik stigningstallet til grafen.

Hvis den plottes på log-log skala og kan tilpasses med en rett linje vil det være indikasjon på at sviktintensiteten kan tilpasses med en weibullfordeling.

# Eksponensialfordelingen vil være lineær
lam = np.sum(t_f * H_f) / np.sum(t_f**2)

NA = NelsonAalen(failures=failures, right_censored=right_censored, label='Nelson-Aalen', print_results=False, plot_type='CHF')
ax = plt.gca()
for coll in ax.collections:          
    if coll.get_label().startswith('_'):
        coll.set_label('95% KI')
plt.plot(t_f, lam * t_f, '-', label=f'Eksponensial fit (λ={lam:.4g})')
plt.xlabel('Tid')
plt.ylabel('CHF')
plt.legend()
plt.title('Lineær skala (eksponensialfordeling)')
plt.xlim(t_f.min() * 0.8, t_f.max() * 1.1)
plt.ylim(H_f.min() * 0.5, H_f.max() * 1.5)
plt.show()

#Weibullfordelingen vil være lineær i log-log
beta, intercept = np.polyfit(np.log(t_f), np.log(H_f), 1)   # intercept = -beta*ln(alpha)
alpha = np.exp(-intercept / beta)

NA = NelsonAalen(failures=failures, right_censored=right_censored, label='Nelson-Aalen', print_results=False, plot_type='CHF')
t_line = np.linspace(t_f.min(), t_f.max(), 200)
plt.plot(t_line, (t_line / alpha) ** beta, '-',
         label=f'Weibull fit (β={beta:.3g}, α={alpha:.3g})')

ax = plt.gca()
for coll in ax.collections:          
    if coll.get_label().startswith('_'):
        coll.set_label('95% KI')

plt.xlabel('t')
plt.ylabel('CHF')
plt.xscale('log')
plt.yscale('log')
plt.xlim(t_f.min() * 0.8, t_f.max() * 1.2)
plt.ylim(H_f.min() * 0.5, H_f.max() * 1.5)
plt.legend()
plt.title('Log-log skala (weibullfordeling)')
plt.show()

Sammenligning¶

Modellsammenligning med AICc

For å avgjøre om Weibull- eller eksponensialfordelingen passer best til dataene, sammenligner vi de to modellenes Akaike-informasjonskriterium med liten-utvalgs-korreksjon (AICc):

\text{AICc} = -2\ell(\hat\theta) + 2k + \frac{2k^2+2k}{n-k-1}

(1)

hvor $\ell(\hat\theta)$ er den maksimerte log-likelihoodverdien, $k$ er antall estimerte parametre (1 for eksponensial, 2 for Weibull), og $n$ er utvalgsstørrelsen. AICc balanserer modellens tilpasning til dataene mot dens kompleksitet, slik at en modell med flere parametre (Weibull) Får et ekstra stort straffeledd, og dermed ikke automatisk foretrekkes over de med færre parametre. Vi tolker differansen $\Delta\text{AICc} = \text{AICc}_{\text{Exp}} - \text{AICc}_{\text{Weibull}}$ etter standard tommelfingerregler: under 2 gir ingen reell støtte for weilbull over eksponensialfordeling, 2–7 gir moderat støtte, og over 7–10 gir sterk støtte for at weibull bør foretrekkes og at sviktraten endres med alder.

def modelltilpasning(failures, right_censored):
    fit_exp = Fit_Exponential_1P(failures=failures, right_censored=right_censored, show_probability_plot=False, print_results=False)
    fit_weib = Fit_Weibull_2P(failures=failures, right_censored=right_censored, show_probability_plot=False, print_results=False)
    delta_AICc = fit_exp.AICc - fit_weib.AICc 
    if delta_AICc < 2:
        print(f'''Differansen mellom AICc for eksponensialfordelingen og AICc for Weibullfordelingen er {delta_AICc}.
Ettersom dette er mindre enn 2 er det ikke statistisk signfikant grunnlag til å si at sviktraten øker med alder.''')
    elif delta_AICc < 7:
        print(f'''Differansen mellom AICc for eksponensialfordelingen og AICc for Weibullfordelingen er {delta_AICc}.
Ettersom dette er mellom 2 og 7 gir det svak til moderat støtte for å si at sviktraten øker med alder.''')
    else:  # delta_AICc >= 7
        print(f'''Differansen mellom AICc for eksponensialfordelingen og AICc for weibullfordelingen er {delta_AICc}.
Dette gir sterk støtte for å si at sviktraten øker med alder.''')

modelltilpasning(failures, right_censored)

Differansen mellom AICc for eksponensialfordelingen og AICc for weibullfordelingen er 10.038254824502872.
Dette gir sterk støtte for å si at sviktraten øker med alder.