ポアソン分布のグラフ表示がうまく表示されません

Question

ポアソン分布のグラフを実装していたところ平均値であるλが10の時だけグラフの点がうまく表示されません。具体的には下の二つの図が違うという事です。(最後の方が違う)
内部的な点の値はあっているのですが点がプロットされない原因が分かる方、教えて下さると光栄です。プログラムも載せておきます。

これが実行結果です。

また次のグラフが本当のグラフです。

ここからはコードを載せます。

#λは平均, kはその時の値。
#つまりポアソン分布は平均λの場合にkの値が出る確率
def poason(k, λ):
    if type(k)==int:
        sum=1
        for j in range(1,k+1):
            sum = sum*j
        return np.e**(-λ)*(λ**k)/sum
    else:
        a = np.zeros(len(k))         #aを作らずにkに代入しても良いがkはnp.arangeで作っているため整数しか入らない.　そのためaを作成
        for i in range(len(k)):
            sum = 1
            for j in range(1,k[i]+1):
                sum = sum*j
            a[i]=np.e**(-1*λ)*(λ**k[i])/sum
        return a


import matplotlib.pyplot as plt
import numpy as np 

fig, ax = plt.subplots(figsize=(10,6))

k = np.arange(21)
p1 = poason(k, 1)  #λが1の場合のグラフ
p2 = poason(k, 4)  #λが4の場合のグラフ
p3 = poason(k,10)  #λが10の場合のグラフ

c1, c2, c3 = "blue", "green", "red"   #各プロットの色
l1, l2, l3 = "λ=1", "λ=4", "λ=10"  #各ラベル

ax.set_xlabel("k", fontsize=20)     #x軸ラベル
ax.set_ylabel("p(x=k)", fontsize=20)          #y軸ラベル
ax.set_title("Poisson distribution", fontsize=20)   #グラフタイトル

ax.grid()    #罫線(けいせん)。縦と横の細い線

ax.plot(k, p1, marker = "o", color = c1, label = l1)
ax.plot(k, p2, marker = "o", color = c2, label = l2)
ax.plot(k, p3, marker = "o", color = c3, label = l3)

#xの目盛りを作る(まあ作らなくても良いがfigsizeなど指定するには必要)目盛りの大きさのみ変えるだけなら一行で書ける。yは一行。
i=0
xtick = np.zeros(1)
while (i<21):
    if i%5==0:
        xtick = np.append(xtick,i)
    i +=1

ax.set_xticks(xtick)  #これでx軸の目盛りの個数と大きさは分かった。z残念ながらここにfontizeはない
ax.set_xticklabels(xtick, fontsize=15)    #第一引数はラベルであるので"A"などの名前でも良い。ここでfontsize指定できる

ax.tick_params(axis='y',labelsize=15)   #このように目盛りの大きさだけならこれで変えられる

ax.legend(loc=0, fontsize=20)   #凡例表示。線の説明。locは場所。0の場合は右上
fig.tight_layout()

plt.show()

Answer 1

32ビット整数を使って計算しているためです。
poason関数の定義を以下に変更すれば、期待している結果になります。

def poason(k, λ):
    if type(k)==int:
        sum=1
        for j in range(1,k+1):
            sum = sum*j
        return np.e**(-λ)*(λ**k)/sum
    else:
        a = np.zeros(len(k))         #aを作らずにkに代入しても良いがkはnp.arangeで作っているため整数しか入らない.　そのためaを作成
        for i in range(len(k)):
            sum = 1
            for j in range(1,k[i]+1):
                sum = sum*j
            a[i]=np.e**(-1*λ)*(λ**int(k[i]))/sum
        return a

Comments

[kaggle777]

丁寧な解説ありがとうございます。たしかに32bit整数のまま計算すると小数点が多い場合に格納できないですね。とても勉強になりました。ありがとうございました。