matplotlibでポワソン分布のヒストグラムを描画

　最近matplotlibというPythonのグラフ描画ライブラリをインストールしたので、ちょっと遊んでみた。まずは簡単なものからということで、ポワソン分布のヒストグラムを描画してみた。

1. numpy.random.poisson
2. 自作のMCMC（メトロポリス法）

の2通りの方法で乱数を発生させてヒストグラムにしてみた。

自作MCMCはランダムウォークするときの次の数の選び方がかなりテキトーだったのだが、それっぽい結果になっていた。メトロポリス法はダメで、メトロポリス・ヘイスティングを使わないといけないというシナリオを期待していたのだけど...

グラフ描画の処理

公式ドキュメントのヒストグラム作成サンプルを参考に書いた。
rpoisson(lam, n)がポワソン分布に従う乱数を生成する関数。
hist()のビンの数、幅はうまく整数値がはまるように調整した。

import math
import matplotlib.pyplot as plt
from metropolis import rpoisson

# pdf of poisson distibution (mean = lam) at k
def p(lam, k):
    ret = math.exp(-lam)
    for i in xrange(k):
        ret *= 1.0 * lam / (i + 1)
    return ret

n = 1000000  # number of samples
lam = 10      # mean 
samples = rpoisson(lam, n)  # get n samples

maxv = max(samples)
plt.hist(samples, range=(-0.5, maxv+0.5), bins=maxv+1, normed=1, facecolor='blue', alpha=0.5)

y = [p(lam, x) for x in range(maxv+1)]
plt.plot(range(maxv+1), y, 'r')
plt.xlabel('x')
plt.ylabel('p(x)')
plt.title(r'Histogram of x ~ Poisson distribution($\lambda$ = %d)' % lam)

plt.subplots_adjust(left=0.15)
plt.show()

numpyの関数を使って乱数生成

import numpy as np

def rpoisson(lam, n):
    return np.random.poisson(lam, n)

lam = {1, 4, 10}と変化させて乱数を生成した。サンプル数は1,000,000。
以下がヒストグラム。赤色の線は理論値。

ポワソン分布にしたがって乱数が生成されていることが分かる。

自作MCMCを使って乱数生成

import math
from random import randint, random
from itertools import islice

def rpoisson(lam, n):
    
    def p(lam, k):
        if k < 0:
            return 0.0
        ret = math.exp(-lam)
        for i in xrange(k):
            ret *= 1.0 * lam / (i + 1)
        return ret

    def metropolis(lam):
        x = lam
        while True:
            yield x
            xt = (x-1) if randint(0, 1) == 0 else (x+1)
            alpha = p(lam, xt) / p(lam, x)
            if random() < alpha:
                x = xt

    return list(islice(metropolis(lam), n))

次の乱数候補には、50%ずつの確率で1大きい数 or 1小さい数を採用している。
おそらくこのような単純なやり方ではダメだろうと思っていたが、意外とそれらしい感じで乱数を生成できた。