津田尚胤, 貝戸清之, 山本浩司, 小林潔司 (2006) ワイブル劣化ハザードモデルのベイズ推計法. 土木学会論文集Ｆ, 62(3), 473-491.
　土木施設の劣化を統計的に予測するためにハザードモデルをベイズ推定します、という話。ひょんな事情がありまして、昼飯のお茶漬けを啜りながら目を通した。えーっと、筆頭著者の方の修士論文とかかしらん？

　故障解析なんてもう全くの門外漢なので、ごく気楽に読み始めたんだけど、数ページ目で、これって要は生存モデルだよな、と気が付いた。しかし私は物事を簡単に忘れてしまう。というわけで途中で躓いてしまい、ほうじ茶をいれて生存時間解析の本をめくり復習する羽目になった。なにをやっているんだか。

　論文のメモに入る前に、あらかじめ復習しておきます。
　生存時間\(\zeta\)を確率変数ととみて、時点\(t\)における生存確率を\(S(t) = P(\zeta \geq t) \)とします。\(\zeta\)の累積分布関数を\(F(t)\)、密度関数を\(f(t)\)とします。
　\(F(t) = 1-S(t)\)である(死ぬ年齢\(\zeta\)が\(t\)より左にある確率とは、つまり\(t\)の時点で死んでる確率だから)。\(f(t) = \lim_{\Delta t \rightarrow 0} \frac{F(t + \Delta t) – F(t)}{\Delta t} \)である。ここまではいいっすね。
　\(t\)まで生きている人がその瞬間にばたりと死ぬ条件付き確率、つまりハザードを\(h(t)\)とする。$$ h(t) = \lim_{\Delta t \rightarrow 0} \frac{P(t \leq \zeta < t + \Delta t | \zeta \geq t) }{\Delta t} $$ $$ = \lim \frac{S(t) – S(t+\Delta t)}{S(t) \Delta t} = \lim \frac{F(t + \Delta t) – F(t)}{S(t)\Delta t} = \frac{f(t)}{S(t)}$$ である。
　\(S(t)=1-F(t)\)の微分は当然ながら\(-f(t)\)になる。ってことは、\(\log S(t)\)を微分しますと、関数の対数の微分だから $$ \frac{d \log S(t)}{dt} = \frac{1}{S(t)} \times -f(t) = -h(t)$$ となる。ってことは、累積ハザード \( H(t) = \int_0^t h(u)du \) というのを考えれば\(\log S(t) = -H(t)\)という関係があるわけよ。
　なんだか、あれだな… 俺これ年に一回くらい復習してる気がするな…

　ここから論文のメモ。
—–
　土木施設の点検で得られる情報を次の2つにわけて考える。

• 完全モニタリング情報。施設の使用開始時点と故障時点がわかるような情報。カレンダー時点を\(\tau\)と書くとして、\(\tau_1\)から使い始めて\(\tau_a\)に故障しました、修理して\(\tau_2\)から使用再開したら\(\tau_b\)に故障しました…というような情報である。完全モニタリング情報を得るには故障を常時監視しないといけないけど、これは土木施設では無理がある。
• 定期モニタリング情報。たとえば\(\tau_W\)と\(\tau_T\)に一斉点検しました、というような情報。

　以下、使用開始-故障発生の区間を完全サンプル、現在使用中でまだ故障してない区間を不完全サンプルと呼ぶ。[おおお… 施設をサンプルと呼ぶんじゃなくて、イベントまでの期間をサンプルと呼ぶのね。へえええ。こういうことがあるから、全然知らない分野の話を読むのは面白い]
　その施設の直近の使用開始からの経過時間を\(t\)とする。
　定期モニタリング情報の場合、施設\(i\)の直近使用開始から点検までの時間(つまり\(t\)上での点検時点)を\(W_i, T_i\)と書くとして、もし「\(t=W_A\)では故障なし、\(t=T_A\)では故障あり」だったら、少なくとも\( (W_A, T_A] \)のどこかで故障が発生したという完全サンプルが手に入ったことになり、「どちらにおいても故障なし」だったら不完全サンプルが手に入ったことになる。
　
　ワイブル劣化ハザードモデルについてのご紹介。[この分野では常識に属することなのだろうが、私にとっては目新しい話題だ… ちょっと楽しい]
　土木施設を\(i (=1, \ldots, n)\)とする[たぶん、正確には土木施設そのものじゃなくて”サンプル”のことであろう]。それぞれの\(i\)について\(t=0\)から\(t=T_i\)までの完全モニタリング情報が手に入るとする。使用時間長を\(y_i\)、寿命を\(\zeta_i\)とする。故障によってモニタリングが終了した場合は\(y_i = T_i = \zeta_i\)であり、故障前にモニタリングが打ち切られた場合は\(y_i = T_i < \zeta_i\)である。[はいはい。\(y_i\)は観察打ち切りつきの生存時間のことね]
　\(\zeta_i\)を確率変数とみて、その密度関数を\(f_i(\zeta_i)\)、分布関数を\(F_i(\zeta_i)\)とする。
　任意の時点\(t\)における生存確率を\(S_i(t) = 1 – F_i(t) \) と書く。[論文では任意の時点\(t\)も\(\zeta_i\)という記号を使い、\(S\)を\(\tilde{F}\)と書いているんだけど、私には読みにくいので勝手に書き換えた。すいません]

　劣化ハザード関数、つまり時点\(t\)まで生存しそこで初めて故障する確率密度を\(h(t)\)とする[論文では\(h(t)\)じゃなくて\(\lambda(\zeta_i)\)。すいません]。$$ S_i (t) = \exp \left(-\int_0^{t} h_i(u) du \right)$$ ってことになりますね。
　ワイブル劣化ハザード関数を用いる。劣化の加速度を表すパラメータを\(m\)として $$ h_i(t) = \gamma_i m t^{m-1}$$ ただし、施設の特性ベクトルを\(x_i\)として $$ \gamma_i = \exp( x_i \beta^\top) $$ とする。
　[へー。ワイブル分布のパラメータのうちワイブル係数\(m\)のほうはそのままで、尺度パラメータが共変量で変わると考えるわけね。
　なお、原文では\(t\)の関数じゃなくて\(y_i\)の関数として書いている。この論文は「検査した施設のなかにいま\(y_i\)年目の施設があって…」という発想なのでそう書いているんだと思うけど、私からみると「観察された生存時間が\(y_i\)の人の瞬間死亡率？ その人のどの時点での瞬間死亡率のことよ？？？」って感じになるので、失礼ながら\(t\)の関数に勝手に書き換えた。度重なる無礼をお許しください]
　密度関数は $$ f_i(t) = \gamma_i m y^{m-1} \exp (-\gamma_i t^m) $$ 生存関数は $$ S_i(t) = \exp (-\gamma_i t^m) $$ となる。[好き勝手に書き換えてたら原文との乖離がだんだん激しくなってきちゃった… 原文の左辺は\(f(y_i), \tilde{F}(y_i) \)である]

　どうやって推定するか。以下、実測値には上にバーをつけることにする。未知パラメータをまとめて\(\theta = (m, \beta)\)と書く。
　まずは完全モニタリング情報が手に入っている場合。モニタリング期間中に施設寿命が終了したときに1になるダミー変数[いわば「死亡観察フラグ」ね、観察打ち切りフラグじゃなくて] を\( d_i \)とする。完全モニタリングによって得られる情報をまとめて\(\bar{\xi}_i = (\bar{d}_i, \bar{y}_i, \bar{x}_i), \bar{\xi} = (\bar{\xi}_1, \ldots, \bar{\xi}_n) \) と書く。
　密度関数と生存関数を\( f(\bar{y}_i, \bar{x}_i:\theta), S(\bar{y}_i, \bar{x}_i : \theta) \) と書こう。施設の劣化の生起が互いに独立だとすれば、尤度関数は$$ L(\theta | \bar{\xi}) = \prod_i^n f(\bar{y}_i, \bar{x}_i:\theta)^{\bar{d}_i} S(\bar{y}_i, \bar{x}_i : \theta)^{1-\bar{d}} $$ である。ワイブル劣化ハザードモデルならば, モニタリング中に寿命を迎えた施設数を\( \bar{d} = \sum_i^n \bar{d}_i \)として $$ L(\theta | \bar{\xi} ) = m^\bar{d} \exp \left( \sum_i^n (\bar{d}_i \bar{x}_i \beta^\top + \bar{d}_i (m-1) \log \bar{y}_i – \exp (\bar{x} \beta^\top) \bar{y}_i^m ) \right) $$ である。これを最大にする\(\hat{\theta} \) を求めればよい。
　定期モニタリング情報の場合。\(\hat{S}(W_i, T_i, \tilde{x}_i : \theta) = S(W_i, \tilde{x}_i : \theta) – S(T_i, \tilde{x}_i : \theta) \)と略記する[原文は\(\hat{S}\)ではなく\(\hat{F}\)]。期間\( (W_i, T_i] \)に故障したことを表すダミー変数を \(d^f_i\)として、尤度関数は$$ L^f(\theta | \xi) = \prod_i^n \hat{S}(\bar{W}_i, \bar{T}_i, \bar{x}_i : \theta)^{\bar{d}^f_i} \times \left( 1-S(\bar{W}_i, \bar{T}_i, \bar{x}_i : \theta) \right)^{1-\bar{d}^f_i} $$ となる。

　ようやく本題です。
　データが足りなくて最尤推定がよくなかったり、データが全然なかったりする場合であってもモデルを特定したい。そこでモデルをベイズ更新していくことを考える。
　残念ながら、\(\theta\)が未知であるとき、うまい共役分布はない[Ibrahim, et al.(2002 “Bayesian Survival Analysis”)をreferしている。買ったきり全然開いてない本だ…]。しょうがないのでMCMCでなんとかすることにします。
[ここでGibbsサンプリングの説明。パス]

　事前分布として、\(m\)にはガンマ分布\(G(m_0, \kappa_0)\)を、\(\beta\)には\(N_k(\mu_0, \Sigma_0)\)を与えることにします。
　ここで先験的知識を組み込むことができる。次の4つの場合について考える。

• Case a. すでに利用可能な既存のワイブル劣化ハザードモデルが存在する場合。このときは、\(m\)の期待値と分散、\(b\)の期待値と分散共分散行列が手に入っている。
• Case b. 参照となる類似のワイブル劣化ハザードモデルが存在する場合。\(m, \beta\)の期待値についてはわかっているけど分散・共分散はわからない。
• Case c. 利用可能なモデルはないが技術者が主観的に判断できる場合。\(m\)の期待値についてはわかるし、\(\beta^1\)の期待値についてもわかるがあとはわからない。[←\(\beta^1\)の期待値ってなに？ あっそうか、切片の期待値か]
• Case d. なんもわからん場合。\(m \sim G(1, \kappa_m^{-1}), \beta \sim N_k (0, \kappa_\beta I) \)とするか( \(\kappa_m, \kappa_\beta\)は大きな値)、Jeffreys分布を使うしかない。

[MCMCで得たサンプルをどう要約するかという説明。パス]
[モデル選択の話。Chibの方法で周辺尤度を求めてベイズファクターを使う。パス]
[収束判断の話。パス]

　さて、1回目の点検データ\(\bar{\xi}^0\)が来たので事後分布を求めたとするじゃあないですか。そこに2回目の点検データ\(\bar{xi}_1\)が来たら、\(m, \beta\)の事前分布を\(f, g\)として、パラメータの事後密度は$$ \pi(\theta | \bar{\xi}^0, \bar{\xi}^1) \propto L(\theta|\bar{\xi}^1) \pi(\theta|\bar{\xi}^0) \propto L(\theta | \bar{\xi}^0, \bar{\xi}^1) f(m | m_0, \kappa_0) g(\beta | \mu_0, \Sigma_0)$$じゃないですか。こういう風に、新しいデータを追加したデータベースについて尤度関数を定義しなおせばいいわけよ。[はいはい。1回目の点検データで得たパラメータの事後分布に2回目の追加データの尤度を掛けて更新するというより、毎度毎度、積み上げデータの尤度に事前分布を掛ければええやないの、という話ね]

　適用事例。
　東北自動車道十和田道路事務所では、12年間にわたって、トンネル照明の使用開始日と、それが点かなくなった日を記録している。全部で7145個。うち完全サンプルは4348個。[なぜ完全モニタリング情報が手に入っているのだろうか… 毎日調べて回っているのかなあ…]
　ワイブル劣化ハザード関数を $$ h_i(y_i) = \exp( x_i \beta^\top) m y_i^{m-1} $$ とする。\(x_i\)は、｛切片項、低圧ランプか高圧ランプか、照明種類ごとの日平均点灯時間｝。
　まずはデータを全部使って最尤推定しておく。で、データベースから200個を非復元抽出したデータ、さらに200個を非復元抽出して追加して400個にしたデータ, …, 1000個にしたデータ、を作っておく。それぞれについてGibbサンプリングでパラメータを推定する。[論文ではこれをベイズ更新と呼んでいるので、つい、最初の200個データを使ってサンプリングした事後分布をつかってなにかをするのかと思ったが、そうではなさそう]。
　期待通り、だんだん良い推定値になりました。

　Case c. の場合について考えよう。
　いま、現場の技術者が「経過年数\(y_1\)でおよそ\(1-p_1\)のランプがつかなくなる。\(y_2\)だとおよそ\(1-p_2\)だ」ということを知っていたとしよう。このとき、$$ p_1 = \exp(-\gamma y_1^m) $$ $$ p_2 = \exp(-\gamma y_2^m) $$ より $$ m^o = \frac{\log(-\log p_1) – \log(-\log p_2)}{\log y_1 – \log y_2}$$ $$ \gamma^0 = – \frac{\log p_1}{y_1^{m^o}} $$ と見積もれる。
　\(m\)の事前分布は\(G(m_0, 1)\)としよう。\(\beta^1\)の事前分布を\(N(\log \gamma, 1)\)としよう[\(\gamma\)って？ 上の\(\gamma^o\)のこと？]。\(\beta^2, \beta^3\)にはJeffreys分布を与えよう。で、試してみると[…中略…] 現場の技術者の経験情報が正しければ、データベースが小さいとき、パラメータの推計精度をCase d.よりも改善できる。

　最後に、定期モニタリングスキームのデータを作って試してみると…[パス]
　云々。
—–
　勉強になりましたですー。　世の中には、いろんな研究をしている人がいるんだなあ…