生存時間解析をやってみる2(PyMC3)

masao-ds.hatenablog.com ベイズモデリングによる生存時間解析はこちらを参考にしてください。

1.やりたいこと/課題設定

ある工場ではとある製品を量産している。複数台の生産装置を使用して製品を量産しているが、装置の故障が多発し、製品を生産することができず機会損失が生じていた。そこで生産装置を別の故障時間が長いものに入れ替えることにした。ベイズモデリングによる生存時間により、入れ替え前と後の生産装置での1ヶ月稼働時(装置メンテナンスを1ヶ月と想定)の機会損失額の差を見積もりたい。

2.ベイズモデリングによる生存時間解析を用いた機会損失額の比較

2つの故障時間が異なる生産装置の故障時間データを作成し、ベイズモデリングによる生存時間解析を行い、故障時間をモデル化する。その後、1ヶ月間工場で納入している生産装置で生産を行う場合の機会損失額を2つの異なる生産装置で比較する。この際、工場への生産装置の納入数は50台、１日１台あたりの生産個数は100個、製品1個あたりの売上は100円とする。

2.1.データの生成

まず故障時間の異なる架空の故障データを作成する。

#旧装置
# ワイブル分布のパラメータ
shape_old = 2 # ワイブル分布の形状パラメータ
scale_old = 28.0  # ワイブル分布の尺度パラメータ

# ランダムな故障時間を生成
np.random.seed(42)  # 乱数のシードを設定（再現性のため）
failures_old = np.random.weibull(shape_old, size=100) * scale_old

#新装置
# ワイブル分布のパラメータ
shape_new = 3.0 # ワイブル分布の形状パラメータ
scale_new = 35.0  # ワイブル分布の尺度パラメータ

# ランダムな故障時間を生成
np.random.seed(42)  # 乱数のシードを設定（再現性のため）
failures_new = np.random.weibull(shape_new, size=100) * scale_new

plt.hist(failures_old, alpha=0.5)
plt.hist(failures_new, alpha=0.5)
plt.show()

変更後の生産装置の方が故障までの時間が長くなるように生成した。また、以下のように40日経過後のデータは打切りデータとして扱うことにした。

#40日よりも後を身観測データとする。
df_failures_old=pd.DataFrame(failures_old, columns=['Time'])
df_failures_old['censored'] = df_failures_old['Time'] >= 40
y_old = df_failures_old['Time'].values
censored_old = df_failures_old['censored'].values

df_failures_new=pd.DataFrame(failures_new, columns=['Time'])
df_failures_new['censored'] = df_failures_new['Time'] >= 40
y_new = df_failures_new['Time'].values
censored_new = df_failures_new['censored'].values

2.2.装置の生存時間解析

2.2.1.変更前の装置の生存時間解析

PyMC3を使用したベイズモデリングを行う。生存時間分布にはワイブル分布を仮定する。

def weibull_lccdf(x, alpha, beta):
    """ Log complementary cdf of Weibull distribution. """
    return -((x / beta) ** alpha)

#打ち切りなしデータ+打ち切りありデータ
with pm.Model() as model_2:
    alpha_sd = 10.0
    mu = pm.Normal("mu", mu=0, sigma=100)
    alpha_raw = pm.Normal("a0", mu=0, sigma=0.1)
    alpha = pm.Deterministic("alpha", tt.exp(alpha_sd * alpha_raw))
    beta = pm.Deterministic("beta", tt.exp(mu / alpha))

    y_obs = pm.Weibull("y_obs", alpha=alpha, beta=beta, observed=y_old[~censored_old])
    #事後対数尤度に制約としてつけ加える。この項を含めた事後対数尤度が最大になるようにMCMCでサンプリングする。
    y_cens = pm.Potential('y_cens', weibull_lccdf(y_old[censored_old], alpha, beta))

with model_2:
    trace_2 = pm.sample(target_accept=0.9, init="adapt_diag")

事後サンプリング結果をプロットする。

_ = pm.traceplot(trace_2)
_ = pm.plot_posterior(trace_2)

サマリーを確認。

pm.summary(trace_2)

    mean    sd  hdi_3%  hdi_97% mcse_mean   mcse_sd ess_bulk    ess_tail    r_hat
mu  5.626   0.497   4.644   6.535   0.020   0.014   637.0   651.0   1.01
a0  0.052   0.008   0.036   0.068   0.000   0.000   668.0   642.0   1.01
alpha   1.692   0.140   1.439   1.974   0.005   0.004   668.0   642.0   1.01
beta    27.813  1.812   24.422  31.111  0.034   0.024   2768.0  2268.0  1.00

2.2.2.変更後の装置の生存時間解析

同様にPyMC3を使用したベイズモデリングにより生存時間解析を行う。

with pm.Model() as model_3:
    alpha_sd = 10.0
    mu = pm.Normal("mu", mu=0, sigma=100)
    alpha_raw = pm.Normal("a0", mu=0, sigma=0.1)
    alpha = pm.Deterministic("alpha", tt.exp(alpha_sd * alpha_raw))
    beta = pm.Deterministic("beta", tt.exp(mu / alpha))

    y_obs = pm.Weibull("y_obs", alpha=alpha, beta=beta, observed=y_new[~censored_new])
    #事後対数尤度に制約としてつけ加える。この項を含めた事後対数尤度が最大になるようにMCMCでサンプリングする。
    y_cens = pm.Potential('y_cens', weibull_lccdf(y_new[censored_new], alpha, beta))

with model_3:
    trace_3 = pm.sample(target_accept=0.9, init="adapt_diag")

事後サンプリング結果をプロットする。

_ = pm.traceplot(trace_3)
_ = pm.plot_posterior(trace_3)

サマリーを確認する。

pm.summary(trace_3)

sd   hdi_3%  hdi_97% mcse_mean   mcse_sd ess_bulk    ess_tail    r_hat
mu  8.309   0.755   6.915   9.784   0.031   0.022   591.0   635.0   1.01
a0  0.084   0.009   0.067   0.101   0.000   0.000   585.0   615.0   1.01
alpha   2.319   0.206   1.898   2.693   0.009   0.006   585.0   615.0   1.01
beta    36.039  1.737   32.854  39.437  0.032   0.023   2977.0  2575.0  1.00

2.3.故障時間の比較

以下のようなワイブル分布において、特定の時刻までに生存している確率を算出する関数を定義する。

# cdf
def weibCumDist(x, a, b):
    return  np.exp(-(x / b) ** a)

続いて、故障時間の分布をプロットする。ベイズモデリングではワイブル分布のパラメータを確率的に決定するため、生存時間分布も幅を持った推定が可能になる。

#サンプル数
n_sim = 2000
t = np.linspace(0, 90 , 90)
weibCum_2 = np.zeros((t.shape[0], n_sim))
weibCum_3 = np.zeros((t.shape[0], n_sim))
np.random.seed(42)
rand_val = pd.Series(np.random.randint(1000, 4000, n_sim))

k = 0
for val in rand_val:
    weibCum_2[:, k] = weibCumDist(t, trace_2['alpha'][val], trace_2['beta'][val])
    weibCum_3[:, k] = weibCumDist(t, trace_3['alpha'][val], trace_3['beta'][val])
    k+=1

plt.plot(t, weibCum_2.mean(axis=1))
plt.fill_between(t, np.percentile(weibCum_2, 95, axis=1), np.percentile(weibCum_2, 5, axis=1), alpha=0.5, label='old')

plt.plot(t, weibCum_3.mean(axis=1))
plt.fill_between(t, np.percentile(weibCum_3, 95, axis=1), np.percentile(weibCum_3, 5, axis=1), alpha=0.5, label='new')
plt.legend()
plt.xlabel('Time[Day]')
plt.ylabel('Probability')
plt.show()

2.4.損失額の比較

最後に機会損失額の分布を比較する。

def predict_loss_sample(trace, rand_v, time_to):
    
    alpha_sim = trace['alpha'][rand_v]
    beta_sim = trace['beta'][rand_v]
    t_sim = pd.Series(np.linspace(0, int(time_to), int(time_to) + 1))
    F_t = weibCumDist(t_sim, alpha_sim, beta_sim)

    return F_t.values * 1000 * 100 #生存確率＊生産個数＊１個あたり売り上げ

n_sim = 50
loss_2 = np.zeros(n_sim)
loss_3 = np.zeros(n_sim)
np.random.seed(42)
rand_val = pd.Series(np.random.randint(0, 2000, n_sim))

k = 0 
for val in rand_val:
    loss_2[k] = sum(predict_loss_sample(trace_2, val, 30))
    loss_3[k] = sum(predict_loss_sample(trace_3, val, 30))
    k += 1
    

plt.hist(loss_2, alpha=0.5, label='old')
plt.hist(loss_3, alpha=0.5, label='new')
plt.legend()
plt.xlabel('Opportunity loss amount')
plt.ylabel('Density')
plt.show()

以下、機械損失額の分布の比較。

機会損失額の差の分布を見てみる。

生産装置を入れ替えることで2万円〜7万円の機会損失を抑制できることがわかる。ベイズモデリングによる生存時間解析では,このように確率的に生存時間を推定できるため、機会損失額を分布で予測することが可能である。