ポアソン分布 Poisson distribution

タイトルでポアソン分布、とか言いつつ別にポアソン分布の性質とか何かに深入りしようというでもなく、ちょっと確率分布のグラフを書きたいだけのエントリです
このブログの想定読者はポアソン分布は当然知ってるのでRSSリーダーでもタイトルだけ見て既読扱いにしてるはずです

以下は単なるPython+numpy+scipy+matplotlibコーディング練習場兼コピペ材料置き場です

ポアソン分布 Poisson distribution

定数 $\lambda>0$ に対し、自然数を値にとる確率変数Xが
　　 $P(X=k)=\frac{\lambda^ke^{-\lambda}}{k!}$
を満たすとき、確率変数Xはパラメータλのポアソン分布に従うという。

ここで、e はネイピア数 (e = 2.71828...)であり、k! は k の階乗を表す。また、λ は所与の区間内で発生する事象の期待発生回数に等しい。

P(X=k)は、「単位時間中に平均で λ 回発生する事象がちょうど k 回（k は0を含む自然数、k = 0, 1, 2, ...）発生する確率」に相当する。例えば、事象が平均で2分間に1回発生する場合、10分間の中で事象が発生する回数は、λ = 5 のポアソン分布モデルを使って求められる。

1時間に特定の交差点を通過する車両の台数。

1mlの希釈された水試料中に含まれる特定の細菌の数（細菌数検査における最確法）。

1ページの文章を入力するとき、綴りを間違える回数。

1日に受け取る電子メールの件数。

1分間のWebサーバへのアクセス数。

例えば、1時間あたりのウィキペディアの最近更新したページの編集数もおおよそポアソン分布。

1マイルあたりのある通り沿いのレストランの軒数。

1ヘクタールあたりのエゾマツの本数。

1立方光年あたりの恒星の個数。

ポアソン分布 - Wikipedia

あ、あとポアソンの極限定理について

n が大きく p が十分小さい場合、np は適度な大きさとなるため、パラメータ λ = np であるポアソン分布が二項分布B(n, p) の良好な近似を与える。すなわち、期待値λ = npを一定とし、nを十分大きくしたとき、
　　 $P(X=k)\simeq\frac{\lambda^ke^{-\lambda}}{k!}$
が成り立つ。この結果は数学者シメオン・ドニ・ポアソンが1837年に著書 Recherches sur la probabilite des jugements (Researches on the Probabilities) の中で与えており、ポアソンの極限定理と呼ばれる。
二項分布 - wikipedia

1. 確率関数の実装

1.1 自分で実装してみる

// パラメータ名にλってなんかあれなのでみんなに合わせてμ(mu)を使うよ

Pythonだと

def pmf(x, mu):
  return mu**x / (e**mu * factorial(x))

みたいな感じか。
ちょっといじってみて、factorialはmath.factorial(n)は引数に配列を渡すと整数値じゃないからだーめとか云われるので、scipy.misc.factorial(n) が便利なことが分かった（そういうのユニバーサル関数って言うんだった？）

パラメータmuを先に渡してカリー化というか部分適用というかするなら

import functools

def pmf2(mu, x):
  return mu**x / (e**mu * factorial(x))

def pmf1_gen(mu):
  return functools.partial(pmf2, mu)

みたいな、あるいはちゃんとe**muを明示的に先に計算する形で

def pmf1_gen(mu):
  e_mu = e ** mu
  return lambda x:mu**x / (e_mu * factorial(x))

とか。（これならグラフ描画時にも e^μ の計算は一度で済む）

まあいい。

from pylab import *
import scipy.misc

def poisson_pmf_gen(mu):
  e_mu = e ** mu
  return lambda x:mu**x / (e_mu * scipy.misc.factorial(x))

for mu in arange(1, 11):
  pmf = poisson_pmf_gen(mu)

  axis([0,15, 0,pmf(mu)+0.05])
  xlabel('x')
  ylabel('p(x)')
  title("Poisson distribution (mu = %d)" % mu)

  x = arange(0,15)
  bar(x, pmf(x), 0.8, color="gray")

  savefig("poisson_naive_%02d" % mu)
  clf()

f:id:n4_t:20120102170411p:plain
f:id:n4_t:20120102170417p:plain
f:id:n4_t:20120102170423p:plain
f:id:n4_t:20120102170428p:plain
f:id:n4_t:20120102170433p:plain
f:id:n4_t:20120102170439p:plain
f:id:n4_t:20120102170445p:plain
f:id:n4_t:20120102170452p:plain
f:id:n4_t:20120102170500p:plain
f:id:n4_t:20120102170505p:plain

1.2 scipy.stats.poisson の pmf() を使ってみる

from pylab import *
from scipy.stats import poisson

x = arange(0,16)

axis([0,15, 0,0.16])
bar(x, poisson.pmf(x,7), 0.8, color="gray")

title("Poisson distribution (mu = 7)")
xlabel('x')
ylabel('p(x)')

show()

f:id:n4_t:20120102170842p:plain

2. ポアソン分布に従う乱数

(0, 1) の一様乱数 $u_1, u_2, u_3, ...$ を用い、 $\displaystyle\prod_{\small i=1}^{\small k}{u_i}\lt e^{\tiny -\lambda}$ となる最初のkを求める方式。→see http://hata.cc/docs/sim/poisson.html http://www.ishikawa-lab.com/montecarlo/4shou.html
(0, 1) の一様乱数 u に対し、 $\displaystyle\sum_{\small k=0}^{\small x-1}{e^{\tiny -\lambda}\frac{m^k}{k!}\leq u\lt \sum_{\small k=0}^{\small x}{e^{\tiny -\lambda}\frac{m^k}{k!}$ となる最小のx(>0)を求める方式。→see http://www.sci.kagoshima-u.ac.jp/~ebsa/wakimoto01/pdf/ch04-02.pdf
numpy.random.poisson() で
scipy.stats.poisson.rvs() で

それぞれμ=7で10万回ぐらいやってヒストグラムを取ってみたりスピード比べしてみたり（※４つのヒストグラムを1画面にまとめて表示したいだけ）

from pylab import *
import numpy.random
import scipy.stats
import timeit

def poisson_rand1(mu):
  th = e**(-mu) ## <1
  x = numpy.random.random()
  k = 0
  while x >= th:
    x *= numpy.random.random()
    k += 1
  return k

def poisson_rand2(mu):
  k = 0.0
  a = 1.0
  s = 0.0
  u = numpy.random.random() * e**mu
  while True:
    s += a
    if s > u:
      return k
    k += 1
    a = a * mu / k

N = 100000
mu = 6

s = []
t = []
time = []
time2 = []

s.append( [poisson_rand1(mu) for x in xrange(N)] )
t.append("poisson_rand1 (mu=%d)" % mu)
tm = timeit.Timer("poisson_rand1(%d)" % mu, "from __main__ import poisson_rand1")
time.append(tm.timeit(N))
time2.append(-1)

s.append( [poisson_rand2(mu) for x in xrange(N)] )
t.append("poisson_rand2 (mu=%d)" % mu)
tm = timeit.Timer("poisson_rand2(%d)" % mu, "from __main__ import poisson_rand2")
time.append(tm.timeit(N))
time2.append(-1)

s.append( numpy.random.poisson(mu, N) )
t.append("numpy.random.poisson(mu=%d)" % mu)
tm = timeit.Timer("numpy.random.poisson(%d)" % (mu), "import numpy.random")
time.append(tm.timeit(N))
tm = timeit.Timer("numpy.random.poisson(%d,%d)" % (mu,N), "import numpy.random")
time2.append(tm.timeit(1))

s.append( scipy.stats.poisson.rvs(mu, size=N) )
t.append("scipy.stats.poisson.rvs(mu=%d)" % mu)
tm = timeit.Timer("scipy.stats.poisson.rvs(%d)" % (mu), "import scipy.stats")
time.append(tm.timeit(N))
tm = timeit.Timer("scipy.stats.poisson.rvs(%d, size=%d)" % (mu,N), "import scipy.stats")
time2.append(tm.timeit(1))

for i in xrange(4):
  subplot(2,2,1+i)
  axis([0,20,0,0.18])
  count, bins, ignored = plt.hist(s[i], arange(21), normed=True)
  text(19, 0.16, "%3.3fmsec" % (time[i]*1000), color='gray', fontsize=12, horizontalalignment='right')
  if time2[i] > 0:
    text(19, 0.15, "%3.3fmsec" % (time2[i]*1000), color='gray', fontsize=12, horizontalalignment='right')
  title(t[i])

show()

f:id:n4_t:20120102234328p:plain
右上の数字はtimeitで計測した実行時間。2段になっているのは、1個ずつN回計算したものと、sizeパラメータにNを渡してまとめて作ってもらった場合のそれぞれの実行時間。まとめて作るならライブラリ関数が断然お得！