【ピクミンブルーム】デコピクミンコンプリートまでの苗取得回数計算プログラム【不均一な確率、複数デコゲット対応】

import numpy as np
import collections
import matplotlib.pyplot as plt
import random

pdist = [1/7]*7  #ピクミン入手確率の偏り
want = [2]*7  #欲しいピクミンの種類と数
lim=80  #苗の入手回数の上限値
NN=N=100000  #試行回数

Type=len(pdist)
hantei = [0]*Type
Pik=[0]*Type

for m in range(0,Type):
    Pik[m] = -want[m]
cumpdist=[0]*(Type+1)
count=0
sum=cump=0
c=[0]*(lim+1)
for i in range (0,Type):
    cump += pdist[i]
    cumpdist[i] = cump
while N>0:
    for j in range(1,lim+1):
        for k in range(0,Type):
            if cumpdist[k-1] <= random.random() < cumpdist[k]:
                Pik[k] += 1
                count=count+1
        if min(Pik) >= 0:
            sum=sum+count
            if count <= lim:
                c[count] += 1
            if N < 101:
                print(count,"回、", [Pik[x]+want[x] for x in range(0,Type)])
            break
    for m in range(0,Type):
        Pik[m] = -want[m]
    count=0
    N=N-1

print()
print("●",Type,"種類")
print("平均",sum/NN,"回")
Xaxis=list(range(0,lim+1))
plt.xticks(np.arange(0, lim+1, step=10))  #stepは適宜変更
plt.grid(linestyle="--")
plt.bar(Xaxis,[x/NN for x in c])
plt.show()

y=[x/NN for x in c]
z=[0]*(lim+1)
psum=0
for i in range(0,lim+1):
    psum += y[i]
    z[i] = psum
plt.plot(Xaxis, z)
plt.xticks(np.arange(0, lim+1, step=10))  #stepは適宜変更
plt.yticks(np.arange(0, 1.1, step=0.1))
plt.grid(linestyle="--")
plt.show()

import numpy as np
import collections
import matplotlib.pyplot as plt
import random


def get_days(k):
    fossil_list = [0]*k #最初のクーポン所持数
    days = 0
    while min(fossil_list) < 2: #クーポンが全部x個以上になったら完了
        num_fossil_per_day = 1
        for i in range(num_fossil_per_day):
            fossil_list[random.randint(0, k - 1)] += 1
        days += 1
    return days

# かかった日数をnum回取得して分布を取得する。numが大きいほど分布は正確になる
def days_plot(k=73, num=10000):
    days = [get_days(k) for i in range(num)]

    c = collections.Counter(days)
    key = list(c.keys())
    value = list(c.values())
    left = np.array(key)
    height = np.array(value) / sampling_num
    plt.title('number of samples: {}'.format(num))
    plt.xlabel('number of days for complete')
    plt.ylabel('probability')
    plt.bar(left, height)
    plt.show()
    print("平均：", sum(left*height))


# クーポンの種類数
k = 7
# データサンプル数
sampling_num = 100000

# 横軸が日数、縦軸が N日目に初めてコンプリートする確率分布
days_plot(k, sampling_num)