二項分布の極限がポアソン分布になることをグラフによりイメージする

二項分布の極限はポアソン分布になります。どういうことだか最初意味がわからなかったので、それをイメージできるようにグラフを作ってみました。

本記事は次の3つの記事の3つ目です。

今度はグラフでイメージします。

二項分布
- 母数 (分布のパラメータ)
  - $\displaystyle n$ (0以上の整数)
  - $\displaystyle p$ (0以上1以下の実数)
- 確率質量関数: $\displaystyle P(k) = {}_nC_k p^k (1-p)^{n-k}$
- 期待値: $\displaystyle np$
- 分散: $\displaystyle np(1-p)$
ポアソン分布
- 母数 (分布のパラメータ)
  - $\displaystyle \lambda$ (0より大きい実数)
- 確率質量関数: $\displaystyle P(k) = \frac{\lambda^k e^{-\lambda}}{k!}$
- 期待値: $\displaystyle \lambda$
- 分散: $\displaystyle \lambda$

二項分布がポアソン分布に近づく様子のグラフです。

f:id:suzuki-navi:20200528184108p:plain

このグラフは6個の分布の確率質量関数です。横軸は $\displaystyle k$ 、縦軸が確率質量関数 $\displaystyle P(k)$ です。

6個の分布は次のとおりです。ピーク付近で上から順にです。

青: $\displaystyle n=10,p=\frac{5}{10}$ の二項分布
橙: $\displaystyle n=13,p=\frac{5}{13}$ の二項分布
緑: $\displaystyle n=17,p=\frac{5}{17}$ の二項分布
赤: $\displaystyle n=25,p=\frac{5}{25}$ の二項分布
紫: $\displaystyle n=50,p=\frac{5}{50}$ の二項分布
茶: $\displaystyle \lambda=5$ のポアソン分布

数式で書くと、

青: $\displaystyle P(k) = \left. {}_nC_k p^k (1-p)^{n-k} \right|_{n=10,p=\frac{5}{10}}$
橙: $\displaystyle P(k) = \left. {}_nC_k p^k (1-p)^{n-k} \right|_{n=13,p=\frac{5}{13}}$
緑: $\displaystyle P(k) = \left. {}_nC_k p^k (1-p)^{n-k} \right|_{n=17,p=\frac{5}{17}}$
赤: $\displaystyle P(k) = \left. {}_nC_k p^k (1-p)^{n-k} \right|_{n=25,p=\frac{5}{25}}$
紫: $\displaystyle P(k) = \left. {}_nC_k p^k (1-p)^{n-k} \right|_{n=50,p=\frac{5}{50}}$
茶: $\displaystyle P(k) = \left. \frac{\lambda^k e^{-\lambda}}{k!} \right|_{\lambda=5}$

二項分布の $\displaystyle np$ を一定に保ったまま、 $\displaystyle n$ を大きくして、 $\displaystyle p$ を小さくしていくと、 $\displaystyle \lambda=np$ のポアソン分布に近づいていくことがグラフからもわかります。

このグラフを書いたコードは以下です。たまにはPHPを書かないと忘れてしまって悲しいのでPHPで久しぶりにコードを書きました。以前はPHPもよく書いてましたが、いまは書くことがほとんどないです。

<?php

function combination($n, $k) {
  $ret = 1;
  for ($i = $n - $k + 1; $i <= $n; $i++) {
    $ret *= $i;
  }
  for ($i = 2; $i <= $k; $i++) {
    $ret /= $i;
  }
  return $ret;
}

function factorial($n) {
  $ret = 1;
  for ($i = 2; $i <= $n; $i++) {
    $ret *= $i;
  }
  return $ret;
}

function binomial($lambda, $n, $k) {
  $p = $lambda / $n;
  return combination($n, $k) * pow($p, $k) * pow(1 - $p, $n - $k);
}

function poisson($lambda, $k) {
  return pow($lambda, $k) * exp(-$lambda) / factorial($k);
}

for ($i = 0; $i <= 10; $i++) {
  printf("%d,%f,%f,%f,%f,%f,%f\n", $i,
    binomial(5, 10, $i),
    binomial(5, 13, $i),
    binomial(5, 17, $i),
    binomial(5, 25, $i),
    binomial(5, 50, $i),
    poisson(5, $i));
}

標準出力にCSVで数値を出力しますので、これをファイルに保存します。

$ php binomial-poisson.php > binomial-poisson.csv

Jupyter Notebook上で、PythonのpandasでCSVファイルを読み込んで、matplotlibでグラフにしました。

%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

df = pd.read_csv("binomial-poisson.csv", names=["x", "10", "13", "17", "25", "50", "poisson"])
plt.plot(df["x"], df["10"])
plt.plot(df["x"], df["13"])
plt.plot(df["x"], df["17"])
plt.plot(df["x"], df["25"])
plt.plot(df["x"], df["50"])
plt.plot(df["x"], df["poisson"])

以上。