ワイブル分布の最尤推定｜閉形式と数値計算をわかりやすく解説

ワイブル分布は、寿命や故障時間を表すときによく登場する連続分布です。最尤推定では、形状母数を固定すればスケール母数をきれいな式で求められます。しかし、2つを同時に推定すると、最後は数値計算が必要です。

この記事では、その「途中までは式で解けるけれど、最後は手計算だけでは終わらない」流れを順番に見ていきます 🧮

尤度・対数尤度・スコア関数の基本から確認したい方は、先に「最尤法とは？尤度・対数尤度・スコア関数と推定の手順」をご覧ください。

この記事で使うワイブル分布

本記事では、スケール母数を \(\theta>0\)、形状母数を \(\tau>0\) とし、生存関数を、

\begin{align*} S(x) =P(X>x) =\exp\left[-\left(\frac{x}{\theta}\right)^\tau\right], \qquad x>0 \end{align*}

と置きます。文献によってパラメータの置き方が違うため、式を比べるときはこの定義を先に確認してください。

確率密度関数は \(f(x)=-S'(x)\) より、

\begin{align*} f(x\mid\theta,\tau) =\frac{\tau}{\theta} \left(\frac{x}{\theta}\right)^{\tau-1} \exp\left[-\left(\frac{x}{\theta}\right)^\tau\right] \end{align*}

となります。

まずは形状母数 \(\tau\) が既知の場合

\(X_1,\ldots,X_n\) を独立同分布な標本とし、観測値 \(x_i\) はすべて正、打ち切りのない完全データとします。まず \(\tau\) は既知で固定し、未知の \(\theta\) だけを推定します。

独立性から尤度関数は、

\begin{align*} L(\theta;\mathbf{x},\tau) &=\prod_{i=1}^{n}f(x_i\mid\theta,\tau)\\ &=\tau^n\theta^{-n\tau} \left(\prod_{i=1}^{n}x_i^{\tau-1}\right)\\ &\quad\times \exp\left[ -\sum_{i=1}^{n}\left(\frac{x_i}{\theta}\right)^\tau \right] \end{align*}

です。対数尤度は、

\begin{align*} \ell(\theta,\tau) &=n\log\tau +(\tau-1)\sum_{i=1}^{n}\log x_i\\ &\quad -n\tau\log\theta -\sum_{i=1}^{n}\left(\frac{x_i}{\theta}\right)^\tau \end{align*}

となります。

スケール母数 \(\theta\) の最尤推定量

\(\tau\) を定数として、対数尤度を \(\theta\) で微分します。

\begin{align*} \frac{\partial\ell}{\partial\theta} =-\frac{n\tau}{\theta} +\frac{\tau}{\theta^{\tau+1}} \sum_{i=1}^{n}x_i^\tau \end{align*}

これを0と置くと、

\begin{align*} \hat\theta^\tau &=\frac{1}{n}\sum_{i=1}^{n}x_i^\tau\\ \hat\theta(\tau) &=\left( \frac{1}{n}\sum_{i=1}^{n}x_i^\tau \right)^{1/\tau} \end{align*}

を得ます。停留点での2階微分は、

\begin{align*} \left. \frac{\partial^2\ell}{\partial\theta^2} \right|_{\theta=\hat\theta} =-\frac{n\tau^2}{\hat\theta^2} <0 \end{align*}

なので、この点は最大値です。

形状母数 \(\tau\) も未知なら数値計算へ

実際には \(\theta\) と \(\tau\) の両方が未知であることも多いです。先ほどの \(\hat\theta(\tau)\) を対数尤度へ代入して整理すると、\(\tau\) は、

\begin{align*} \frac{1}{\tau} +\frac{1}{n}\sum_{i=1}^{n}\log x_i -\frac{\sum_{i=1}^{n}x_i^\tau\log x_i} {\sum_{i=1}^{n}x_i^\tau} =0 \end{align*}

という方程式を満たします。ここでは \(\tau\) がべき乗部分と重みの両方に入っているため、一般には \(\hat\tau=\cdots\) というきれいな閉形式にはできません。

標本に少なくとも2つの異なる値があれば、この方程式は一意な正の解を持ちます。すべての観測値が同じ値 \(c\) なら、左辺は常に \(1/\tau>0\) となり、\(\tau\to\infty\) でプロファイル尤度が上限なく大きくなるため、有限の最尤推定値は存在しません。

そこで、Newton法などの根を探す方法、または対数尤度を直接最大化する方法で \(\hat\tau\) を数値的に求めます。最後に、その値を、

\begin{align*} \hat\theta =\left( \frac{1}{n}\sum_{i=1}^{n}x_i^{\hat\tau} \right)^{1/\hat\tau} \end{align*}

へ代入すれば、\(\hat\theta\) も決まります。

実データでの注意点 ⚠️

• パラメータ化を確認する：スケール母数や率母数の置き方が違うと、式の見た目も変わります。
• 正の範囲で計算する：\(\theta>0,\tau>0\) という制約を数値計算へ入れます。
• 初期値を変えて確かめる：数値最適化では、収束したか、最大値かを確認します。
• 打ち切りデータは尤度が異なる：生存時間分析で右打ち切りがある場合、本記事の完全データ用尤度をそのまま使うことはできません。

まとめ

形状母数 \(\tau\) が既知なら、ワイブル分布のスケール母数は、

\begin{align*} \hat\theta =\left( \frac{1}{n}\sum_{i=1}^{n}x_i^\tau \right)^{1/\tau} \end{align*}

と求められます。一方、\(\tau\) も未知なら、形状母数の方程式を数値的に解く必要があります。閉形式で解ける幾何分布との違いを比べると、最尤法における数値計算の役割が見えやすくなります。

別の推定方法も比べたい方は、次の「モーメント法とは？母平均＝標本平均の発想と最尤法との使い分け」へ進んでください。