数理解析研究所講究録
第 1857 巻 2013 年 58-68

58

小修理を伴うブロック取替え問題に対するブートストラップ法の適用について
広島大学大学院工学研究科 齋藤靖洋 (Yasuhiro Saito) 土肥正 (Tadashi Dohi)
Graduate School of Engineering, Hiroshima University

1

はじめに
ブロック取替え問題は予防保全計画問題の一種であり，年齢取替え問題と並ぶ重要かつ代表的な問題として，Barlow and Proschan によって考察された．様々な文献の中で，このモデルに対する拡張研究が行われてきた [2]. 本稿では，ブロック取替え問題から派生した問題である小修理を伴うブロック取替え問題$[1|$に対し，パラメトリックブートストラップ法を適用して最適予防保全間隔の推定量の確率分布を導出し，高次モーメントの算出や区間推定を行うことを目的とする．予防保全計画問題にブートストラップ法を適用した従来研究として，Tokumoto et al. [3] は故障発生時刻がワイブル分布に従うと仮定した年齢取替え問題に対し，ブートストラップ法を適用した．彼らは，定常状態における単位時間当りの総期待費用を最小にするような最適年齢取替え間隔の推定量の確率分布を導出することを目的として研究を行っていた．これに対して本稿では小修理を伴うブロック取替え問題を対象として，ブートストラツプ法の適用に関する研究を行う．具体的に，パラメトリックブートストラップ法を適用し，シミュレーションを通して本研究で行った 7 ートストラップ法の有効性を検証するとともに，本手法を実データに適用した数値例についても報告$\grave{}$する。

2

小修理を伴うブロック取替え問題

2.1

モデル化

1 つのコンポーネントに対する小修理を伴うブロック取替え問題を考える．コンポーネントの故障時刻を表す非負の確率変数をとおく． は確率分布関数 $Pr(T\leq t)=F(t)$ 及び確率密度関数 $dF(t)/dt=f(t)$
$T$

$T$

$F(t)$ はを持つものとし，に間して連続な単調増加関数とする．更に， を予防保全時間間隔とおく．小修理を伴うブロック取替え問題では，時刻 $t_{0}(>0)$ 以前に故障したコンポーネントは故障直前の状態に応急$t$$t_{0}$修理する．この時，コンポーネントは故障直前の状態に戻されることから，修理後は以前の故障率を維持した状態で動作を再開するものと考える．このような故障率が変化しない応急修理のことを小修理と呼ぶ．
また，時刻 が経過した時点で，使用されたコンポーネントは新品に取り替えることとする．こうすることにより，年齢取替えのように過去の取替えの履歴を把握する必要が無く，管理者にとって容易な予防保全の方法であると言える．一方，応急修理したコンポーネントであっても，時刻 の時点で新品に取り替えることとなるため，余分な費用がかかるというデメリットも存在する．本稿では，小修理にかかる費用を$t_{0}$$t_{0}$とし， ごとの周期的な予防保全にかかる費用をで表す．更に，時刻 までに観測される期待累積故障数をとおけば，一周期にかかる総期待費用はとなる．したがって，単位$c_{1}(>0)$$c_{2}(>c_{1})$$t_{0}$$t_{0}$$c_{1}\Lambda(t_{0})+c_{2}$$\Lambda(t_{0})$時間当たりの総期待費用 $C(t_{0})$ は$C(t_{0})= \frac{c_{1}\Lambda(t_{0})+c_{2}}{t_{0}}$(1)59となる．ここでの問題は式 (1) を最小にするような最適予防保全間隔$t_{0}^{*}$を求めることとなる．関数$C(t_{0})$に関する一階の最適性条件より $q(t_{0}^{*})=0$ を得る．ここで，(2)$q(t_{0})=c_{1}\{t_{0}\lambda(t_{0})-\Lambda(t_{0})\}-c_{2}$であり，$\lambda(t_{0})=d\Lambda(t_{0}\rangle/dt_{0}$$(d\lambda(t)/dt>0)$はコンポーネントの故障強度関数を表す．時間の経過とともに強度関数が増加するという仮定の下で，$q(\infty)>0$ ならば最適予防保全間隔$t_{0}^{*}(0<t_{0}^{*}<\infty)$

は唯一つ存在する．

2.2

非同次ボアソン過程

本稿では，システムの故障発生事象の振る舞いを非同次ボアソン過程 (以下，NHPP) を用いてモデル化する．NHPP では，時刻

.
.
.

$t$

までに検出される累積フォールト数 $N(t)$ は以下の性質をもつ計数過程にょって記述されるものとする．

.

$N(0\}=0$$N(t)$ は独立増分を持つ$Pr\{N(t+\Delta t)-N(t)\geq 2\}=o(\Delta t)$$Pr\{N(t+\Delta t)-N(t)=1\}=\lambda\triangle t+o(\Delta t)$

この時，確率過程 $N(t)$ は以下のように表すことができる．
$Pr\{N(t)=n\} = \frac{\{\Lambda(t)\}^{n}}{n!}\exp\{-\Lambda(t)\}, n=0,1,2, \cdots$$\Lambda(t) = \int_{0}^{t}\lambda(x)dx$

,.

(3)

(4)

$N(t)$ は平均値関数 $A(t)$ 及び強度関数 $\lambda(t)=d\Lambda(t)/dt$ をもつ NHPP である．

更に本稿では，次の平均値関数を持つ NHPP のべき法則モデル (Power law model) を仮定する．
(5)$\Lambda(t;\eta, \beta)=-\log\overline{F}(t)=(\frac{t}{\eta})^{\beta}$

ここで $\overline{F}(t)=1-F(t)$ である．べき法則モデルにおいては $\beta>1$ の場合 $d\lambda(t)/dt>0$ となり，かつ$q(\infty)>0$が成り立つことから，最適予防保全間隔くなる現象を扱うことから，以降では$\beta>1$$t_{0}^{*}$は唯一つ定まる．時間の経過とともに故障が発生し易の場合のみを取り扱うこととする．

$t_{i}(i=1,2, \cdots, n)$ のような今，個の故障時刻列データがべき法則モデルに従う場合，尺度パラメータ$\eta(>0)$ と形状パラメータ $\beta(>1)$ の最尤推定値は，次の対数尤度方程式を最大にするようなパラメータの組$n$$(\hat{\eta},\hat{\beta})$として定義される．

$\log L_{NHPP}(\beta, \eta|t_{i})=\sum_{i=1}^{n}\log\lambda(t_{i};\eta, \beta)-\Lambda(t_{n};\eta, \beta)$式 (6) に対する一階の最適性条件により，$(\hat{\eta},\hat{\beta})$.

(6)

は次の同時方程式を解くことにょり得られる．

$\hat{\beta} = \frac{n}{n\log t_{n}-\sum_{i=1}^{n}\log t_{i}}$$\hat{\eta} = (\frac{t_{n}^{\hat{\beta}}}{n})^{\beta}4$

,

(7)

(8)

60
を式 (2) に代入することで，小修理を伴うブロックを算出することと最小総期待費用の点推定値

式 (7) と式 (8) より導出したパラメータの推定値取替え問題に関する最適な予防保全間隔$t_{0}^{*}$$(\hat{\eta},\hat{\beta})$$C(t_{0}^{*})$$C(t_{0}^{*})\wedge\wedge$$(t_{0}^{*}),$

が出来る．

3

パラメトリックブートストラップ法
2.2 節で説明した推定法では，小修理を伴うブロック取替え問題の解となる最適予防保全間隔の点推定値$t_{0}^{*}\wedge$を算出することが出来る．しかし，それは一本の故障時刻列データ $t_{i}(i=1,2, \cdots, n)$ に基づいて算出された推定量の近似的な期待値であり，確率変数としての推定量自身の不確実性を考慮することは出来ていない．そこで，本稿ではパラメトリックブートストラップ法を適用し，推定量の確率分布を導出することを考える．ブートストラップ法は元となる一つの故障時刻列データから複数の擬似故障時刻列データを複製する統計的手法であり，得られた擬似時刻列データをブートストラップ標本と呼ぶ．ブロック取替え問題に対して，本稿では元となる故障時刻列データ
目 $(1\leq k\leq m)$ の複製における$i$$t_{1},$$\cdots,$$t_{n}$から $m(>0)$ 組の故障時刻列データを複製する． 回$k$

番目の擬似故障時刻データ $t_{i}(i=1,2, \cdots, n)$ を組の擬似故障時刻列データから $m$ 組のモデルパラメータ$t_{k,i}^{*}$とおく．得られた$m$の最尤推定値を算出することが出来る．本稿では，ブートストラップ標本を生成するための方法としてシミュレーションベース及びリサンプリングベースの 2 種類の方法を用いる．更に，ブートストラップ法の適用に際し，NHPP に従うデータを直接用$\eta$と$\beta$いる場合と同次ボアソン過程 (HPP) へ変換したデータを用いる場合の二通りの手法について考察を行う．

HPP に従う擬似故障時間間隔データを $w_{k,i}^{*}(i=1,2, \cdots, n, k=1,2, \cdots,m)$ とおく．以下ではブートスト
ラップ法の種類別及びデータ形式の種類別ごとに，4 つのブートストラップ標本生成方法を紹介する．
方法 (i)NHPP データに対するシミュレーションベースのブートストラップ法 (NHPP- $BS$ 法 1) : 元となる故障時刻列データ $t_{i}(i=1,2, \cdots, n)$ により算出されたモデルパラメータの最尤推定値て，NHPP のべき法則モデルに従う擬似故障時刻列データ$i=1,2,$$\cdots,$ $n$$t_{k,i}^{*}$$\hat{\eta}$を生成する．ここで，$k=1,2,$である．具体的なシュミレーション方法としては，最初の故障時刻データ$t_{k,1}^{*}$と$\hat{\beta}$を用い$\cdots,$ $m$

及びが確率分布関数 $F(t_{k,1}^{*})=1-\exp\{-\Lambda(t_{k_{\}}1}^{*})\}$ に従うこと，及び 2 番目以降の故障時刻データ $t_{k,j}^{*}(j=2, \cdots, n)$ が確率分

布関数 $F(t_{k,j}^{*}|t_{k,J-1}^{*})=1-\exp[-\Lambda(t_{k,j}^{*})+\Lambda(t_{k,j-1}^{*})]$ に従うことを利用して，以下の式を逐次的に用いる．
$t_{k,1}^{*} = \hat{\eta}[-\log(U_{k}(0,1))]\not\supset 1$

(9)

,

(10)

$t_{k,j}^{*} = [-\hat{\eta}^{\hat{\beta}}\log(U_{k}(0,1))+t_{k,j-1}^{*\hat{\beta}}]^{\frac{1}{\beta}}$は区間 [0,1] での一様乱数を表す．式 (9) と式 (10) を用いて生成した擬似故障時刻列デー
ここで，を算出する．
タを用いて， 組のパラメータの最尤推定値$U_{k}(0,1)$$(\hat{\beta}_{k},\hat{\eta}_{k})(k=1,2, \cdots, m)$$m$

方法 (ii)NHPP データに対するリサンプリングベースのブートストラップ法 (NHPP- $BS$ 法 2) : 元となる
故障時刻列データ $t_{i}(i=1,2, \cdots, n)$ から，重複を許して無作為に 個の故障時刻データを抽出する．こう$n$して得たブートストラップ標本 $t_{ki}^{*}(i=1,2, \cdots, n, k=1,2, \cdots, m)$ から $m$ 組のパラメータの最尤推定値$(\hat{\beta}_{k,\hat{\eta}_{k}})(k=1,2, \cdots, m)$を算出する．

変換データに対するシミュレーションベースのブートストラップ法 (HPP- $BS$ 法 1) : パラメータの指数分布に従う擬似故障時間間隔列データを指数分布の逆関数を用いて生成し，得られ$m$ 及び $i=1,2,$た時間間隔データを HPP に従う時刻列データへと変換する．ここで，$k=1,2,$
方法$(iii\rangle HPP$$\lambda=1$$w_{ki}^{*}$$\cdots,$ $n$$\cdots,$である．その後，元となる故障時刻列データ $t_{i}(i=1,2, \cdots, n)$ から算出されたモデルパラメータの推定値と平均値関数の逆関数を用いて，NHPP のべき法則モデルに従う擬似故障時刻列データを行う．これら擬似故障時刻列データを用いて， 組のパラメータの最尤推定値及び$\hat{\beta}$$m$を算出する．

$t_{ki}^{*}$$\hat{\eta}$

へと変換$(\hat{\beta}_{k},\hat{\eta}_{k})(k=1,2, \cdots, m)$

61
方法 (iv)HPP 変換データに対するリサンプリングベースのブートストラップ法 (HPP- $BS$ 法 2) : 元となる故障時刻列データ $t_{i}(i=1,2, \cdots, n)$ を算出されたモデルパラメータの推定値$\hat{\eta}$

及びとべき法則モデ$\hat{\beta}$

ルの平均値関数を用いて，平均値関数列データ (HPP データ) へと変換する．それを時間間隔データへと変換した後，重複を許して無作為に 個の故障時間間隔データを抽出する．こうして得られたブートスト$(i=1,2, \cdots , n, k=1,2, \cdots, m)$ を逆の手順で NHPP のべき法則モデルに従う擬似故障時
ラップ標本
$n$

$w_{ki}^{*}$

刻列データ$t_{ki}^{*}$へと再変換を行う．これら擬似故障時刻列データを用いて， 組のパラメータの最尤推定値$m$$(\hat{\beta}_{k,\hat{\eta}_{k}})(k=1,2, \cdots, m)$を算出する．

以上 4 つの方法により得られた $m$ 組のパラメータの最尤推定値と式 (2) を用$m$ 組の最適予防保全間隔
いて，の推定量を求めることが出来る．更に式 (1) を用いて，対応する最小総$(\hat{\beta}_{k},\hat{\eta}_{k})(k=1,2, \cdots, m)$$t_{0}^{*}$

期待費用$C(t_{0}^{*})$の推定量を $m$ 組算出し，これらを用いて経験分布を構成することでの確率分布を導出する．続いて，得られた最適予防保全間隔$k=1,2,$
こで，待費用$C(t_{0}^{*})$

は，$t_{0}^{*}[m/2]$$\cdots,$ $m$

であり，$t_{0}^{*}[k]$$t_{0}^{*}$$t_{0}^{*}$

及びの $m$ 組の推定量をそれぞれ$C(t_{0}^{*})$$t_{0}^{*}[k]$の推定量とおく．こは $t_{0}^{*}[1]\leq t_{0}^{*}[2]\leq\cdots\leq t_{0}^{*}[m]$ を満たすものとする．同様に，最小総期の $m$ 組の推定量をそれぞれ $C(t_{0}^{*}\}[k](C(t_{0}^{*})[1]\leq C(t_{0}^{*})[2]\leq\cdots\leq C(t_{0}^{*})[m])$を $BS$ 中央値，$\overline{t}_{0}^{*}=(\sum_{k=1}^{m}t_{0}^{*}[k])/m$

とおく．本稿でを $BS$ 平均値と呼ぶ．更に，優位水準を $\alpha\in(0,1\}$ とした

場合の 100 $(1-\alpha$ %信頼区間を $[t_{0}^{*}[m(\alpha/2)\cdot], t_{0}^{*}[m(1-\alpha/2)]]$ によって算出する．最適予防保全間隔の推定量の分散 V

歪度

$S$

尖度 $K$ はそれぞれ$V = \frac{\sum_{k=1}^{m-1}(t_{0}^{*}[k]-\overline{t}_{0}^{*})^{2}}{m-1}$$S = \frac{\sum_{k=1}^{m}(t_{0}^{*}[k]-\overline{t}_{0}^{*})^{3}}{mV^{\frac{3}{2}}}$$K = \frac{\sum_{k=1}^{m}(t_{0}^{*}[k]-\overline{t}_{0}^{*})^{4}}{mV^{2}}$,,(11)(12)(13)によって求める．

4

数値例

4.1

シミュレーション結果

NHPP のべき法則モデルのパラメータを $(\beta, \eta)=(3,0.2)$ と設定し，元データとなる故障時刻列データを生成した．費用の係数はそれぞれ $c_{1}=1,$ $c_{2}=1000$ とする．この場合，小修理を伴うブロック取替え問題の最適予防保全間隔と最小総期待費用の真の値はそれぞれ $t_{0}^{*}=1.5S7401,$ $C(t_{0}^{*})=944.9407$ となる．更に，元データに対して 3 節で説明した 4 通りのブートストラップ法を適用し，$m=10000$ 組のブートストラップ標本を生成する．10000 組のブートストラップ標本それぞれに対し，最尤法を適用してパラメータ推定を行い，10000 個の最適予防保全間隔及び最小総期待費用の推定量を算出する．それら推定量に基づいて経験分布を構成し，最適予防保全間隔及び最小総期待費用の推定量の確率分布に関する分散・尖度・歪度を求める．また，優位水準 $\alpha=0.05$ として，最適予防保全間隔及び最小総期待費用の 95 %信頼区間を導出する．
よく知られているように，尖度が 3 に近く，かつ歪度が 0 に近い場合，推定量の確率分布を近似的に正規分布とみなすことが出来る．
元となる故障時刻列データの個数 の変化による影響を調べるため，$n=5$ から $n=245$ まで 5 刻みで増加させた場合の $BS$ 推定値 (平均と中央値) や高次モーメント (分散、 尖度、 歪度) 及び 95 %信頼区間$n$$BS$ l 及び $BS$ 2 を適用した場合の最適を計算する．表 1-表 8 は NHPP データ及び HPP 変換データに対し，予防保全間隔と対応する最小総期待費用に関する結果を示している．歪度と尖度の値に注目してみれば，最小総期待費用の推定量の分布では，元データ数が少ない場合でも歪度が 0 に近く，尖度が 3 に近い値をとっている．これに比べて，最適予防保全間隔の推定量の分布は歪度・尖度ともに 0 と 3 よりも大きな値を取る

62
傾向がある．このことから，最小総期待費用の推定量に関しては漸近的に正規分布として近似することが可能であるが，最適予防保全間隔の推定量に関しては正規近似が適用出来ないという結論が得られる．更に，最適予防保全間隔及び最小総期待費用の推定量に関する漸近的変化の様子を図 1-図 8 に示す．これらの図から，元データ数を増やすほど信頼区間はより狭くなること，及び，分布の平均値と中央値は最尤推定値に近づいていくことが確認できる．しかし，データ数を増やしたとしても分布の平均値や中央値が最尤推定値と完全に一致することはない，すなわち単一系列の故障データから一致推定量を構成することは出来ないことが知られている．更に分布の平均値及び中央値は，真の値に近い値付近で推移していることも確認できるが，こちらも完全に一致することはなく，4 つの手法いずれが点推定の観点から優れているかを結論付けることは出来ない．続いて，各手法における 95 %信頼区間の大きさを図 9-図 10 に示す．これら及び最小総期待費用に関する 95の図では，縦軸が最適予防保全間隔に関する 95 %信頼区間の大きさ%信頼区間の大きさを表している．結果から，いずれの推定量に関しても 95 %信頼区間の幅が狭いという意味で NHPP- $BS$ 法 2 が最も優れた手法であると言える．
$t_{[95|}^{*}$$C_{[95]}(t^{*})$

4.2

実データ解析

ディーゼルエンジンの小修理データに対して，本研究の手法を適用した例を示す．本研究では，Meekerand Escobar [4] において紹介されていた実故障データ ( $n=71$ 個，平均故障時間間隔 0. $359hr$ ) を用いた．
これらのデータを用いて算出したモデルパラメータの最尤推定値は $(\hat{\beta},\hat{\eta})=(2.76,5.45)$ であった．費用の係数はそれぞれ $c_{1}=1,$ $c_{2}=1000$ とする．元データに対して，95%信頼区間の観点で最も優れていた$m=10000$ 組のブートストラップ標本を生成して経験分布を構成する．表 9 は，NHPP- $BS$ 法 2 を適用し，ブートストラップ法を用いた場合の最適予防保全間隔と対応する最小総期待費用に関する結果を示している．この表から，いずれの推定量に関しても，分布の平均値及び中央値が最尤推定値とかなり近い値をとっていることが確認できる．更に図 11-図 12 は，それぞれの場合における経験分布を表している．表に示した歪度及び尖度の値と図の形状から，最小総期待費用の推定量の分布は比較的正規分布に近い形状を保っているが，最適予防保全間隔の推定量の分布は正規分布とはかなり離れた右に裾の長い分布形状となっていることが分かる．

5

まとめと今後の展望

本稿では，小修理を伴うブロック取替え問題の解である最適予防保全間隔と最小総期待費用の推定量の確率分布を導出し，高次モーメントや信頼区間を得るためにブートストラップ法を適用した．適用手法として具体的には，NHPP に従うデータを直接用いる場合や HPP へのデータ変換を用いる場合を提案し，それぞれに対してシミュレーションベースのブートストラップ法及びリサンプリングベースのブートストラップ法を適用した．シミュレーション実験や実データ解析を通じて，様々な推定量に関する統計的性質を確認することが出来た．今後の課題としては，ノンパラメトリックブートストラップ法を用いて，同様に信頼区間の導出を行う予定である．

参考文献[1] R. E. Barlow, and F. Proschan (1965), Mathematical Theory of Reliability, Wiley.
[2] T. Nakagawa (2005), Maintenance Theory of Reliability, Springer.

63[3] S. Tokumoto, T. Dohi and W. Y. Yun (2012), Confidence interval of optimal age replacement policy-parametric bootstrapping approach-, Advanced Reliability and Maintenance Modeling V ( $H.$Yamamoto, C. Qian, L. Cui and T. Dohi, eds.), pp.494-501, McGraw Hill.
[4] W. Q. Meeker, and L. A. Escobar (1998), Statistical Methods for Reliability Data, Wiley.

表 1: NHPP- $BS$ 法 1 を用いた場合の最適予防保全間隔に関する結果．

表 2: HPP- $BS$ 法 1 を用いた場合の最適予防保全間隔に関する結果．

64

表 3:NHPP BS 法 2 を用いた場合の最適予防保全間隔に関する結果．

表 4: HPP- $BS$ 法 2 を用いた場合の最適予防保全間隔に関する結果．

表 5: NHPP- $BS$ 法 1 を用いた場合の最小総期待費用に関する結果．

65

表 6: HPP- $BS$ 法 1 を用いた場合の最小総期待費用に関する結果．

表 7: NHPP- $BS$ 法 2 を用いた場合の最小総期待費用に関する結果．

表 8: HPP- $BS$ 法 2 を用いた場合の最小総期待費用に関する結果．

66

図 1: NHPP- $BS$ 法 1 を用いた場合の最適予防保全間図 2: HPP- $BS$ 法 1 を用いた場合の最適予防保全間隔隔の推定量の漸近的振る舞いの推定量の漸近的振る舞い．

図 3: NHPP- $BS$ 法 2 を用いた場合の最適予防保全問図 4: HPP- $BS$ 法 2 を用いた場合の最適予防保全間隔隔の推定量の漸近的振る舞いの推定量の漸近的振る舞い．

67

図 5: NHPP- $BS$ 法 1 を用いた場合の最小総期待費用図 6: HPP- $BS$ 法 1 を用いた場合の最小総期待費用のの推定量の漸近的振る舞い推定量の漸近的振る舞い．

図 7: NHPP- $BS$ 法 2 を用いた場合の最小総期待費用図 8: HPP- $BS$ 法 2 を用いた場合の最小総期待費用のの推定量の漸近的振る舞い

推定量の漸近的振る舞い．

68

図 9: 最適予防保全間隔の推定量に関する 95 %信頼図 10: 最小総期待費用の推定量に関する 95 %信頼区区間の大きさ間の大きさ．

図 11: NHPP- $BS$ 法 2 を用いた場合の最適予防保全問図 12: NHPP- $BS$ 法 2 を用いた場合の最小総期待費用隔の推定量の分布の推定量の分布．

表 9: NHPP- $BS$ 法 2 を用いた場合の実データ解析結果．