Rayleigh分布のあれ

きむらさんがすっきりしないようなので、シミュレーションしてみた。

  • スケールパラメータ 5*sqrt(2) *1で、Layleigh分布*2 の乱数 t を496回生成。
  • n - 1 <= t < n の場合に、n週目発見の欠陥としてカウント
  • 各週の発見欠陥数から K_i を計算

ということを1,000回実施して、K_i の平均と標準偏差を計算。

import random
import math

NUM_TRIAL = 1000
TOTAL_DEFECTS = 496
ESTIMATED_TM = 5
T_MAX = 10

def mean(seq):
    return sum(seq) / len(seq)

def std(seq):
    m = mean(seq)
    return math.sqrt(sum((x - m) ** 2 for x in seq) / len(seq))

def trial(total_defects, e_tm, t_max):
    counts = [0 for _ in xrange(t_max)]
    scale = e_tm * math.sqrt(2)
    for _ in xrange(total_defects):
        t = random.weibullvariate(scale, 2)
        idx = int(t)
        if idx < t_max:
            counts[idx] += 1
    k_estimates = []
    for i, c in enumerate(counts, 1):
        k = 1.0 * e_tm * e_tm / i * math.exp(0.5 * i * i / e_tm / e_tm) * c
        k_estimates.append(k)
    return { 'counts': counts, 'k_estimates': k_estimates }

def main():
    k_all = [[] for _ in xrange(T_MAX)]

    for _ in xrange(NUM_TRIAL):
        result = trial(TOTAL_DEFECTS, ESTIMATED_TM, T_MAX)
        for i, c in enumerate(result['k_estimates']):
            k_all[i].append(c)

    for i, k in enumerate(k_all, 1):
        print i, mean(k), std(k)

if __name__ == '__main__':
    main()

実行結果

1 251.352105149 79.4676230281
2 376.753701049 69.2076414754
3 432.474818805 62.0288908122
4 470.504307734 58.1793029766
5 490.379167544 60.2312676994
6 507.273800261 61.3508169168
7 526.239623696 68.082631596
8 543.75572901 72.5164253587
9 552.920371528 84.6362785681
10 576.050813473 95.2159227826

なんだが、週によって全然違う分布を示すようなので、スライドのように K_i の平均や標準偏差計算しても意味ないんじゃ無いですかね。

*1:こうすれば tm=5 になる

*2:つまり形状パラメータ2の Weibull 分布