二項分布・ポアソン分布・指数分布
回数と時間に合った確率モデルを選び、固定seedの試行で理論値と推定値を比較します。
説明
二項分布は、独立な一定回数の試行で成功する回数を表します。ポアソン分布は一定時間内の到着回数などを表し、一定率のポアソン過程では到着までの待ち時間が指数分布になります。回数は整数、待ち時間は連続量です。
成功確率0.3の試行を10回繰り返す例と、1秒当たり平均2件の到着モデルを扱います。確率式と、seedを固定した乱数による試行を比較します。最後に試行回数を増やしたときの平均の揺れを調べます。
1成功回数の確率
10回中k回成功する確率は、成功する位置の組合せ数×0.3のk乗×0.7の残り回数乗です。
import math, random確率式と乱数; n = 10; p = 0.3回数と成功確率; binomial = [math.comb(n,k)*p**k*(1-p)**(n-k) for k in range(n+1)]二項分布の確率; [(k, round(prob,4)) for k,prob in enumerate(binomial)]成功回数と確率成功回数0から10までの分布です。小数第4位に丸めているため、10回成功の表示は0.0ですが、確率そのものは0ではありません。各回の成功確率が一定で、試行同士が独立という条件を使っています。成功位置の組合せを数えるためcombが必要です。
2理論上の平均と分散
確率を足して1になるか、値×確率から平均が3になるかを確認します。二項分布の分散はnp(1-p)です。
round(sum(binomial),12)確率の合計, sum(k*prob for k,prob in enumerate(binomial))分布から期待値, n*p, n*p*(1-p)平均と分散の公式平均は3、分散は2.1です。浮動小数点の計算では末尾に小さな誤差が出るため、確率合計の表示を丸めています。
3同じ試行を1000組実施
一様乱数が0.3未満なら成功とし、10回分の成功を数えます。その組を1000回作ります。
rng = random.Random(2026)独立した乱数器と固定seed; trials = [sum(rng.random() < p for _ in range(n)) for _ in range(1000)]10回試行を1000組; from collections import Counter; counts = Counter(trials)成功回数を集計; sorted(counts.items()), sum(trials)/len(trials)度数と平均各成功回数の度数と平均が表示されます。seedを固定すると同じ環境で比較しやすくなりますが、一回の試行列が理論分布そのものになるわけではありません。
4一定時間内の到着回数
平均到着数λ=2のポアソン分布では、k件の確率はexp(-λ)λ^k/k!です。0から8件までを表示します。
rate = 21秒当たり平均2件; poisson = [math.exp(-rate)*rate**k/math.factorial(k) for k in range(9)]0〜8件の確率; [(k, round(prob,4)) for k,prob in enumerate(poisson)]回数と確率, 1-sum(poisson)表示範囲外の確率最後の値は9件以上の確率です。最大8件というモデルではなく、表示範囲の外にも確率が残っています。到着率が一定で独立な増分を持つなど、ポアソン過程の仮定が合う場合に使います。
5次の到着までの待ち時間
同じ率2のポアソン過程では、待ち時間の平均は1/2秒です。指数分布の乱数を1000個生成します。
arrival_rng = random.Random(2026)待ち時間用の乱数器; waits = [arrival_rng.expovariate(rate) for _ in range(1000)]指数分布の待ち時間; [round(x,4) for x in waits[:8]]最初の8間隔, round(sum(waits)/len(waits),4), 1/rate標本平均と理論平均待ち時間は整数回数ではなく秒数です。最初の数個と平均を表示しています。平均0.5秒でも毎回0.5秒おきになるわけではなく、短い間隔も長い間隔もあります。
6待ち時間から1秒内の件数を作る
待ち時間を足し、累積時刻が1秒を超えるまでに何件到着したかを数えます。
def arrivals(rng, rate, duration):乱数器・率・観測時間
time, count = 0.0, 0累積時刻と件数
while True:到着を追加
time += rng.expovariate(rate)次の間隔を加える
if time > duration: return count観測窓の外で終了
count += 1窓内の到着を数える1秒を超えた到着は件数に入れません。時間を区切って件数を数えることで、指数分布の間隔とポアソン分布の到着数を同じ過程として結び付けます。
7到着回数の分布を比較
1秒間の観測を1000回作り、各件数の度数と平均を求めます。
count_rng = random.Random(2026)件数用の乱数器; observed = [arrivals(count_rng, rate, 1.0) for _ in range(1000)]1秒の窓を1000回; sorted(Counter(observed).items()), sum(observed)/len(observed)度数と平均到着数平均到着数は理論値2の近くになります。0件の窓も複数件の窓もあり、平均だけでは各窓の値は決まりません。有限試行の度数は理論確率と完全には一致しません。
8試行回数と推定の揺れ
先ほど作った1000組の二項試行を、先頭10・100・1000組まで使って平均を推定します。
[(size, sum(trials[:size])/size, round(math.sqrt(n*p*(1-p)/size),4)) for size in [10, 100, 1000]]試行数・推定平均・理論上の揺れ最後の列は独立な標本平均の標準偏差√(2.1/試行数)です。試行数を100倍にすると、この幅は10分の1になります。ただし一つの実現列の誤差が、試行を増やすたびに必ず単調に小さくなるわけではありません。
まとめ
二項分布は一定回数の独立試行の成功回数、ポアソン分布は一定区間の到着回数を表します。
一定率のポアソン過程では到着間隔が指数分布となり、率2なら平均間隔0.5です。
試行による推定は揺れます。seedで再現できることと、理論値に一致することは別です。
操作履歴と解説を振り返れます。
Python 3
import math, random; n = 10; p = 0.3; binomial = [math.comb(n,k)*p**k*(1-p)**(n-k) for k in range(n+1)]; [(k, round(prob,4)) for k,prob in enumerate(binomial)][(0, 0.0282), (1, 0.1211), (2, 0.2335), (3, 0.2668), (4, 0.2001), (5, 0.1029), (6, 0.0368), (7, 0.009), (8, 0.0014), (9, 0.0001), (10, 0.0)]
round(sum(binomial),12), sum(k*prob for k,prob in enumerate(binomial)), n*p, n*p*(1-p)(1.0, 2.9999999999999987, 3.0, 2.0999999999999996)
rng = random.Random(2026); trials = [sum(rng.random() < p for _ in range(n)) for _ in range(1000)]; from collections import Counter; counts = Counter(trials); sorted(counts.items()), sum(trials)/len(trials)([(0, 32), (1, 123), (2, 219), (3, 265), (4, 198), (5, 110), (6, 40), (7, 10), (8, 2), (9, 1)], 3.033)
rate = 2; poisson = [math.exp(-rate)*rate**k/math.factorial(k) for k in range(9)]; [(k, round(prob,4)) for k,prob in enumerate(poisson)], 1-sum(poisson)([(0, 0.1353), (1, 0.2707), (2, 0.2707), (3, 0.1804), (4, 0.0902), (5, 0.0361), (6, 0.012), (7, 0.0034), (8, 0.0009)], 0.0002374473282610312)
arrival_rng = random.Random(2026); waits = [arrival_rng.expovariate(rate) for _ in range(1000)]; [round(x,4) for x in waits[:8]], round(sum(waits)/len(waits),4), 1/rate([0.0634, 0.3491, 0.3585, 0.9831, 0.0541, 0.1263, 0.4594, 0.4066], 0.4945, 0.5)
def arrivals(rng, rate, duration):
time, count = 0.0, 0
while True:
time += rng.expovariate(rate)
if time > duration: return count
count += 1count_rng = random.Random(2026); observed = [arrivals(count_rng, rate, 1.0) for _ in range(1000)]; sorted(Counter(observed).items()), sum(observed)/len(observed)([(0, 139), (1, 247), (2, 283), (3, 168), (4, 98), (5, 46), (6, 14), (7, 4), (8, 1)], 2.059)
[(size, sum(trials[:size])/size, round(math.sqrt(n*p*(1-p)/size),4)) for size in [10, 100, 1000]][(10, 2.2, 0.4583), (100, 2.86, 0.1449), (1000, 3.033, 0.0458)]
演習完了