概要
陰陽座の楽曲『廿弐匹目は毒蝮』の状況を、Pythonでシミュレートしてみました。
前回の投稿では、「ほんのわずかな確率でも、繰り返していれば起こる」ことをシミュレーションしました。理論上の数値と実際の結果を比較してみるのが楽しかったです。
他のネタでもやってみようと思い、この曲のことを思い出しました。
シチュエーション
「この茂みから100回に1回、もしくは1％の確率で毒蛇が出る」と言われたとします。人によりますが、これは概ね「ではこの1回は大丈夫だろう」と思える確率だと言えるでしょう。しかし、確率1％が持つ本当の意味は“100回やれば必ず1回出る”ではありませんし、逆に“100回やるまでは出ない”という保証もまったくありません。そういうときに限って、中途半端な“22回目”くらいのところで出たりもする…
こんなことを音楽で表現してしまう、瞬火兄上の発想力には敬服します。
再現するにあたって
Pythonでシミュレーションするにあたって、以下のように解釈しています。
- 100本の枝があり、そのうち1本が毒蝮
- ランダムに1本引いて、毒蝮でない場合は、残りの枝の中から再度引く
- 毒蝮に当たるまで引き続ける(99本目まで引けたら完走＝ほんほんえーい)
発生確率
1本目で毒蝮に当たる確率は、瞬火兄上の記載どおり 1/100＝1％ です。
2本目を引く時に毒蝮に当たる確率は、全体が99本になっているので、 1/99＝1.0101…％ にです。
22本目を引く時に毒蝮に当たる確率は、全体が79本になっているので、 1/79＝1.2658…％ です。
（n本目を引く時に毒蝮に当たる確率は、 1/(101-n) と表現できそうです）
22本目を引く時に毒蝮に当たる確率は1.3％にも満たない、と考えると低い確率のように見えますが、これは「21本目まで毒蝮を当たっていない」ことが前提にあります。
100本の中から毒蝮に当たるまで引き続けたときに、22本目で当たる確率は、
「21本目まで毒蝮に当たっていない確率」×「79本の中から毒蝮に当たる確率」になります。
「21本目まで毒蝮を当たっていない確率」をさらに分解すると、
「100本の中から毒蝮に当たらない確率」×「99本の中から毒蝮に当たらない確率」× … × 「80本の中から毒蝮に当たらない確率」です。
(99/100) * (98/99) * … * (79/80) = 79/100
つまり、「21本目まで毒蝮に当たっていない確率」×「79本の中から毒蝮に当たる確率」は、
(79/100) * (1/79) = 1/100
22本目で毒蝮に当たる確率は、1/100＝1％であることが確認できました。
結局、22本目で当たる確率も、1本目でいきなり当たる確率も、99本目まで当たりを引かない確率も、等しく1/100です。
また、「22本目までに毒蝮に当たる確率」は、
1 - (79/100) + (1/79) = 22/100
あるいは、1～22本目それぞれで毒蝮に当たる確率が1/100なので、
(1/100) * 22 = 22/100
と計算できます。2割強と、まあまあ起こりえそうな感じになってきます。
コード Ver.1
0から99のリストからrandom.choiceしていき、0＝毒蝮が出るまでループさせています。
progress関数は、is_progress変数がTrueのときに、途中の実行結果を表示させるために作成しました。
import random
from statistics import mean
def progress(*args):
is_progress = True
if is_progress:
for i in range(len(args)):
print(args[i], end=' ')
print('')
done = 0
results = []
print('～～～ 『廿弐匹目は毒蝮』シミュレータ ～～～')
while True:
try:
turn = int(input('サンプル数を入力してください。>>>'))
break
except:
print('半角数字を入力してください。')
progress('--------')
for t in range(turn):
pick_list = [num for num in range(0, 100)]
for i in range(100):
pick = random.choice(pick_list)
if pick == 0:
progress(i + 1,': 毒蝮！')
results.append(i + 1)
break
else:
pick_list.remove(pick)
if i == 98:
progress('ほんほんえーい！')
done += 1
break
result_22 = results.count(22)
result_1_to_22 = sum(i <= 22 for i in results)
print('--------')
print(f'サンプル数 : {turn}')
print(f'失敗時の平均回数 : {mean(results):.2f}')
print(f'ほんほんえーい : {done} -> {(done/turn):.2%}')
print(f'『廿弐匹目は毒蝮』 : {result_22} -> {(result_22/turn):.2%}')
print(f'『廿弐匹目以内で毒蝮』 : {result_1_to_22} -> {(result_1_to_22/turn):.2%}')
実行結果 Ver1
～～～ 『廿弐匹目は毒蝮』シミュレータ ～～～
サンプル数を入力してください。>>>20
--------
95 : 毒蝮！
89 : 毒蝮！
45 : 毒蝮！
95 : 毒蝮！
3 : 毒蝮！
3 : 毒蝮！
48 : 毒蝮！
68 : 毒蝮！
61 : 毒蝮！
79 : 毒蝮！
22 : 毒蝮！
12 : 毒蝮！
45 : 毒蝮！
66 : 毒蝮！
16 : 毒蝮！
24 : 毒蝮！
27 : 毒蝮！
74 : 毒蝮！
79 : 毒蝮！
84 : 毒蝮！
--------
サンプル数 : 20
失敗時の平均回数 : 51.75
ほんほんえーい : 0 -> 0.00%
『廿弐匹目は毒蝮』 : 1 -> 5.00%
『廿弐匹目以内で毒蝮』 : 5 -> 25.00%
コード Ver.2
リストをシャッフル(全体をsample)して、0＝毒蝮の位置を調べるバージョンです。
シャッフルしたリストの2番目(内部的には1)以降から順に引いていくイメージなので、シャッフル後の先頭が「0」だったら、完走＝ほんほんえーい になります。
ループが1か所減ってるので、choiseより処理は軽いはず。
import random
from statistics import mean
def progress(*args):
is_progress = False
if is_progress:
for i in range(len(args)):
print(args[i], end=' ')
print('')
done = 0
results = []
print('～～～ 『廿弐匹目は毒蝮』シミュレータ ～～～')
while True:
try:
turn = int(input('サンプル数を入力してください。>>>'))
break
except:
print('半角数字を入力してください。')
progress('--------')
pick_list = [num for num in range(0, 100)]
for t in range(turn):
shuffled_list = random.sample(pick_list, len(pick_list))
pick = shuffled_list.index(0)
if pick == 0:
progress('ほんほんえーい！')
done += 1
else:
progress(pick, ': 毒蝮！')
results.append(pick)
result_22 = results.count(22)
result_1_to_22 = sum(i <= 22 for i in results)
print('--------')
print(f'サンプル数 : {turn}')
print(f'失敗時の平均回数 : {mean(results):.2f}')
print(f'ほんほんえーい : {done} -> {(done/turn):.2%}')
print(f'『廿弐匹目は毒蝮』 : {result_22} -> {(result_22/turn):.2%}')
print(f'『廿弐匹目以内で毒蝮』 : {result_1_to_22} -> {(result_1_to_22/turn):.2%}')
実行結果 Ver2 （試行回数 100000回）
十数秒で完了しています。is_progress変数はFalseに変更。
かなり理論値に近づいてますね。
～～～ 『廿弐匹目は毒蝮』シミュレータ ～～～
サンプル数を入力してください。>>>100000
--------
サンプル数 : 100000
失敗時の平均回数 : 49.97
ほんほんえーい : 985 -> 0.98%
『廿弐匹目は毒蝮』 : 1012 -> 1.01%
『廿弐匹目以内で毒蝮』 : 21978 -> 21.98%
所感
- Pythonの勉強以上に、確率の再勉強になりました。
- choiceを使って1本ずつ引いていく直感的なVer.1、「0」の位置を検索するよう発想を切り替えて高速化を図ったVer.2を作成しました。目的に対してどんな手段を使えばよいか、過去事例を調べたり、確率を計算して「つまり何をしてるか」を考えたりすることが重要ですね。
- 途中結果を表示させる関数を作成してみました。デバッグ系のモジュールとして同じものはありそうですが、関数定義の勉強のため。
- 処理時間の確認をまた後回しにしちゃいました。timeitを使うようです。次回の課題。
- 『迦陵頻伽』は名盤。