不均一な確率分布のコンプガチャに必要な平均回数計算プログラム【クーポンコレクター問題】

偏った確率のクーポンコレクター問題のクーポンコンプに必要な平均回数計算プログラムを書きます。コピペして自由に遊んでください。

※この記事ではクーポンをデコ（ピクミンブルームにおけるクーポン）と呼ぶことがあります。

平均回数を求めるシミュレーションプログラム
1. プログラム
2. デコの種類によるpdistなどの書き換え方
参考サイト
その他
別の記法のプログラム

平均回数を求めるシミュレーションプログラム

プログラム

言語はPythonです。pやpdistをお好みの値に変えて（書き換え方は後述）、ご自由に実行してください。

import itertools

#★  確率を記述するためのパラメータ。必須ではない
p = 0.14

#★★  各クーポンの取得確率。取得済みクーポンの要素は書かなくてよい
pdist = [1/7,1/7,1/7,1/7,1/7,1/7,1/7] 

a=0

for n in range(1,len(pdist)+1):
    for team in itertools.combinations(pdist, n):
        a += ((-1)**(len(team)+1))*(sum(team)**-1)
print()
print("平均：", a)

つかいかた

プログラム中の★★部分のpdist（各クーポンの出現確率）を書き換えます。例えば7種なら

pdist = [1/7,1/7,1/7,1/7,1/7,1/7,1/7]。

7種デコで0,1,7番目をコンプ済みなら、先頭から1番目、2番目、7番目のセルを削除して

pdist = [1/7,1/7,1/7,1/7,1/7,1/7,1/7]　→　pdist = [1/7,1/7,1/7,1/7]。

pは必要なら要変更。詳細は後述の「デコの種類によるpdistなどの書き換え方」参照。

出力結果

上のプログラムのpdistは、7種のデココンプまでにかかる苗取得回数の平均を

平均： 18.150000000000006

と出力します。理論値は18.15回で、小数第15位はコンピュータでの計算による誤差です。

デコの種類によるpdistなどの書き換え方

デコ対象のピクミンの種類、デコの出現率などによってpやpdistをどのように書き換えるべきかの例を下に書きます。

pはそのデコの出現率を表します。好きに書き換えて実験してください。

#レアレストラン　[赤,黄,青,白,紫,岩,羽,レア赤,レア黄,レア青]
p=0.15
pdist = [(1-p)/7,(1-p)/7,(1-p)/7,1/7,1/7,1/7,1/7,p/7,p/7,p/7]


#7種デコ　[赤,黄,青,白,紫,岩,羽]
pdist = [1/7,1/7,1/7,1/7,1/7,1/7,1/7]


#3種デコ
pdist = [1/3,1/3,1/3]


#4種デコ
pdist = [1/4,1/4,1/4,1/4]


#14種、デコ1とデコ2の偏りなし　[赤1,黄1,青1,白1,紫1,岩1,羽1,赤2,黄2,青2,白2,紫2,岩2,羽2]
pdist = [1/14]*14
#※ [1/k]*k は、k個のセルがすべて1/k のリスト


#14種、デコの偏りあり（公園など）　[赤1,黄1,青1,白1,紫1,岩1,羽1,赤2,黄2,青2,白2,紫2,岩2,羽2]
pdist = [(1-p)/7,(1-p)/7,(1-p)/7,(1-p)/7,(1-p)/7,(1-p)/7,(1-p)/7,p/7,p/7,p/7,p/7,p/7,p/7,p/7]

pdistのうち、入手済みのデコにあたるセルの部分の数字は削除します。

例えばレアレストランのうち通常7種を入手済みの場合、pdistを

pdist = [(1-p)/7,(1-p)/7,(1-p)/7,1/7,1/7,1/7,1/7,p/7,p/7,p/7]
↓
pdist = [p/7,p/7,p/7]

と書き換えてプログラムを実行します。

参考サイト

・Wikipedia（クーポンコレクター問題）

・Flajolet, Philippe; Gardy, Danièle; Thimonier, Loÿs (1992), “Birthday paradox, coupon collectors, caching algorithms and self-organizing search”

・リスト内のべき集合を取るためのプログラム：

【Python】リスト内の要素の全ての組み合わせを出力

その他

m種類デコのうち未収集デコ k 種類のときのコンプ回数倍率

同じ確率の未収集デコ n 種類をコンプする期待値は、デコ出現確率 p の逆数 1/p に H(n)（1～nの逆数の和：1/1 + 1/2 + 1/3 + ⋯ + 1/n）をかければOK。

例えば7種類等確率なら、出現確率1/7の逆数7に、1～k（kは残りのデコの数）の逆数の和をかければ、未収集が0,1,2,3,4,5,6,7のときのコンプまでの期待値はそれぞれ

0, 7, 10.5, 12.83, 14.58, 15.98, 17.15, 18.15

プログラム実行にかかる時間

O(2^len(pdist)) 　（len(pdist)：pdistの長さ（配列サイズ）、つまり収集したいピクミンの数）

長さ20までなら1秒以内に結果が出るものの、長さ25で20秒。そこから長さが1増えるたびに2倍の実行時間がかかってしまう

プログラムの計算原理について

・不均一な確率分布のクーポンコレクター問題の期待値：

$E(T)=\int _{0}^{\infty }{\big (}1-\prod _{i=1}^{n}(1-e^{-p_{i}t}){\big )}dt$

出典：Wikipedia（クーポンコレクター問題）

展開して積分すれば、不均一な確率分布のリストを構成してそのPower Set（べき集合）のうち空集合を除くもののうち、それぞれに確率の和を求める→逆数を求める→集合の要素数が奇数なら＋、偶数ならーをかけて全部足す、というのをシミュレーションプログラムで記述しました。

※デコ入手済みのものがある場合、それに対応するセルの値を削除してしまってそのままプログラムを実行しても正しい値が出ているような気がするが、本当に正しいのかは未証明

別の記法のプログラム

出力内容は最初のものと同じです。

引数：p, pdist, want

pdistはすべての種類のデコ出現確率を記載。例えばレア含めたレストランなら7種＋レア3種の順に[(1-p)/7,(1-p)/7,(1-p)/7,1/7,1/7,1/7,1/7,p/7,p/7,p/7]

配列wantは欲しいデコピクミンの種類を入力。want = [0,0,0,0,0,0,0,0,0,1]なら、レアレストランの青だけゲットするまでの回数の期待値。この配列の長さはpdistの長さと同じであること。

import itertools

p=0.14 #レアデコになる確率
pdist = [(1-p)/7,(1-p)/7,(1-p)/7,1/7,1/7,1/7,1/7,p/7,p/7,p/7] #各色ピクミンのデコ苗取得確率
want = [0,0,0,0,0,0,0,0,0,1] #欲しいデコピクミンの種類。1なら欲しい、0なら欲しくない

a=0
wantpik=[]

for j in range(0,len(pdist)):
    if want[j] != 0:
        wantpik.append(pdist[j])

for n in range(1,len(pdist)+1):
    for team in itertools.combinations(wantpik, n):
        a += ((-1)**(len(team)+1))*(sum(team)**-1)
print()
print("平均：", a)