クーポンコレクター問題の理論的な解やグラフ、分布関数出力プログラムや推移行列など

from sympy.functions.combinatorial.numbers import stirling
import matplotlib.pyplot as plt
import math
import numpy as np


k = CouponKind = 7  #クーポンの数
numTrial = 50  #入手回数の上限値


x = []
y = []
for n in range(1, numTrial+1):
    stir = stirling(n-1, k-1, kind = 2)
    prob = math.factorial(k) * stir / (k**n)
    x.append(n)
    y.append(prob)
plt.bar(x, y, label='k = {}'.format(k))
plt.xticks(np.arange(0, numTrial+1, step=10))  #stepは適宜変更
plt.grid(linestyle="--")
plt.legend()
plt.show()

z=[0]*numTrial
psum=0
yy=[0]+y
for i in range(0,numTrial):
    psum += y[i]
    z[i] = psum
plt.plot(x, z, label='k = {} cdf'.format(k),)
plt.xticks(np.arange(0, numTrial+1, step=10))  #stepは適宜変更
plt.yticks(np.arange(0, 1.1, step=0.1))
plt.grid(linestyle="--")
plt.legend()
plt.show()
print()


hk=0
for i in range(1,CouponKind+1):
    hk += 1/i
var=0
for j in range(1,CouponKind+1):
    var += (1-j/CouponKind)/((j/CouponKind)**2)
print("●",CouponKind,"種類")
print("平均：",CouponKind*hk,"回")
print("分散：",var)
print("標準偏差：",math.sqrt(var))
print()

zz=[0]*numTrial
for i in range(0,numTrial):
    zz[i]=float(z[i])
    print(i+1,"回目まで：",100*zz[i],"%")

from sympy.functions.combinatorial.numbers import stirling
import matplotlib.pyplot as plt
import math
import numpy as np


CouponKind = [3,4,7,14]  #デコピクミンの種類数。[]の中にコンマ付きで入力
numTrial = 70  #苗の入手回数の上限値


for k in CouponKind:
  x = []
  y = []
  for n in range(1, numTrial+1):
    stir = stirling(n-1, k-1, kind = 2)
    prob = math.factorial(k) * stir / (k**n)
    x.append(n)
    y.append(prob)
  plt.plot(x, y, label='k = {}'.format(k))
  plt.xticks(np.arange(0, numTrial+1, step=5))  #stepは適宜変更
  plt.grid(linestyle="--")
plt.legend()
plt.show()


for k in CouponKind:
  x = []
  y = []
  for n in range(1, numTrial+1):
    stir = stirling(n-1, k-1, kind = 2)
    prob = math.factorial(k) * stir / (k**n)
    x.append(n)
    y.append(prob)
  z=[0]*numTrial
  psum=0
  yy=[0]+y
  for i in range(0,numTrial):
    psum += y[i]
    z[i] = psum
  plt.plot(x, z, label='k = {} cdf'.format(k),)
  plt.xticks(np.arange(0, numTrial+1, step=5))  #stepは適宜変更
  plt.yticks(np.arange(0, 1.1, step=0.1))
  plt.grid(linestyle="--")
plt.legend()
plt.show()

import matplotlib.pyplot as plt
import numpy as np
import math


Coupon = 24 #デコ数
have = 0 #すでに持っているデコ数。ただしk + have ≦ Couponであること。超えると必ずダブる；Daburanai = 0
k = Coupon - have #何個の苗を取るか

def Daburanai(C,get,havekind):
    return math.factorial(C-havekind)/math.factorial(C-get-havekind)/C**get

def Daburi(C,get,havekind):
    return 1-Daburanai(C,get,havekind)


z=[0]*(k+1)
left=[x for x in range(0,k+1)]
for i in range(0,k+1):
    #print(i,"個：",100-100*Daburanai(365,i,have),"%")
    z[i] = Daburanai(Coupon,i,have) #★　ダブり確率。適宜ダブらない"Daburanai(Coupon,i,have)"に変更可


plt.xticks(np.arange(0, k+1, step=2)) #stepはx軸の目盛り幅。適宜変更
plt.yticks(np.arange(0, 1.1, step=0.1))
plt.xlabel('Number of Coupon')
plt.ylabel('Probability')
plt.grid(linestyle="--")
plt.plot(left, z)
plt.show()

for i in range(1,k+1):
    print(i,"個：",100*z[i],"%")