PythonでPoisson分布を扱う便利なライブラリとしては，scipy.stats.poissonがある．
docs.scipy.org
早速使ってみる．
まずは，pmf関数を用いて確率密度を計算する．
ポアソン分布は平均回数muを指定すると，muを期待値として，もっともらしい，発生回数，訪問者数や購入数の分布を作ってくれる．
以下は，muを10にした場合の0～50までの確率密度を計算し，プロットした例である．
from scipy.stats import poisson import matplotlib.pylab as plt import numpy as np mu = 5 xs = np.arange(50) pdfs = [poisson.pmf(x,mu) for x in xs] plt.bar(xs, pdfs) plt.xlabel('x') plt.ylabel('density')
ポアソンは，回数muが大きくなるほど正規分布に近づいていくが，mu=10の場合は，若干右の方（x=20あたり）が密度が高くなっているのがわかる．もう少し非対称性を確認するために，mu=5のときを試してみる。
より右の方が裾野が広がっているのがわかる。
次に，確率（確率密度関数の面積）を求める関数cdfを用いて5%点と95%点を近似的に求めてみる．
ちなみに、poissonクラスにどのようなメソッドがあるかは，以下のようおにdictで調べることができる．
> dir(poisson) ['__call__', '__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__getstate__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__setstate__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_argcheck', '_argcheck_rvs', '_cdf', '_cdf_single', '_cdfvec', '_construct_argparser', '_construct_default_doc', '_construct_doc', '_construct_docstrings', '_ctor_param', '_entropy', '_isf', '_logcdf', '_logpmf', '_logsf', '_munp', '_nonzero', '_open_support_mask', '_parse_args', '_parse_args_rvs', '_parse_args_stats', '_pmf', '_ppf', '_ppfvec', '_random_state', '_rvs', '_sf', '_stats', '_stats_has_moments', '_support_mask', '_updated_ctor_param', 'a', 'b', 'badvalue', 'cdf', 'entropy', 'expect', 'extradoc', 'freeze', 'generic_moment', 'inc', 'interval', 'isf', 'logcdf', 'logpmf', 'logsf', 'mean', 'median', 'moment', 'moment_tol', 'name', 'numargs', 'pmf', 'ppf', 'random_state', 'rvs', 'sf', 'shapes', 'stats', 'std', 'var', 'vecentropy']
本題に戻って，以下のコードでパーセント点を近似的に求めてみた．
以下，mu=5，信頼区間としてp=95%を設定した場合の
mu=5 p = 0.95 #累積確率の計算 cmfs=np.array([poisson.cdf(x,mu) for x in xs]) # パーセント点の計算 lowerX = np.max(xs[cmfs<(1-p)/2]) higherX = np.min(xs[cmfs>(1+p)/2]) print(f"{(1-p)/2:.2f}点:{lowerX}，{(1+p)/2:.2f}点:{higherX}") #プロット pdfs = [poisson.pmf(x,mu) for x in xs] plt.bar(xs, pdfs) plt.bar(xs[:lowerX], pdfs[:lowerX],color="red") plt.bar(xs[higherX:], pdfs[higherX:],color="red") plt.xlabel('x') plt.ylabel('density')
