ワイブル分布パラメータの推定手法の比較

前回の記事で、データの規模がある程度大きくても、生データから累積ハザード関数や不信頼度関数を推定してそこにワイブル分布に当てはめる手法(累積ハザード法とメジアンランク法) と生データを直接ワイブル分布に当てはめる方法(最尤法とベイズ推定)で、結果に差があることがわかりました。求めたハイパーパラメータのワイブル分布とEmprical CDFとの比較をすると前者のほうが生データに近そうです。理由を少し深掘りします。

おさらい

どんなデータ？

Engelhardt et al. 1981より孫引きしたアルミA6061-T6の疲労破壊試験データを使いました。

サンプルデータセット
($ 10^3 $ サイクル単位の破壊寿命(サイクル数); N=101)

70, 90, 96, 97, 99, 100, 103, 104, 104, 105,
107, 108, 108, 108, 109, 109, 112, 112, 113, 
114, 114, 114, 116, 119, 120, 120, 120, 121,
121, 123, 124, 124, 124, 124, 124, 128, 128,
129, 129, 130, 130, 130, 131, 131, 131, 131,
131, 132, 132, 132, 133, 134, 134, 134, 134,
134, 136, 136, 137, 138, 138, 138, 139, 139,
141, 141, 142, 142, 142, 142, 142, 142, 144,
144, 145, 146, 148, 148, 149, 151, 151, 152,
155, 156, 157, 157, 157, 157, 158, 159, 162,
163, 163, 164, 166, 166, 168, 170, 174, 196, 212

懸案のCDF比較結果

累積ハザード法とメジアンランク法がEmpirical CDFによく当てはまっている一方で、最尤法とベイズ推定は分布の端っこで乖離している傾向でした。

仮説と検証

ECDFへの当てはまりのよい手法はチート(Match fixing)では？

累積ハザード法とメジアンランク法がフィットしているのは、各手法で求めた累積ハザード関数や不信頼度関数がEmpirical CDFとよく一致しているからでは？ ⇔ 実質ECDFに対する最小二乗法をした状態になっているからでは？

累積ハザード法で求めた累積ハザード関数 $ \hat{H}_{cum} $, メジアンランク法で求めた不信頼度 $ \hat{F}_{med} $ およびEmprical CDF $ \hat{F}_{emp} $ をプロットします。$ F(t) = 1 - exp(H(t)) $ の関係を使います。

plt.plot(x_emp, y_emp, marker='.', linestyle='none', label='ecdf')
plt.plot(df['k-cycles'], 1-np.exp(-df['H']), marker='+', linestyle='none', label='ecdf')
plt.plot(df['k-cycles'], df['med-rank-F'], marker='^', linestyle='none', label='ecdf')
plt.xlabel('k-cycles')
plt.ylabel('CDF')
plt.legend(['Empirical CDF','Cum-Harzard CDF','Med-Rank CDF'])

なんかもう… 言うまでもなく仮説どおりっぽい。

この手法で最大化している尤度関数は\[ \mathcal{L}(m, \eta, \sigma^2; ln ln \frac{1}{1-\hat{F}(X)},ln X) = (2\pi\sigma^2)^{-N/2}exp(-\frac{1}{2\sigma^2}\sum_{i=1}^N(ln ln \frac{1}{1-\hat{F}(x_i)}-(m ln x_i+m ln \eta))^2) \] であって、これがEmpirical CDFへの回帰問題($ \mathcal{L}(m, \eta, \sigma^2; \hat{F}(X), X) = (2\pi\sigma^2)^{-N/2}exp(-\frac{1}{2\sigma^2}\sum_{i=1}^N(\hat{F}(x_i)-(1-exp(-(x_i/\eta)^m)))^2) $ の最大化) に近いということだったわけですね。

最尤法では何が起きたか？

最尤法で直接分布を推定する場合、尤度関数は\[ \mathcal{L}(m, \eta; X) = \prod_{i=1}^N \frac{m}{\eta}(\frac{x_i}{\eta})^{m-1}exp(-(\frac{x_i}{\eta})^m) \] です。Log-Likelihood関数は \[ \ell(m, \eta; X) = \sum_{i=1}^N ln(\frac{m}{\eta})+(m-1)ln(\frac{x_i}{\eta})-(\frac{x_i}{\eta})^m \] から、ちょっと絵を描いてみます。

特性寿命に比べてワイブル係数のほうが穏やかな勾配になっている印象。ワイブル係数の"誤差"は尤度に与える影響としては低そうなので、特性寿命を143に固定してワイブル係数が6の場合(≒最尤法やベイズで求めた値)と7の場合(≒メジアンランク法等で求めた値)のLog-Likelihood関数の値を眺めます。

def LLF_Weibull(x, m, eta):
  return -(x/eta)**m+(m-1)*np.log(x)+np.log(m)-m*np.log(eta)

df['MLE(143,6)']=LLF_Weibull(df['k-cycles'],6,143)
df['MLE(143,7)']=LLF_Weibull(df['k-cycles'],7,143)
ax = df[["MLE(143,6)","MLE(143,7)"]][::4].plot.bar(figsize=(13,4))
ax.xaxis.tick_top()

尤度関数に与える影響は裾の端っこの方が大きく、特に212k-cyclesや196k-cyclesといったデータの影響でワイブル係数が7から6に引っ張られていると見ても良さそうです。

最後にサンプルと確率密度関数を比べてみます。

def weibull_pdf(x, eta, m):
  return m/eta*(x/eta)**(m-1)*np.exp(-(x/eta)**m)

num_bins=20

fig, ax = plt.subplots()
n, bins, patches = ax.hist(df['k-cycles'],num_bins, density=1)
y1 = weibull_pdf(bins, 143, 6)
y2 = weibull_pdf(bins, 143, 7)
ax.plot(bins, y1, '--')
ax.plot(bins, y2, '--')
plt.legend(['6','7','df'])
plt.show()

分布の裾野「起こる確率の低いはずのことが起こったのだから、その分、『裾が高い』方が正しいのだ」とでも言われているような気がしてきます。

お粗末！！