読了:Hallinan (1993) ワイブル分布七変化 (正確には五変化)

Hallinan, A.J. (1993) A Review of the Weibull Distribution. Journal of Quality Technology. 25(2), 85-93.

 先日、ひょんな事情により目を通した論文。
 ワイブル分布について調べていて、資料によって式が違うので大変混乱し、解説を探して見つけた次第である。ワイブル分布って工学とかで使う奴でしょ? 私は文系も文系なのにそんなの知るわけないじゃん? と思うわけだが、生きていく上ではいろんなことが起きる。

 いわく。
 Weibull分布の定式化が本によって異なるので実務家は混乱する。なんと5つも定式化があるのだ。歴史的背景とともに整理いたします。

 Waloddi Weibullは1887年生まれ。彼の最初の論文は1914年[←第一次大戦がはじまった年だ]、彼はスウェーデンの沿岸警備隊の中尉であった。艦長を務めて1917年に退役、1924年に学位取得、スウェーデン王立工科大学の正教授となる。
 [ここからこの論文は、Weibullが折々の論文に書いたオリジナルの数式をまず紹介し、次にそれを統一的な記法で書き換えて示している。以下では後者の、統一的記法による表現のみをメモする]

定式化1
 1939年、彼はストレス下の材料の強度の分布についての論文を書く。ストレスを\(x\)として、material functionを$$ \phi(x) = \left( \frac{x-c}{a} \right)^b $$ と定義した。
 [material functionというのがわからなくて困惑したんだけど… 生存分析で言うと、たとえばハザードが定数\(\lambda\)だったら、生存時間関数は\(S(t) =\exp(-\lambda t)\)じゃないですか。この\(\lambda t\)の部分をmaterial functionと呼んでいるらしい。実質的にどういう概念だと思えばいいのかよくわからんが、この分野ではわかりやすい概念なんでしょうね]
 \(a, b, c\)は定数で、スケール、形状、位置と呼ばれている。スケールは「故障の間の平均時間」、位置は「故障確率が非ゼロとなる最小の時間」に相当する。

 時間\(t\)までに故障する確率は$$ F(t) = 1 – \exp \left( -\left(\frac{t-c}{a} \right)^b \right)$$ となる。なお、時間\(t\)までに故障しない確率を信頼性\( R(t) = 1-F(t) \)と呼ぶ。[生存時間関数のことね]
 \(F(t)\)は累積分布関数だが、確率密度関数は $$ f(t) = ba^{-b} (t-c)^{b-1} \exp \left( -\left(\frac{t-c}{a} \right)^b \right)$$ となる。ハザード \( h(t) = \frac{f(t)}{1-F(t)} \)は $$ h(t) = ba^{-b} (t-c)^{b-1}$$ となる。

定式化2
 第二次大戦中、彼は材料の強度の研究を周期的疲労の研究へと拡張し、次のように定式化する。$$ \phi(x) = k(x-c)^b$$ ここで\(x\)は負荷、\(u\)は故障確率が0となる負荷を指していた。Weibullは\(k\)に名前をつけなかったが、スケールと形状が一緒になったパラメータである。

 この定式化だと、$$ F(t) = 1-\exp \left( -k (t-c)^b \right)$$ $$ f(t) = bk(t-c)^{b-1} \exp \left(-k(t-c)^b \right)$$ $$h(t) = bk(t-c)^{b-1}$$ となる。

定式化3
 さて、1951年に彼は自分がつくった分布を米国に紹介するんだけど、この論文に誤植があり、定式化1のつもりで書いたのに、かっこが分子だけにかかる形になってしまった。\(k’ = 1/k\)として$$ \phi(x) = \frac{(x-c)^b}{k’} $$ と、意図せず定式化2の形で書いちゃったことになる。

 この定式化だと、$$ F(t) = 1-\exp \left( -\frac{(t-c)^b}{k’} \right)$$ $$ f(t) = b(k’)^{-1} (t-c)^{b-1} \exp \left( – \frac{(t-c)^b}{k’} \right)$$ $$h(t) = b(k’)^{-1} (t-c)^{b-1}$$ となる。

定式化4
 のちの研究者らは、計算しやすくするために、Weibullの定式化1を次のように書き換えることがある。\(a’ = 1/a\)として$$ \phi(x) = \left( a’ (x-c) \right)^b $$ \(a’\)は「単位時間あたり故障回数」に相当する。ややこしいことに、この\(a’\)のことをスケールと呼ぶことがある。
 この定式化だと $$ F(t) = 1 – \exp \left( – \left( a’ (t-c) \right)^b \right)$$ $$ f(t) = b(a’)^{b} (t-c)^{b-1} \exp \left( – \left( a’ (t-c) \right)^b \right)$$ $$ h(t) = b(a’)^{b} (t-c)^{b-1}$$ となる。

定式化5
 1953年、Weibullは王立工科大学を退職するが、仕事は続ける。1961年に本を出版するんだけど、そこでの定式化がこれまでと違っているもんだからややこしい。整理すると、\(b’ = 1/b\)として $$ \phi(x) = \left( \frac{x-c}{a} \right)^{1/b’} $$ と書いたことになる。なお、この本でWeibullは\(a\)を\(\beta\)と書いたので、\(a\)はbetaパラメータと呼ばれることがある。
 この定式化だと $$ F(t) = 1 – \exp \left( – \left( -\frac{t-c}{a} \right)^{1/b’} \right)$$ $$ f(t) = (1/b’)a^{-1/b’} (t-c)^{(1/b’)-1} \exp \left( – \left( \frac{t-c}{a} \right)^{1/b’} \right)$$ $$ h(t) = (1/b’)a^{-1/b’} (t-c)^{(1/b’)-1}$$ となる。 
 Weibullは数々の名誉ある賞を受け、1979年にフランスで死去した。

[論文の後半は、ワイブル用紙を使った推定、数値例、先行するFisherの論文との関連について。読まずにとばした]
—–
 なるほど… 私が混乱した理由がわかった。

  • たまたま机の横にあった中村「Cox比例ハザードモデル」(朝倉書店)、西川「カプラン・マイヤー法」(共立出版)をめくると、どちらも、Weibull分布に従う生存時間関数を$$S(t) = \exp( -(\lambda t)^p )$$と書いている。生存分析ではこの形で書くことが多いのではなかろうか。私もこの式をイメージしてました。
     この論文の表記で言うと定式化4だ。この論文でいう\(a’\)を\(\lambda\)と書き、この論文でいう形状\(b\)を\(p\)と書き、この論文でいう位置\(c\)をゼロにしているわけだ。よく見ると、西川本には\(\lambda\)のことをスケールと呼んだり\(1/\lambda\)のことをスケールと呼んだりすると書いてある。なるほど、こういう経緯だったのか。
  • webで検索したところ、日本機械学会による解説では、Weibull分布関数を$$ F(x) = 1 – \exp \left( – \left( \frac{x – \gamma}{\beta} \right)^\alpha \right) $$ と表現している。定式化1ですね。
  • RのWeibull分布の関数のマニュアルをみると、累積分布関数は$$ F(x) = 1 – \exp \left( – \left( \frac{x}{b} \right)^a \right) $$ とし、\(a\)を形状, \(b\)をスケールと呼んでいる。そうか、これも定式化1だ。この論文でいうスケール\(a\)を\(b\)と書き、この論文でいう形状\(b\)を\(a\)と書いているのだ。
  • Wikipediaのワイブル分布の説明では、確率分布を $$ f(t) = \frac{m}{\eta} \left( \frac{t}{\eta} \right)^{m-1} \exp \left( – \left( \frac{t}{\eta} \right)^m \right)$$ とし、\(m\)をワイブル係数、\(\eta\)を尺度パラメータと呼んでいる。これも定式化1だ。右辺前半のハザードの部分が見慣れなくて目が回っちゃうけど、要するに、この論文でいうスケール\(a\)を\(\eta\)と書き、この論文でいう形状\(b\)を\(m\)と書き、この論文でいう位置\(c\)をゼロにしているのである。

 こうしてみると、Weibull分布の式を理解するコツは、分布関数の\(\exp\)の内側をみることだな。\(t\)を一次変換してほにゃらら乗して負をとっているやつ(定式化1,4,5)と、\(t\)から定数を引いてほにゃらら乗して定数を掛けてから負をとっているやつ(定式化2,3)があるわけだ。ハザード関数や確率密度関数を真面目に理解しようとするとかえって混乱する(ハザードの書き方はいろいろあるから)。また、「スケール」とか「形状」といったパラメータ名称にこだわるのも混乱のもとになる可能性がある(\(a\)をスケールを呼ぶこともあれば\(1/a\)をスケールと呼ぶこともあるから)。