大学数学

HOME
大学数学
【生存時間解析】ワイブル分布をわかりやすく｜確率分布関数とパラメータ推定方法

07/27/2024 / 最終更新日時 : 10/30/2024 bluest 大学数学

【生存時間解析】ワイブル分布をわかりやすく｜確率分布関数とパラメータ推定方法

こんにちは、青の統計学です！

今回は、ワイブル分布について解説します。

青の統計学では、noteで統計検定やG検定に関するチートシートを掲載しております。
こちらをクリック！

ワイブル分布

ワイブル分布は、材料強度や寿命データの分析などで広く用いられる連続確率分布です。

製造業でいう故障に関わる解析や、寿命に関するデータのモデリングに使われている印象があります。

指数分布とともに信頼性工学で使われる分布ですね。

確率密度関数は次のように表されます。

$$f(x; k, \lambda) = \begin{cases}
\frac{k}{\lambda}\left(\frac{x}{\lambda}\right)^{k-1}e^{-(x/\lambda)^k} & x \geq 0\\
0 & x < 0 \end{cases}$$

$k > 0$は形状パラメータ（シェイプパラメータともいいます）、$\lambda > 0$はスケールパラメータです。

形状パラメータ$k$の値によって分布の形状が変化します。

$k < 1$のとき分布は狭く左に裾を引く形状、$k = 1$のとき指数分布、$k > 1$のとき狭く右に裾を引く形状となります。

スケールパラメータ$\lambda$は尺度を決めるパラメータで、$\lambda$が大きいほど分布が右に広がります。

ワイブル分布の確率分布関数$F(x)$は次式で表されます。

$$F(x) = 1 – e^{-(x/\lambda)^k}$$

また、$r$次の積率母関数$M_r(t)$は以下のようになります。

$$M_r(t) = \lambda^r\Gamma\left(1 + \frac{r}{k}\right)$$

$\Gamma(x)$はガンマ関数です。

$r = 1$とすると、平均$\mu = \lambda\Gamma(1 + 1/k)$が得られます。

指数分布との関わりについて

さて、形状パラメータが${m=1}$であるとワイブル分布は指数分布になります。

${m=1}$を代入すると${t}$が消滅することに気がつくはずです。

指数分布をわかりやすく解説｜無記憶性

具体的な使い方

先ほども申し上げましたが、

ワイブル分布は、高信頼性製品の寿命分布モデルや、材料強度データの解析などで幅広く活用されています。

実際の応用例を見ていきましょう。

製造業におけるボールベアリングの寿命データを考えてみます。

ある種類のボールベアリングの寿命データ(単位: 時間)は次のようになっているとします。

$$15.8, 16.4, 17.7, 18.1, 19.5, 20.9, 21.4, 23.0, 24.2, 25.3$$

最尤法によってワイブル分布のパラメータを推定してみましょう。

対数尤度関数は以下のようになります。

$$\begin{aligned}
\ell(k, \lambda) &= n\log k – n\log\lambda + (k-1)\sum_{i=1}^n\log x_i – \sum_{i=1}^n\left(\frac{x_i}{\lambda}\right)^k\\
&= 10\log k – 10\log\lambda + (k-1)\sum_{i=1}^{10}\log x_i – \sum_{i=1}^{10}\left(\frac{x_i}{\lambda}\right)^k
\end{aligned}$$

この関数を最大化する$k, \lambda$を数値的に求める必要があります。

さて、こんな結果が得られました、

$$\hat{k} = 3.518, \quad \hat{\lambda} = 20.966$$

つまり、このボールベアリングの寿命データは、シェイプパラメータが$k \approx 3.5$、スケールパラメータが$\lambda \approx 21$のワイブル分布によってモデル化できることがわかります。

一方、別の製品の寿命データに対してワイブル分布をあてはめた結果、$\hat{k} = 0.7$という推定値が得られたとします。

この場合、$k < 1$なので分布形状は左に裾を引いた形になります。

つまり、製品の初期不良が比較的多いことを示唆しています。

$$\begin{align*}
f(x; 0.7, 10) &= \begin{cases}
\frac{0.7}{10}\left(\frac{x}{10}\right)^{-0.3}e^{-(x/10)^{0.7}} & x \geq 0\\
0 & x < 0 \end{cases}\\ F(x) &= 1 – e^{-(x/10)^{0.7}} \end{align*}$$

危険率と生存関数、指数分布との関わりについて

ここで、危険率と生存関数という新しい概念を紹介します。

あるものが故障する確率の分布関数を${F(t)}$とすると、信頼性（つまり故障しない確率みたいなもの）は、以下のように表すことができます。

$${S(t)=1-F(t)}$$

この$${S(t)}$$を生存関数と呼びます。

そして、確率密度関数を生存関数で割ったものを危険率${h(t)}$と呼びます。

ワイブル分布の確率密度関数は

$${f(x;m,\lambda)=\frac{m}{\lambda}(\frac{x}{\lambda})^{m-1}e^{-(\frac{x}{\lambda})^m}}$$

でした。

積分して累積分布関数を取り、1との差分を見たものが生存関数でした。

$${S(t;\lambda,m)=1-F(t;\lambda,m)=1-(1-e^{-(\frac{t}{\lambda})^m})=e^{-(\frac{t}{\lambda})^m}}$$

密度関数を生存率で割ります。

$${h(t;\lambda,m)=\frac{f(t;\lambda,m)}{S(t;\lambda,m)}=\frac{\frac{m}{\lambda}(\frac{x}{\lambda})^{m-1}e^{-(\frac{x}{\lambda})^m}}{e^{-(\frac{t}{\lambda})^m}}=\frac{m}{\lambda}(\frac{t}{\lambda})^{m-1}}$$

この式から、ワイブル分布の危険率は、時間${t}$が増加するとともに変化し、${m}$（形状パラメータ）の値によってその振る舞いが変わることがわかります。

また、時刻 ${t}$までに故障しないもとで，${t}$から ${t+\Delta t}$の間に故障するという条件付き確率は${h(t;\lambda,m)\Delta t}$で表すことができます。