概要

陰陽座の楽曲『廿弐匹目は毒蝮』の状況を、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を使うようです。次回の課題。
『迦陵頻伽』は名盤。