数理解析研究所講究録 1297 巻 2002 年 179-185

179
Comparison of Estimation Algorithms for NHPP-Based
Software Reliability Models
岡村寛之, 渡部保博, 土肥正

Hiroyuki Okamura, Yasuhiro Watanabe and Tadashi Dohi
Department of Information Engineering, Graduate School of Engineering,
Hiroshima University, Japan

1

はじめに

現在, 高度情報化社会あるいは情報ネットワーク社会であると言われているように, コンピュータシステムはあらゆる分野で重要な役割を担っている. 特に金融や通信といった社会の基盤となる分野において, 業

務支援・情報戦略などの主要業務におけるコンピュータシステムの必要性が加速度的に増大している.

コンピュータシステムへの依存性が高まるにつれて, コンピュータシステムの構成要素であるソフトウエアの信頼性確保に対する十分な検討が必要とされている.

現在までにソフトウェアの信頼性を定量的に評価する方法が数多く提案されている. 一般的によく知られている手法は, テスト段階あるいは運用段階におけるフオールト発見事象を確率過程を用いて表現する手法であり, フォールト除去作業によってソフトウェアの信頼度が向上することから, これらの確率過程はソフトウェア信頼度或長モデルと呼ばれている. 代表的なソフトウェア信頼度成長モデルとしては,
Jelinski-Moranda モデル [1] や指数形信頼度成長モデル [2] 等がある.
本論文では, 非同時ポアソン過程 (NHPP) に基づいたソフトウエア信頼性モデルにおけるパラメータ
アルゴリズム [3,4,5,6] に基づいた推定手法を提案し, 他の推定アルゴリズムとの性推定手法として,能比較を行う.
$\mathrm{E}\mathrm{M}$2NHPP に基づいたソフトウエア信頼性モデル
NHPP によるソフトウェア信頼度成長モデルでは, 任意の時刻までに発見されるフオールト数が次の確率関数に従う計数過程 $\{N(t), t\geq 0\}$ によって表現される.
$\mathrm{P}\mathrm{r}\{N(t)=n\}=\frac{\{\Lambda(t)\}^{n}}{n!}\exp\{-\Lambda(t)\}$,$n=0,1,2,$$\ldots$.

(1)

は NHPP の平均値関数と呼ばれ, 時間区間 $(0, t]$ において発見される総期待フオールト数を表す. さらに, $\lambda(t)=dAQ)/dt$ ま NHPP の強度関数と呼ばれ, テスト時刻 \dagger こおける瞬間
フォールト発見率 (ソフトウェア故障率) を表す. NHPP で記述されるソフトウエア信頼度成長モデルでの設定により様々なパターンを表現することが可能である. 特に有限のフオールトには, 平均値関数対するデバッギング過程を記述するモデルは, 一般化順序統計量モデルによって統一的な取扱いが可能で
ここで$\Lambda(t)$$t$$\mathrm{E}[N(t)]$$\Lambda(t)$あることが知られている $[7,8]$ .

一般化順序統計量モデルは, 次の仮定の下で議論される.

(仮定 A) ソフトウェアのテスト中にソフトウェア故障が発生した場合, その原因となるフオールトは瞬間的に発見・除去される.

(仮定 B) プログラム中に含まれる初期フオールト数$N_{0}$

は平均$\mathrm{t}d(>0)$のポアソン分布に従う.

(仮定 C) ソフトウェア故障はおのおの独立かつ時間に関してランダムに発生し, 各ソフトウエアフオー
ルトの発見時刻は確率分布 $F(t)$ , 密度関数 $f(t),$ $t\geq 0$ に従う非負で連続の確率変数によって記述される.

いま,をソフトウエア故障の発生時刻を表す確率変数列として定義する. このとき,が互いに独立で同一の確率分布 $F(t)$ に従う確率変数であることを意味する. 初期時刻 $t=0$$\{X:, i=1, \ldots, N_{0}\}$は仮定においてソフトウェア内にトが発見される確率は$\mathrm{C}$$X_{:}$$n$個のフォールトが潜在している条件の下で, 時刻$\mathrm{P}\mathrm{r}\{N(t)=m|N_{0}=n\}$

$=$

$t$

までに $m$ 個のフオール$(\begin{array}{l}nm\end{array})\{F(t)\}^{m}\{1-F(t)\}^{n-m}$(2)

180となる. 仮定より初期フォールト数はパラメータれるフォールト数の確率分布は$\mathrm{B}$$\omega$のポアソン分布に従うため, 時刻$\mathrm{P}\mathrm{r}\{N(t)=m\}=\frac{\{\omega F(t)\}^{m}}{m!}\exp$$t$までに発見さ{-。F(t)}(3)をもつ NHPP と等価であり, 有限フォールトをもつ既存のソフと表現される. これは平均値関数
トウェア信頼性モデルに対する一般的な表現となっている.
$\omega F(t)$3パラメータ推定ソフトウェア信頼性モデルを実際のソフトウェア開発工程で用いるためには, ソフトウェアフォールトデータを用いてモデル内部のパラメータ ( で表記する) を推定する必要がある. 最も一般的な推定手法は最尤法であり, 最尤法によるパラメータの推定値は時刻 までに観測された $m(>0)$ 個のフォールトデータ$D_{t}=(x_{1}, \ldots, x_{m})$ に対する対数尤度関数$\theta$$t$(4)$\log L_{\mathrm{N}\mathrm{H}\mathrm{P}\mathrm{P}}(\theta|D_{t})=\sum_{i=1}^{m}\log\lambda(x_{i})-\Lambda(x_{m})$を最大にするパラメータとして決定される. つまり, 以下の非線形連立方程式で表現される同時尤度方程式を解くことによってパラメータの推定値を算出する.
$\frac{\partial}{\partial\theta}\log L_{\mathrm{N}\mathrm{H}\mathrm{P}\mathrm{P}}(\theta|D_{t})=0$.

(5)

上記の非線形連立方程式を解く方法として, 最も有名なものは Netwon 法である. これは, 対数尤度関数の 1 階微分と 2 階微分をそれぞれ$u( \theta)=\frac{\partial}{\partial\theta}\log L_{\mathrm{N}\mathrm{H}\mathrm{P}\mathrm{P}}(\theta|D_{t})$

,

(6)

(7)

$H( \theta)=\frac{\partial^{2}}{\partial\theta^{2}}\log L_{\mathrm{N}\mathrm{H}\mathrm{P}\mathrm{P}}(\theta|D_{t})$とおくと, 現在の推定値$\theta_{0}$

から

(8)

$\hat{\theta}=\theta_{0}-H^{-1}(\theta_{0})u(\theta_{0})$を Fisher の情報行列の逆の更新を繰り返すことで最尤推定値を算出する. Newton 法において一行列にしたものが Fisher’s scoring 法と呼ぼれる方法である.
Newton 法や Fisher’s scoring 法の欠点は式 (5) の非線形連立方程式の解が必ずしもソフトウェア信頼性モデルにおいて有効な解となり得ない点にある. つまり, 推定に用いる観測データの種類によっては, パラメータに対する (暗黙の) 制約を違反する可能性がある.
この問題を回避するため, 文献 [9] では NHPP に基づいたソフトウェア信頼性モデルのパラメータ推定に EM アルゴリズムを用いた手法を提案している. 以下ではその概要を述べる.
EM アルゴリズムとは, 観測不可能な変数を含む推定問題において用いられる手法である. いま, 一般が観に対して確率変数的に確率変数 $X$ の標本値が観測不可能であり, 任意の可測関数をも測可能であるとき, パラメータ集合 を推定する問題を考える. 確率変数 $X$ が確率密度関数
つとき, EM アルゴリズムの $k+1$ 番目のステツプにおける推定値は 番目のステツプで推定されたパラメータを適用したときの期待対数尤度を最大にするパラメータによって決定される. すなわち$H^{-1}$ $()$$u(\cdot)$$X_{\mathrm{o}\mathrm{b}\mathrm{s}}=u(X)$$f(\cdot;\xi)$

$\xi$

$k$

$\{\mathrm{E}[\log f(X;\xi)|x_{\mathrm{o}\mathrm{b}\mathrm{s}};\hat{\xi}^{(k)}]\}$$\hat{\xi}^{(k+1)}=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}\xi$

ここで,.

(9)

番目のステップにおけるパラメータの推定値であはそれぞれ観測された標本値とはパラメータを用いは括弧中の値が最大値になる引数を返す関数である. また,り,たときの期待値演算子を表す.
(は一般化順序統計量モデルのパラメータを推定するために, 総数が個の確率変数’観測不可能な確率変数) の順序統計量 $X_{[1]}<\cdots<X_{[N_{0}]}$ を考える. いま, ソフトウェア内に潜在するフオー$x_{\mathrm{o}\mathrm{b}\mathrm{s}}$と$\hat{\xi}^{(k)}$

$k$

$\xi$

$\mathrm{E}[\cdot;\xi]$$\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}\{\cdot\}$$N_{0}$$X_{1},$$\ldots X_{N_{\mathrm{O}}}\backslash$$N_{0}$181
ルトの総数がという条件の下で, すべてのフォールトに対する発見時刻データが与えられたならば, その尤度関数および対数尤度関数は$N_{0}\ovalbox{\tt\small REJECT} n$$L( \omega, \theta|D_{\infty})=\omega^{n}\exp\{-\omega\}\prod_{i=1}^{n}f(x_{i}; \theta)$$D_{\ovalbox{\tt\small REJECT}}\ovalbox{\tt\small REJECT}(x\mathrm{t},$$\ovalbox{\tt\small REJECT}$

.

(10)

,

$\log L(\omega, \theta|D_{\infty})=n\log\omega-\omega+\sum_{i=1}^{n}\log f(x_{i}; \theta)$$x\mapsto$$,$.

(11)

ソフトウエア内に潜在するすべてのフオールトの発見時刻データが与えられた場合, パとなる. また, 各
ラメータの推定値は通常の最尤法より観測されたフォールトの総数, つまりの推定量は対数尤度関数
ソフトウェアフォールトの発見時間分布に関するパラメータを最大にするものとして与えられる.
EM アルゴリズムは一部のフォールトの発見時刻に対する観測値が得られた下で構築される. すなわ
ち, $m<n$ に対して不完全なフォールトデータ $D_{t}=(x_{1}, \ldots, x_{m})$ が得られた時の EM アルゴリズムを考える. 式 (9) と式 (11) から, 一般化順序統計量モデルにおける期待対数尤度はとなる. 従って,$\hat{\omega}=n$$\omega$$\sum_{=1}^{n}\dot{.}\log f(x_{i}; \theta)$$\theta$$\log L_{E}(\omega, \theta|D_{t})=\mathrm{E}[N_{0}|D_{t}; \omega, \theta]\log\omega-\omega+\mathrm{E}[\sum_{i=1}^{N_{0}}\log f(X_{i};\theta)|D_{t}$であり, これを最大にする
テップの更新式は$\theta$$\omega,$;$\omega,$(12)$\theta]$をパラメータの更新値とする. 従って, EM アルゴリズムによる $k+1$ ス$\ovalbox{\tt\small REJECT}^{(k+1)}=\mathrm{E}[N_{0}|D_{t}; \hat{\omega}^{(k)},\hat{\theta}^{(k)}]$

,

(13)
(14)

$\hat{\theta}^{(k+1)}=\arg\max_{\wedge}(\mathrm{E}|\dot{\sum}\log f(X_{i};\theta)|D_{t};\hat{\omega}^{(k)},\hat{\theta}^{(k)}|\}$$\theta$$($$\lfloor_{\overline{i=1}}$$|$$\rfloor/$

として与えられる. ここで, 各フオールト発見時刻分布に具体的な確率分布を適用することでパラメータ推定のための更新式を得ることができる.

4

推定手法の評価

前節で紹介した各推定手法の有効性を, 実際のソフトウェア開発工程で得られたフォールトデータを用いて数値的に検証する. 実際のソフトウェアフォールトデータとして, 文献 [10] に記載されている 6 種類のデータ(データ数 86), DS2(データ数 136), DS3(データ数 207), DS4(データ数 394), DS5(データ数123), DS6(データ数 97) を用いる.
はじめに, 単一のフオールト発見時刻が指数分布に従う場合に関して, Newton 法および Fisher’s scoring法との解の安定性の比較を行う. 図 1 および図 2 は DS1 において異なる 2 種類の初期値から推定を行ったにお場合の, 各推定手法による解の収束の振る舞いを示したものである. 同様に, 図 3 および図 4 はいて異なる 2 種類の初期値から推定を行った場合の, 各推定手法による解の収束の振る舞いを示したものである. これらの図から, Newton 法および Fisher’s sconng 法による推定は初期値に依存しており, 常にアルゴリズムは全ての初期値に対して常に安定して解安定して解を得ることができていない. 一方,を得ていることがわかる.
$500,$,
次に各推定手法の性能を定量的に評価するために, 100 個の初期値 $(\omega=50,100,$, I.OE-03) からそれぞれのアルゴリズムを 10 回試行したときの推定値を算出し, 実際の最尤
0.2E-03,推定値との相対誤差を比較した. 相対誤差は以下の式で定義される.
$\mathrm{D}\mathrm{S}1$$\mathrm{D}\mathrm{S}2$$\mathrm{E}\mathrm{M}$$\lambda=0.1\mathrm{E}- 03$$\ldots,$$\ldots$(相対誤差)$= \frac{(\mathrm{a}\mathrm{e}_{\acute{i\Xi}\mathrm{f}\mathrm{f}\mathrm{i}\mathrm{L})-(\text{最}\lambda:\text{推_{}\acute{i\mathrm{E}}\mathrm{f}_{\mathrm{L}}^{\mathrm{g})}}}}{(\text{最}\lambda;3\mathbb{E}\acute{j\mathrm{E}}\{\mathrm{i}\Xi)}\mathrm{x}100(\%)$

.

(15)

表 1-3 は, 単一のフオールト発見時刻が指数分布, 二次のガンマ分布, 形状パラメータ 2 のワイブル分布に従う場合の結果を示している. 特に推定値を得る過程でパラメータの制約違反が発生した場合, その時点で計算を終了し, 相対誤差の評価には用いていない. また, 表中の convergence は制約違反が発生しなかった割合を示している. これらの表から, 前述の結果と同様に, EM アルゴリズムによる推定手法
では安定して解が得られるが, 全体的に収束スピードが遅い (相対誤差が大きい) という性質があることがわかる. 他方 Newton 法と Fisher’s scoring 法に関しては, 最終的に推定値が算出できるかどうかは初期値に大きく依存してることがわかる.

182
5

まとめ

本論文では, 非同時ポアソン過程 (NHPP) に基づいたソフトウェア信頼性モデルにおけるパラメータ推定手法として, Newton 法, Fisher’s scoring 法および EM アルゴリズムに基づいた推定手法を紹介し, それらの性能を数値的に検証した.

6

参考文献

[1] Jelinski, Z. and Moranda, P. B., Software reliability research, In Statistical Computer Performance
Evaluation (Ed. by Freiberger, W.), PP. 465-484, Academic Press, New York (1972).
[2] Goel, A. and Okumoto, K., Time-dependent error-detection rate model for software reliability andother performance measures, IEEE Ikans. Reliab., Vol. R-28, pp. 206-211 (1979).
[3] Dempster, A. P., Laird, N. M. and Rubin, D. B., Maximum likelihood from incomplete data via thealgorithm, J. $Roy$ . Statist. Soc. Ser. , Vol. 39, PP. 1-38 (1977).
$\mathrm{E}\mathrm{M}$$B$[4] Wu, C. F. J., On the convergence properties of the(1983).

$\mathrm{E}\mathrm{M}$algorithm, Ann. Statist., Vol. 11, PP. 95-103[5] Redner, R. A. and Walker, H. F., Mixture densities, maximum likelihood and the EM algorithm,
SIAM Review, Vol. 26, PP. 195-239 (1984).

[6] 宮川雅巳, EM アルゴリズムとその周辺, 応用統計学, Vol. 16,$\mathrm{P}\mathrm{P}\cdot 1-19$(1987).

[7] Langberg, N. and Singpurwalla, N. D., Unification of some software reliability models, SIAM J. Sci.
Comput., Vol. 6, pp. 781-790 (1985).
[8] Joe, H., Statistical inference for general order statistics and nonhomogeneous Poisson process software
reliability models, IEEE Trans. Software Eng., Vol. SE-15, pP. 1485-1490(1989).
アルゴリズムに基づいたソフトウエア信頼性モデルの推[9] 岡村寛之, 渡部保博, 土肥正, 尾崎俊治,定, 電子情報通信学会論文誌 (A), Vol. J85-A,. 442-450(2002).
$\mathrm{E}\mathrm{M}$$\mathrm{P}\mathrm{P}$[10] Lyu, M. R. (ed.), Handbook of Software Reliability Engineering, McGraw-Hill, New York, NY(1996).

183

Figure 1: Behavior of estimats for DS1 $(\{\omega, \lambda\}=\{30,0.000015\})$Figure 2: Behavior of estimats for DS1 $(\{\omega, \lambda\}=\{250,0.00003\})$Figure 3: Behavior of estimats for $\mathrm{D}\mathrm{S}2(\{\omega, \lambda\}=\{50,0.00005\})$ .

184

Figure 4: Behavior of estimats for$\mathrm{D}\mathrm{S}2$$(\{\omega, \lambda\}=\{275,0.00012\})$.

Table 1: Comparison with estimation methods (Exponential SRM).

data set, )(MLE:$\omega$

methodconvergence

relative errors

EM
NewtonFisherEM
NewtonFisherEM
NewtonFisherEMNewton
Fisher
EMNewton
FisherEM
NewtonFisher

100%
19%
66%
100%
44%
82%
100%

(6.60%, 10.17%)(8.54%, 6.81%)(0.24%, 0.23%)(2.54%, 3.81%)(0.04%, 0.65%)(0.05%, 2.85%)(33.50%, 27.51%)$\lambda$

DS1(107.59, 1.56E-05)

DS2(142.88, 3.42E-05)

DS3(256.76, 9.86E-05)

DS4(401.06, 4.22E-05)

DS5(133.13, 3.90E-05)

DS6(118.23, 1.37E-04)$(\omega, \lambda)$

97%

72%
100%
84%
64%
100%
45%
90%
100%
81%
78%

(0.71%, 1.06%)(0.07%, 1.07%)(0.00%, 0.00%)(2.11%, 1.58%)(1.03%, 1.40%)(0.00%, 0.00%)(1.67%, 1.93%)(0.03%, 9.61%)(77.79%, 48.69%)(0.33%, 0.24%)(0.23%, 1.43%)

185
Table 2: Comparison with estimation methods (Gamma SRM).
data set, )(MLE:$\omega$methodconvergence

DS1(90.06, 4.74E-05)

DS2(136.99, 7.90E-05)

DS3(216.36, 2.95E-04)

DS4(397.28, 8.83E-05)

DS5(129.51, 8.66E-05)

DS6(106.58, 3.64E-04)

relative errors$(\omega, \lambda)$$\lambda$

EM
NewtonFisherEMNewton
FisherEM
NewtonFisher
EMNewton
Fisher
EMNewton
Fisher
EMNewton
Fisher100%
39%
83%100%50%
84%
100%
65%
26%
100%
100%
77%
100%
50%
85%
100%
37%
46%

(17.61%, 24.90%)(34.53%, 2.87%)(0.00%, 0.23%)(0.63%, 0.04%)(0.00%, 0.00%)(0.00%, 0.00%)(5.20%, 8.78%)(0.00%, 0.00%)(0.03%, 1.97%)(0.06%, 0.00%)(0.00%, 0.00%)(0.00%, 0.00%)(0.35%, 0.01%)(0.06%, 0.01%)(0.00%, 0.00%)(10.69%, 16.36%)(0.42%, 0.04%)(0.00%, 0.00%)

Table 3: Comparison with estimation methods (Weibull SRM).

data set, )(MLE:$\omega$

methodconvergenceDS1(87.25, 4.03E-l0)$\mathrm{D}\mathrm{S}2^{\cdot}$(136.17, 8.48E-l0)

DS3(209.28, 1.63E-08)

DS4(397.02, 8.22E-l0)

DS5(129.08, 9.24E-l0)

DS6(105.13, 1.92E-08)

relative errors$(\omega, \lambda)$$\lambda$

EM
NewtonFisherEM
NewtonFisherEM

NewtonFisherEM
NewtonFisherEM
NewtonFisherEM
NewtonFisher

100%
32%
0%
100%
46%
0%
100%
45%
0%
100%
86%
0%
100%
48%
0%
100%
32%0%

(0.00%, 0.00%)(0.04%, 0.10%)(–, –)(0.00%, 0.00%)(3.34%, 5.47%)(–, –)(0.00%, 0.00%)(1.09%, 2.49%)(–, –)(0.00%, 0.00%)(0.44%, 1.77%)(–, –)(0.00%, 0.00%)(2.72%, 5.08%)(–, –)(0.00%, 0.00%)(1.39%, 3.50%)(–, –)