確実に再現しないバグの修正確認

今日は午後から現実逃避のプログラミング。本当は人事評価とか事務的作業があるのですが、他人の評価は気が重い。

スレッド間の協調に関する障害等タイミングや他のプロセスとの通信や負荷によって発生したりしなかったりするバグの場合、修正してから動作確認の際に充分な試行回数を取って障害が収まったことを保証しなくてはいけません。

このようなバグによる障害発生が確率過程であると仮定すれば、障害の発生確率(p)が決まれば、ベイズ推定により n回の試行で障害再現しない時にバグが残っている確率(p_bug(n))の収束計算ができます。こんな感じかな。

 p_bug(n) = (1 - p)*p_bug(n-1) / (1 - p*p_bug(n-1))

これで求める品質水準(ε)を定めて、 p_bug(n) < ε となる最小の n を求めればよいでしょう。

ところでこの障害の発生確率は、従来のバージョンでの試行毎に障害発生 n1 回、正常処理 n2 回として n1/(n1+n2) としてよいのでしょうか。
試行回数が充分あればほぼ実際の確率に近づくと思われますが、検査環境の利用に制限がある場合などであまり試行回数が取れなかった場合は障害発生確率も確率分布から望む水準での最低確率を推定することができます。

障害発生確率を p (0 <= p <= 1) として、試行結果(n1, n2)を得る確率 P(n1, n2) は次の性質を持ちます。

 P(n1,n2) ∝ p**n1 * (1-p)**n2

実際は比例定数として(n1+n2)からn1を選ぶ組み合わせ (n1+n2)Cn1 を掛けることになりますが、今は必要ありません。

従って試行結果 (n1, n2) を得た時に障害発生確率 p = x である確率が f(x)*dx であるような確率分布は

 f(x) = A*(p**n1*(1-p)**n2)

f(x) の不定積分 F(x) とすると p が a < p < b の範囲である確率は F(b) - F(a) なので、F(0) = 0 および F(1) = A = 1 とします。
従って、求める水準を ε として、 F(x) < ε となる x=xmin が p の最低値になります(試行結果(n1 n2) において p が xmin以下である確率が F(xmin) < ε)。この時、試行結果が得られる確率がε以下であるとは限りません。

というわけで最低確率を推定するスクリプトを作成。

 # p**n1*(1-p)**n2 の積分値を返す Proc を返す
 def integral(n1, n2)
   func = Proc.new do |x| (x**n1)*((1.0-x)**n2) end
 
   trape = Proc.new do |x, n|
     h = x / n.to_f
     # func.call(0) == 0 なので省く
     sum = func.call(x) * 0.5
     (n-1).times do |i|
       sum += func.call(h*(i+1).to_f)
     end
     sum*h
   end
   Proc.new do |x|
     prev = sum = trape.call(x, 1)
     20.times do |i|
       sum = trape.call(x, 2**(i+1))
       if((sum - prev).abs < prev*1e-5)
         break
       end
       prev = sum
     end
     sum
   end
 end
 
 def combine(n1, n2)
   m = 1
   ((n1+1)..(n1+n2)).each do |i|
     m *= i
   end
   n = 1
   (1..n2).each do |i|
     n *= i
   end
   m/n
 end
 
 if($0 == __FILE__)
   n1 = ARGV.shift.to_i
   n2 = ARGV.shift.to_i
   acc = ARGV.shift.to_f
 
   poly = integral(n1, n2)
   total = poly.call(1.0)
 
   # f(p) = p**N*(1-p)**L とすると
   # 0 <= p <= 1 において
   # f(p) >= 0
   # f(0) = f(1) = 0
   # f'(p) = p**(N-1)*(1-p)**(L-1)*(N - (N+L)p)
   # より上に凸で変曲点を p = N / (N+L) の一つだけ持つ
   # 上限をこの変曲点とする
   bottom = 0.0
   top = n1.to_f / (n1+n2).to_f
   middle = top/2.0
 
   20.times do |i|
     tmp = poly.call(middle)/total
     if(tmp > acc)
       top = middle
     elsif(tmp < acc)
       bottom = middle
     else
       top = bottom = middle
       break
     end
     middle = (bottom+top)/2
   end
 
   puts "推定障害発生確率: #{middle} (誤差幅: #{top-bottom})"
   puts "この推定値での #{n1}回障害 #{n2}回正常の確率は: #{middle**n1*(1.0-middle)**n2*combine(n1, n2)}"
 end

実行結果

 $ ruby estimate.rb 3 5 0.05
 推定障害発生確率: 0.168750464916229 (誤差幅: 3.57627868652344e-07)
 この推定値での 3回障害 5回正常の確率は: 0.106802160491759

8回の試行で 3回障害、5回正常の結果が得られた場合、障害発生確率 p < 0.168750464916229 である確率が約 5% と言えます。
この最低確率を採用した時に品質水準として 95% とすると、16回の試行で再現しなければ OKと言えます(最初のベイズ推定の数列より)。

数値積分は単純な台形則を利用。最初は真面目に多項式の積分をしていたのですが、n1,n2 が大きいと係数がとても大きくなってしまうので計算結果が不安定になってしいまいます。後で二分探索しているので(F(x)は[0, 1]の範囲で単調増加)最初に適当な精度h を決めて x = 0 から累積していった方が早いかも。次数(試行回数)が増えると結構時間かかります。でも今日はもうここまで。

ちなみにバグ修正の確認だけじゃなくて、他にも2つの状態を取る確率的事象について、確率を推定するのにも使えます。ただし試行毎に独立した確率過程でないといけませんけどね。