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')