Help us understand the problem. What is going on with this article?

🍺Python3でポアソン分布、ガンマ分布

More than 1 year has passed since last update.

卒論でポアソン分布とガンマ分布を扱ったのでメモしておく。

ポアソン分布

ポアソン分布では単位時間あたり平均$λ$回発生するようなランダムな事象があるときに、時刻$t$までに$n$回起こる確率を表現することができる。
その確率を以下のように表せる。

  P(n,t)=\frac{\mathrm{e}^{-{\lambda}t}({\lambda}t)^n}{n!}

実装は以下のようになる

python3を用いて実装を行なった。

import math
import matplotlib.pyplot as plt

# ポアソン分布を計算する
# 返り値は単位時間あたり平均λ回発生するようなランダムな事象があるときに、時刻tまでにn回起こる確率
def poisson_probability(n, t, lambda_poisson):
    return math.e**(-lambda_poisson * t) * (lambda_poisson * t)**n / math.factorial(n)

# 累積ポアソン分布を計算する
# nは離散であるので積分は用いずΣを用いて計算する
# 返り値は単位時間あたり平均λ回発生するようなランダムな事象があるときに、時刻tまでに0,1,2,3,4...n回起こる確率
def cumulative_poisson_probability(n, t, lambda_poisson):
    cumulative_probability = 0
    for i in range(0, n + 1):
        cumulative_probability += math.e**(-1 * lambda_poisson * t) * (lambda_poisson * t)**i / math.factorial(i)
    return cumulative_probability

# ポアソン分布をプロットする
# x軸:n, y軸:確率
# 単位時間あたり平均λ回発生するようなランダムな事象があるときに、時刻tまでにn回起こる確率をnを変動させて出力する
def plot_poisson(time, lambda_poisson):
    x_axis = np.linspace(0, 2 * time * lambda_poisson, 2 * time * lambda_poisson + 1)
    y_axis = [poisson_probability(x, time, lambda_poisson) for x in x_axis]
    plt.title('poisson time: {0} lambda: {1}'.format(time, lambda_poisson))
    plt.xlabel('n')
    plt.ylabel('probability')
    plt.plot(x_axis, y_axis)
    plt.savefig("poisson.png")
    plt.show()

# 累積ポアソン分布をプロットする
# x軸:n, y軸:確率
# 単位時間あたり平均λ回発生するようなランダムな事象があるときに、時刻tまでに0,1,2,3...n回起こる確率をnを変動させて出力する
def plot_cumulative_poisson(time, lambda_poisson):
    x_axis = np.linspace(0, 2 * time * lambda_poisson, 2 * time * lambda_poisson + 1)
    y_axis = []
    for x in x_axis:
        y_axis.append(cumulative_poisson_probability(int(x), time, lambda_poisson))
    plt.title('cumulative poisson time: {0} lambda: {1}'.format(time, lambda_poisson))
    plt.xlabel('n')
    plt.ylabel('probability')
    plt.plot(x_axis, y_axis)
    plt.show()

poisson.png
上図はポアソン分布である。単位時間あたり0.5回発生するランダムな事象は時刻50までに発生する回数の確率を表している。直感通り、$n=25$あたりでピークとなっている。
c_poisson.png
上図は累積ポアソン分布である。累積して確率の和は1に収束する。

ガンマ分布

ガンマ分布では単位時間あたり平均$λ$回発生するようなランダムな事象があるときに、時刻$tに$n$回目が起こる確率を表現することができる。
その確率は以下のように表せる。

f(n,t) = \frac{{\lambda}^n}{(n-1)!}t^{n-1}\mathrm{e}^{-{\lambda}t}

実装は以下のようになる

python3を用いて実装を行なった。

import math
import matplotlib.pyplot as plt

# ガンマ分布を計算する
# 返り値は単位時間あたり平均λ回発生するようなランダムな事象があるときに、時刻tにn回目が起こる確率
def gamma_probability(n, t, lambda_poisson):
    return lambda_poisson**n * t**(n - 1) * math.e**(-lambda_poisson * t) / math.factorial(n - 1)

# 累積ガンマ分布を計算する
# 時間tは連続であるので積分を用いて計算する
# 返り値は単位時間あたり平均λ回発生するようなランダムな事象があるときに、時刻0~tにn回目が起こる確率
def cumulative_gamma_probability(n, T, lambda_poisson):
    cumulative_probability, abserr = integrate.quad(lambda t: gamma_probability(n, t, lambda_poisson), 0, T)
    return cumulative_probability!

# ガンマ分布をプロットする
# x軸:t, y軸:確率
# 単位時間あたり平均λ回発生するようなランダムな事象があるときに、時刻tにn回目が起こる確率をプロットする
def plot_gamma(n, lambda_poisson):
    x_axis = np.linspace(0, 2 * n / lambda_poisson , 10 * n / lambda_poisson)
    y_axis = [gamma_probability(n, x, lambda_poisson) for x in x_axis]
    plt.title('gamma n: {0} lambda: {1}'.format(n, lambda_poisson))
    plt.xlabel('time')
    plt.ylabel('probability')
    plt.plot(x_axis, y_axis)
    plt.savefig("gamma.png")
    plt.show()

# 累積ガンマ分布をプロットする
# x軸:t, y軸:確率
# 単位時間あたり平均λ回発生するようなランダムな事象があるときに、時刻0~tにn回目が起こる確率をプロットする
def plot_cumulative_gamma(n, lambda_poisson):
    x_axis = np.linspace(0, 2 * n / lambda_poisson , 2 * n / lambda_poisson + 1)
    y_axis = []
    for x in x_axis:
        y_axis.append(cumulative_gamma_probability(n, x, lambda_poisson))
    plt.title('cumulative gamma n: {0} lambda: {1}'.format(n, lambda_poisson))
    plt.xlabel('time')
    plt.ylabel('probability')
    plt.plot(x_axis, y_axis)
    plt.show()

gamma.png
上図はガンマ分布である。単位時間あたり0.5回発生するランダムな事象が30回目に発生する時の確率を表している。直感通り、$t=60$あたりでピークとなっている。
c_gamma.png
上図は累積ガンマ分布である。$t$は連続であるので積分を行うと確率の和は1に収束する。

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away