Sannsynlighetsmaksimeringsestimator* - AssetLife verktøykasse

Dette kapittelet ser nærmere på sannsynlighetsmaksimeringsestimator (SME) for estimering av sviktsannsynlighet ved bruk av sensurert tilstandshistorikk og/eller ekspertvurderinger. Vi fokuserer på oppholdstid i TK1 og antar gammafordeling. SME brukes her for å estimere forventningsverdi $\mu$ og scale parameter $\theta$ . Metoden kan brukes på forskjellige typer data. Her demonstreres følgende input data:

Intervallsensurert oppholdstid. Intervallsensurert oppholdstid kan regnes ut fra tidsserier av tilstandskarakterer som vist her.
Ekspertvurderinger i form av persentiler P $_{10}$ og P $_{50}$ . Se forklaring her.
Kombinasjon av disse to datatypene.

Metoden beskrevet under er basert på Welte, T. (2008) Estimation of Sojourn Time Distribution Parameters Based on Expert Opinion and Condition Monitoring Data

* Maximum Likelihood Estimator (MLE) på engelsk

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gamma
from scipy.optimize import minimize
%config InlineBackend.figure_format = 'retina'
import matplotlib as mpl
mpl.rcParams.update({
    "figure.facecolor": "none",
    "axes.facecolor": "none",
    "savefig.facecolor": "none",})

Fontconfig error: No writable cache directories

Intervallsensurert data som input¶

Intervallsensurert data for oppholdstid i TK1 er gitt under for 10 observerte komponenter. Her har, for eksempel, den første komponenten blitt observert i TK1 i minimum 15 år og maksimum 20 år.

# Eksempeldata for min og max oppholdstid i TK1
min_TK1 = [15, 25, 35, 50, 10, 45, 30, 30, 40, 25]
max_TK1 = [20, 30, 40, 55, 15, 50, 35, 35, 45, 30]

Loglikelihood data¶

En loglikelihood-funksjon svarer på følgende spørsmål: gitt en gammafordeling med parametre $\mu$ og $\theta$ , hvor sannsynlig er den observerte dataen (min_TK1 og max_TK1)? Funksjonen er definert under og baserer seg på differansen i den kumulative sannsynlighetskurven mellom henholdsvis minimum og maximum tid.

def negloglikelihood_data(par, xmin, xmax):
    mu, theta = par

    if mu <= 0 or theta <= 0:
        return np.inf

    shape = mu/theta

    diff = (
        gamma.cdf(xmax, a=shape, scale=theta)
        - gamma.cdf(xmin, a=shape, scale=theta)
    )

    diff = np.clip(diff, 1e-300, None)

    return -np.sum(np.log(diff))

Estimering av parametre¶

Ved å minimere den negative loglikelihooden kan man finne de mest sannsynlige verdiene for $\mu$ og $\theta$ gitt ved mu_star og theta_star under.

res = minimize(
    negloglikelihood_data,
    x0=[25, 5],
    args=(min_TK1, max_TK1),
    bounds=[(1e-6, None), (1e-6, None)],
    method="L-BFGS-B"
)
mu_star, theta_star = res.x

Plotting av resultater¶

Videre kan vi plotte loglikelihood som funksjon av $\mu$ og $\theta$ og vise 95% konfidensintervallet, for å se usikkerheten i mu_star og theta_star. Vi plotter også den resulterende gammafordelingen med 95% konfidensintervall.

LLmax = res.fun

mu_grid = np.linspace(0.25*mu_star, 3*mu_star, 300)
theta_grid = np.linspace(0.25*theta_star, 3*theta_star, 300)

MU, TH = np.meshgrid(mu_grid, theta_grid, indexing="ij")
data = np.empty_like(MU)

for i in range(MU.shape[0]):
    for j in range(MU.shape[1]):
        data[i, j] = negloglikelihood_data(
            [MU[i, j], TH[i, j]],
            min_TK1,
            max_TK1,
        )

fig, ax = plt.subplots(1, 1, figsize=(5, 4), layout='constrained')
im = ax.imshow(
    np.exp(-(data-data.min())),
    origin='lower',
    aspect='auto',
    cmap='viridis',
    extent=[theta_grid.min(), theta_grid.max(), mu_grid.min(), mu_grid.max()]
)
LLmax = data.min()
levels = [
    LLmax + 0.5 * 5.99,  # 95%
]

cs = ax.contour(
    theta_grid,
    mu_grid,
    data,
    levels=levels,
    colors=['w', 'r'],
    linewidths=2,
)

ax.clabel(cs, fmt={
    levels[0]: '95% CI'})
cbar = fig.colorbar(im, ax=ax);cbar.set_label('Likelihood norm')
ax.plot([theta_star],[mu_star],'*r',markersize=15,label = f'(theta,mu) = ({theta_star:.1f},{mu_star:.1f})')
ax.legend(facecolor="white");ax.set_title('SME data');ax.grid()
ax.set_xlabel('theta');ax.set_ylabel('mu')

# Relative likelihood on grid
weights = np.exp(-(data - data.min()))
weights /= weights.sum()
# Sample parameter pairs
N = 5_000

idx = np.random.choice(
    weights.size,
    size=N,
    p=weights.ravel()
)

i_mu, j_theta = np.unravel_index(idx, weights.shape)
mu_samples = mu_grid[i_mu]
theta_samples = theta_grid[j_theta]

t = np.linspace(0,120,100)
pdfs = np.array([
    gamma.pdf(
        t,
        a=m/th,
        scale=th
    )
    for m, th in zip(mu_samples, theta_samples)
])
pdfs = pdfs[np.all(np.isfinite(pdfs), axis=1)]

# Pointwise confidence band
pdf_lo = np.percentile(pdfs, 2.5, axis=0);
pdf_hi = np.percentile(pdfs, 97.5, axis=0);

fig, ax = plt.subplots(figsize=(5, 3))
ax.fill_between(
    t,
    pdf_lo,
    pdf_hi,
    alpha=0.3,
    label="95% CI"
)

ax.plot(
    t,
    gamma.pdf(
        t,
        a=mu_star/theta_star,
        scale=theta_star
    ),
    'k',
    lw=2,
    label=f'SME data'
)

ax.set_xlabel("Tid (år)");ax.set_ylabel("f(t)");ax.set_ylim([0,0.035]);
ax.grid();ax.set_title('Oppholdstid i TK1')
ax.legend()
plt.show()

Usikkerhet i SME og gammafordeling¶

95% konfidensintervall for parametrene i SME er regnet ut med Wilks’ teorem som sier at loglikelihood $l$ for en fordeling med $n=2$ parametre går mot $\chi^2$ når antall data går mot uendelig

2(l_{max}-l)\sim\chi_{95\%,n=2}^2 = 5.99.

(1)

95% konfidensintervall for den resulterende gammafordelingen er regnet ut ved å trekke et stort antall parametre som vektlegges med den relative likelihooden i forhold til $l_{max}$ . Deretter beregnes 2.5% og 97.5% persentiler per tidssteg.

Ekspertvurderinger som input¶

Ekspertvurderinger kan angis som persentiler for oppholdstiden i TK1. Det vil si, en eller flere eksperter oppgir tidspunkt for når de tror 10% og 50% har gått fra TK1 til TK2 (P $_{10}$ , P $_{50}$ ). Usikkerhet i disse anslagene, sigma, må også spesifiseres. Under ser vi et eksempel med tre eksperter:

P_10 = [10,15,15]
P_50 = [30,35,29]
sigma = [2,2,2]

Loglikelihood ekspert¶

En egen loglikelihood funksjon defineres for ekspertvurderinger. Her er det antatt at estimatene for persentiler er normalfordelte rundt den sanne verdien, gitt $\mu$ og $\theta$ , med standardavvik sigma.

def negloglikelihood_exp(par, P_10, P_50, sigma):
    mu, theta = par

    if mu <= 0 or theta <= 0:
        return np.inf

    shape = mu / theta

    P_10_model = gamma.ppf(0.10, a=shape, scale=theta)
    P_50_model = gamma.ppf(0.50, a=shape, scale=theta)

    P_10 = np.asarray(P_10)
    P_50 = np.asarray(P_50)
    sigma = np.asarray(sigma)

    return 0.5 * np.sum(
        ((P_10 - P_10_model) / sigma) ** 2
        + ((P_50 - P_50_model) / sigma) ** 2
    )

Plotting av resultater¶

res = minimize(
    negloglikelihood_exp,
    x0=[25, 5],
    args=(P_10, P_50,sigma),
    bounds=[(1e-6, None), (1e-6, None)],
    method="L-BFGS-B"
);
mu_star, theta_star = res.x;
LLmax = res.fun;

mu_grid = np.linspace(0.1*mu_star, 3*mu_star, 300);
theta_grid = np.linspace(0.1*theta_star, 3*theta_star, 300);

MU, TH = np.meshgrid(mu_grid, theta_grid, indexing="ij");
data = np.empty_like(MU);

for i in range(MU.shape[0]):
    for j in range(MU.shape[1]):
        data[i, j] = negloglikelihood_exp(
            [MU[i, j], TH[i, j]],
            P_10,
            P_50,
            sigma,
        );

fig, ax = plt.subplots(1, 1, figsize=(5, 4), layout='constrained')
im = ax.imshow(
    np.exp(-(data-data.min())),
    origin='lower',
    aspect='auto',
    cmap='viridis',
    extent=[theta_grid.min(), theta_grid.max(), mu_grid.min(), mu_grid.max()]
)

LLmax = data.min();
levels = [
    LLmax + 0.5 * 5.99,  # 95%
]

cs = ax.contour(
    theta_grid,
    mu_grid,
    data,
    levels=levels,
    colors=['w', 'r'],
    linewidths=2,
);

ax.clabel(cs, fmt={
    levels[0]: '95% CI'})
cbar = fig.colorbar(im, ax=ax);cbar.set_label('Likelihood norm')
ax.plot([theta_star],[mu_star],'*r',markersize=15,label = f'(theta,mu) = ({theta_star:.1f},{mu_star:.1f})')
ax.legend(facecolor="white");ax.set_title('SME eksperter');ax.grid()
ax.set_xlabel('theta');ax.set_ylabel('mu')

# Relative likelihood on grid
weights = np.exp(-(data - data.min()));
weights /= weights.sum();
# Sample parameter pairs
N = 5_000

idx = np.random.choice(
    weights.size,
    size=N,
    p=weights.ravel()
);

i_mu, j_theta = np.unravel_index(idx, weights.shape);
mu_samples = mu_grid[i_mu];
theta_samples = theta_grid[j_theta];

t = np.linspace(0,120,100)
pdfs = np.array([
    gamma.pdf(
        t,
        a=m/th,
        scale=th
    )
    for m, th in zip(mu_samples, theta_samples)
])
pdfs = pdfs[np.all(np.isfinite(pdfs), axis=1)]

# Pointwise confidence band
pdf_lo = np.percentile(pdfs, 2.5, axis=0);
pdf_hi = np.percentile(pdfs, 97.5, axis=0);

fig, ax = plt.subplots(figsize=(5, 3))
ax.fill_between(
    t,
    pdf_lo,
    pdf_hi,
    alpha=0.3,
    label="95% CI"
);

ax.plot(
    t,
    gamma.pdf(
        t,
        a=mu_star/theta_star,
        scale=theta_star
    ),
    'k',
    lw=2,
    label=f'SME eksperter'
);

ax.set_xlabel("Tid (år)");ax.set_ylabel("f(t)");ax.set_ylim([0,0.035]);
ax.grid();ax.set_title('Oppholdstid i TK1');
ax.legend();
plt.show();

Kombinasjon av data og ekspertvurderinger¶

Loglikelihood kan enkelt legges sammen for intervallsensurert data og ekspertvurderinger. For å vite om det er den intervallsensurerte dataen eller ekspertvurderingene som dominerer resultatet, ser vi på loglikelihood av disse hver for seg for de optimerte parametrene. MERK! Justering av sigma har mye å si for vekting av eksperter og intervallsensurert data.

min_TK1 = [15, 25, 35, 50, 10, 45, 30, 30, 40, 25]
max_TK1 = [20, 30, 40, 55, 15, 50, 35, 35, 45, 30]


P_10 = [10,15,15]
P_50 = [30,35,29]
sigma = [2,2,2]

def negloglikelihood_all(par, P_10, P_50, sigma, xmin, xmax):
    return negloglikelihood_exp(par, P_10, P_50, sigma)+negloglikelihood_data(par, xmin, xmax)

res = minimize(
    negloglikelihood_all,
    x0=[20, 5],
    args=(P_10, P_50, sigma, min_TK1, max_TK1),
    bounds=[(1e-6, None), (1e-6, None)]
)

mu_star, theta_star = res.x

print(
    "Expert logL:",
    negloglikelihood_exp(
        [mu_star, theta_star],
        P_10,
        P_50,
        sigma,
    )
)

print(
    "Data logL:",
    negloglikelihood_data(
        [mu_star, theta_star],
        min_TK1,
        max_TK1,
    )
)


mu_grid = np.linspace(0.1*mu_star, 4*mu_star, 100)
theta_grid = np.linspace(0.1*theta_star, 4*theta_star, 100)
MU, TH = np.meshgrid(mu_grid, theta_grid, indexing="ij")
data = np.empty_like(MU)

for i in range(MU.shape[0]):
    for j in range(MU.shape[1]):
        data[i, j] = negloglikelihood_all( [MU[i, j], TH[i, j]],
                                          P_10, 
                                          P_50, 
                                          sigma, 
                                          min_TK1, 
                                          max_TK1)

fig, ax = plt.subplots(1, 1, figsize=(5, 4), layout='constrained')
im = ax.imshow(np.exp(-(data-data.min())),
    origin='lower',
    aspect='auto',
    cmap='viridis',
    extent=[theta_grid.min(), theta_grid.max(), mu_grid.min(), mu_grid.max()]
)


LLmax = data.min()
levels = [
    LLmax + 0.5 * 5.99,  # 95%
]

cs = ax.contour(
    theta_grid,
    mu_grid,
    data,
    levels=levels,
    colors=['w', 'r'],
    linewidths=2,
)

ax.clabel(cs, fmt={
    levels[0]: '95% CI'})
cbar = fig.colorbar(im, ax=ax);cbar.set_label('Likelihood norm')
ax.plot([theta_star],[mu_star],'*r',markersize=15,label = f'(theta,mu) = ({theta_star:.1f},{mu_star:.1f})')
ax.legend(facecolor="white");ax.set_title('SME data+eksperter');ax.grid()
ax.set_xlabel('theta');ax.set_ylabel('mu')

# Relative likelihood on grid
weights = np.exp(-(data - data.min()))
weights /= weights.sum()
# Sample parameter pairs
N = 5_000

idx = np.random.choice(
    weights.size,
    size=N,
    p=weights.ravel()
)

i_mu, j_theta = np.unravel_index(idx, weights.shape)
mu_samples = mu_grid[i_mu]
theta_samples = theta_grid[j_theta]

t = np.linspace(0,120,100)
pdfs = np.array([
    gamma.pdf(
        t,
        a=m/th,
        scale=th
    )
    for m, th in zip(mu_samples, theta_samples)
])

pdfs = pdfs[np.all(np.isfinite(pdfs), axis=1)]

# Pointwise confidence band
pdf_lo = np.percentile(pdfs, 2.5, axis=0);
pdf_hi = np.percentile(pdfs, 97.5, axis=0);

fig, ax = plt.subplots(figsize=(5, 3))
ax.fill_between(
    t,
    pdf_lo,
    pdf_hi,
    alpha=0.3,
    label="95% CI"
)

ax.plot(
    t,
    gamma.pdf(
        t,
        a=mu_star/theta_star,
        scale=theta_star
    ),
    'k',
    lw=2,
    label=f'SME data+eksperter'
)

ax.set_xlabel("Tid (år)");ax.set_ylabel("f(t)");ax.set_ylim([0,0.035]);
ax.grid();ax.set_title('Oppholdstid i TK1')
ax.legend()
plt.show()

Expert logL: 4.792634804924227
Data logL: 24.06392905060703