「入門はじめての統計解析」をPythonで写経 Vol.3 ～ ２章「はじめての確率分布」①二項分布、超幾何分布、ポアソン分布

２章「はじめての確率分布」

書籍の著者　石村貞夫 先生

この記事は、書籍「入門はじめての統計解析」２章「はじめての確率分布」の Python写経活動 を取り扱います。

書籍の図・表・計算を淡々とPython化する写経シリーズです。
２章の前編では離散型確率分布を中心に取り組みます。
二項分布、超幾何分布、ポアソン分布などを実践します。
では書籍を開いて統計解析の旅に出発です🚀

はじめに

このブログシリーズは、書籍「入門はじめての統計解析」（東京図書、「テキスト」と呼びます）の Python 写経を通じて得た「統計解析の楽しさ」をご紹介します。

書籍の紹介と引用表記はリンク先の記事に掲載しています。

２章 はじめての確率分布

この記事は２章の以下のSectionを取り扱います。

記事に用いるデータは、テキストに掲載されたデータそのものを引用しています。
データ件数の少ないものはコード上でデータを登録し、データ件数の多いものはCSVファイル化してデータを読み込みしています。

２章で用いるライブラリをインポートします。

### インポート

# 数値計算
import math
import numpy as np
import pandas as pd
from scipy.special import gamma  # ガンマ関数

# 統計
import scipy.stats as stats

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'

Section 2.1 確率変数と確率分布

p.66 表2.1.1「身長の度数分布表」データを作成します。
相対度数に注目します。

### データの作成 p.66 表2.1.1

# 階級値と度数でデータフレームを作成
data1 = pd.DataFrame(
    {'階級値': np.arange(142.5, 167.6, 5),
     '度数': [1, 2, 13, 24, 14, 6]}
)
# 相対度数列の追加
data1['相対度数'] = data1['度数'] / data1['度数'].sum()

# 結果の表示
print('data1.shape: ', data1.shape)
data1

【実行結果】

テキストによると、度数分布表の相対度数は合計１になるので「確率の条件を満たしているようですね！」とされています。

合計をチェックしましょう。

# 合計チェック
data1.sum(axis=0).rename('合計').to_frame().T

【実行結果】
相対度数の合計は１です。
度数分布表は確率とつながっているのですね！

■ 離散型確率分布の平均と分散 p.68
度数分布表の「階級値」は連続していなくて「とびとび」の値になっています。
とびとびの値に関する確率を取り扱うのが「離散型確率分布」です。
階級値のように確率値と個別に関連する変数を「確率変数」と呼びます。
ひとまず、確率分布は確率変数の値と確率値のマッピング、と捉えておきましょう。

p.68 表2.1.4「確率分布の平均$${\mu}$$と分散$${\sigma^2}$$」のデータを作成します。

### 例 表2.1.1の平均μと分散σ² p.68 表2.1.4

# data1の確率変数と相対度数をコピー、列名をXとPに変更
data1_stats = data1.drop(columns=['度数']).copy()
data1_stats.columns = ['X', 'P']

# 列xp, x-μ, (x-μ)², (x-μ)²*pを追加
data1_stats['xp'] = data1_stats['X'] * data1_stats['P']
data1_stats['x-μ'] = data1_stats['X'] - data1_stats['xp'].sum()
data1_stats['(x-μ)²'] = data1_stats['x-μ']**2
data1_stats['(x-μ)²×p'] = data1_stats['(x-μ)²'] * data1_stats['P']

# 結果の表示
data1_stats

【実行結果】
確率$${P}$$の値の小数点の丸めを行っていないため、テキストの数値とズレています。

表の項目内容は以下のようになっています。

• X：確率変数の値（最初の表の階級値が相当）

• P：確率

• xp：確率変数の値 × P ⇒ 合計が「期待値（平均）」

• x-μ：確率変数の値－平均

• (x-μ)²：x-μ の二乗

• (x-μ)²×p：(x-μ)²×確率 ⇒ 合計が「分散」 

度数分布表の確率値を元にした離散確率分布の「期待値（平均）」「分散」を確認します。
各列の合計を算出します。

# 合計値を表示
data1_stats.sum(axis=0).rename('合計').to_frame().T

【実行結果】
確率変数の値とその確率を掛け算した xp が $${158.0}$$ が期待値（平均）になります。
(x-μ)²×p の $${28.08}$$ が分散になります。

平均（期待値）$${\mu}$$ と分散 $${\sigma^2}$$ を別の方法で算出します。

### 期待値μの算出 p.68
mu = sum(data1_stats['X'] * data1_stats['P'])
print(f'期待値 {mu}')

【実行結果】
表の合計による期待値の値と同じになりました。

### 分散σ²の算出 p.68
sigma2 = sum((data1_stats['X'] - mu)**2 * data1_stats['P'])
print(f'分散 {sigma2:.2f}')

【実行結果】
表の合計による分散の値と同じになりました。
なお、テキストの分散と若干異なります。

■ 連続型確率変数と連続型確率分布 p.69
連続型確率変数は「とびとび」ではなく「連続値」である確率変数です。
連続型確率変数に対応する確率分布が連続型確率分布です。

p.69 図2.1.2「確率変数の区間とその確率」の図を描画します。
連続型確率分布には正規分布 $${N(0, 3)}$$ を用います。
確率計算には scipy.stats を利用します。

### 確率変数の区間とその確率 p.69 図2.1.2

## 設定と準備
a, b = -2.5, 3.5                        # 区間a,bの値
norm_dist = stats.norm(loc=0, scale=3)  # 平均0, 標準偏差3の正規分布

x_val = np.linspace(-10, 10, 1001)      # x軸の値
y_val = norm_dist.pdf(x_val)            # xに対応する正規分布の確率密度関数
ab_x_val = np.linspace(a, b, 1001)      # [a, b]のx軸の値
ab_y_val = norm_dist.pdf(ab_x_val)      # [a, b]のxに対応する正規分布の確率密度関数

## 描画
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 3))
# y=0の水平線（黒）の描画
ax.axhline(0, color='black', lw=0.5)
# 区間[a, b]確率Pの塗りつぶし描画
ax.fill_between(ab_x_val, 0, ab_y_val, alpha=0.3)
# 正規分布の確率密度関数（青）の描画
ax.plot(x_val, y_val)
# X=a, bの点の描画
ax.plot([a, b], [0, 0], 'o', color='tab:blue')
# 文字a, b、Xの表示
ax.text(a, -0.012, s='a', ha='center')
ax.text(b, -0.012, s='b', ha='center') 
ax.text(0, -0.012, s='$X$', fontsize=12, ha='center')
# P(a≦X≦b)の表示
ax.text(-2.1, 0.04, s=r'$P(a \leq X \leq b)$', fontsize=14)
# 枠線の消去
ax.axis('off');

【実行結果】
横軸が確率変数 $${X}$$、縦軸が確率密度です（確率密度って？）。
青い線がつながっています。連続しているのです！
また、青く塗りつぶした領域が区間 $${a, b}$$ の面積であり、確率 $${P(a \leq X \leq b)}$$です（こちらは確率です）。

■ 分布関数と確率密度関数 p.70
p.70 図2.1.4「確率密度関数の確率 $${P(a \leq X \leq b)}$$」を描画します。
連続型確率分布には正規分布 $${N(0, 3)}$$ を適用します。
確率計算には scipy.stats を利用します。

### 確率密度関数の確率 P(a≦X≦b) p.70 図2.1.4

## 設定と準備
a, b = -2.5, 3.5                        # 区間a,bの値
norm_dist = stats.norm(loc=0, scale=3)  # 平均0, 標準偏差3の正規分布

x_val = np.linspace(-10, 10, 1001)      # x軸の値
y_val = norm_dist.pdf(x_val)            # xに対応する正規分布の確率密度関数
ab_x_val = np.linspace(a, b, 1001)      # [a, b]のx軸の値
ab_y_val = norm_dist.pdf(ab_x_val)      # [a, b]のxに対応する正規分布の確率密度関数

## 描画
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 3))
# y=0の水平線（黒）の描画
ax.axhline(0, color='black', lw=0.5)
# 区間[a, b]確率Pの塗りつぶし描画
ax.fill_between(ab_x_val, 0, ab_y_val, alpha=0.3)
# 正規分布の確率密度関数（青）の描画
ax.plot(x_val, y_val)
# X=a, bの点の描画
ax.plot([a, b], [0, 0], 'o', color='tab:blue')
# 文字a, b、Xの表示
ax.text(a, -0.012, s='a', ha='center')
ax.text(b, -0.012, s='b', ha='center') 
ax.text(0, -0.012, s='$X$', fontsize=12, ha='center')
# P(a≦X≦b)の表示
ax.text(-1.6, 0.04, s=r'$\int_a^b f(x) dx$', fontsize=14)
# 枠線の消去
ax.axis('off');

【実行結果】
青い曲線が「確率密度関数 $${f(x)}$$」です。
曲線上の特定の点が「確率密度」です。
そして、青く塗りつぶした領域である確率は、区間 $${a, b}$$ の面積は定積分で算出できます。
この面積が連続型確率分布の確率です。
連続型確率分布の確率は定積分で求められるのです！
連続型確率分布の確率密度は点、確率は区間で示される面です。

続いて「分布関数」を探ります。
p.70 図2.1.3「分布関数 $${F(x)=P(X \leq x)}$$」の図を描画します。
連続型確率分布には正規分布 $${N(0, 3)}$$ を用います。
確率計算には scipy.stats を利用します。

### 分布関数 F(x)=P(X≦x)` p.70 図2.1.3

## 設定と準備
x = 3.5                                 # xの値
norm_dist = stats.norm(loc=0, scale=3)  # 平均0, 標準偏差3の正規分布

x_val = np.linspace(-10, 10, 1001)      # x軸の値
y_val = norm_dist.pdf(x_val)            # x軸の値に対応する正規分布の確率密度関数
x_x_val = np.linspace(-10, x, 1001)    # [-10, x]のx軸の値
x_y_val = norm_dist.pdf(x_x_val)      # [-10, x]のxに対応する正規分布の確率密度関数

## 描画
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 3))
# y=0の水平線（黒）の描画
ax.axhline(0, color='black', lw=0.5)
# 区間[-10, x]確率Pの塗りつぶし描画
ax.fill_between(x_x_val, 0, x_y_val, alpha=0.3)
# 正規分布の確率密度関数（青）の描画
ax.plot(x_val, y_val)
# X=xの点の描画
ax.plot([x], [0], 'o', color='tab:blue')
# 文字a, b、Xの表示
ax.text(x, -0.012, s='x', ha='center')
ax.text(0, -0.012, s='$X$', fontsize=12, ha='center')
# P(X≦x)の表示
ax.text(-2.1, 0.04, s=r'$P(X \leq x)$', fontsize=14)
# 枠線の消去
ax.axis('off');

【実行結果】
青く塗りつぶした領域が $${X \leq x}$$ の面積であり、「分布関数 $${F(x)}$$ の値」です。
分布関数は、確率変数 $${X}$$ の具体的な値 $${x}$$ 以下の確率を示しています。

分布関数 $${F(x)}$$ と確率密度関数 $${f(x)}$$ の関係に関するテキストの数式をお借りします。

Section 2.2 二項分布 － 比率のときは

この Section から離散型確率分布の具体例を見ていきます。
確率変数が「とびとび」の場合の確率分布です。
二項分布から始めましょう。

二項分布の定義・平均・分散に関するテキストの数式をお借りします。
■ 二項分布 $${B(n,p)}$$ の定義 p.72
確率変数 $${X}$$ が $${0, 1, 2, \cdots, n}$$ の値を取るとき、

■ 二項分布の平均・分散の公式 p.72

二項分布の定義の数式を用いて、二項分布関数を【Python標準縛り】で作成します。

### 二項分布の関数定義

# 階乗関数
def factorial(x):
    return math.prod([i for i in range(1, x + 1)])

# 二項分布の確率算出関数
def binom_prob(x, n, p):
    # 組み合わせnCxの計算
    c = factorial(n) / (factorial(x) * factorial(n - x))
    # 戻り値: 確率の算出
    return c * p**x * (1 - p)**(n - x)

# テスト
binom_prob(x=3, n=8, p=0.3)

【実行結果】
二項分布 $${B(8, 0.3)}$$ に従う確率変数 $${X=3}$$ の確率を算出しました。

scipy.stats で検算します。

# scipy.statsで答え合わせ
stats.binom.pmf(k=3, n=8, p=0.3)

【実行結果】
検算結果は一致しました。

■ 二項分布の確率とグラフ p.73
p.73 例１）の $${n=10,\ p=0.5}$$ の場合の、表2.2.1「二項分布」（確率分布の表）と図2.2.1「二項分布」（棒グラフ）を作成します。
確率の計算には自作関数を利用します。

### 例1) n=10, p=0.5の場合 p.73 表2.2.1

# 設定
n, p = 10, 0.5     # 試行回数n, 成功確率p
x_val = range(11)  # Xのとりうる値
# 二項分布の確率計算を行い、データフレーム化
data2 = pd.DataFrame(
    {'X': x_val, 'P(X=x)': [binom_prob(x, n, p) for x in x_val]})
# 結果確認
data2

【実行結果】
確率変数 $${X}$$ の値はとびとびです。
$${P(X=x)}$$ は確率変数の具体的な値に対応する確率です。

### 例1) n=10, p=0.5の場合 p.73 図2.2.1
sns.barplot(data=data2, x='X', y='P(X=x)', width=0.98, alpha=0.7);

【実行結果】
棒グラフに描画したので、確率変数がつながって連続のように見えますが、離散型ですのでご注意下さい。
$${n=10,\ p=0.5}$$ の二項分布の確率のチャートは左右対称です！

scipy.stats で確率を算出して可視化してみましょう。

### scipy.statsでBin(10, 0.5)を可視化

# 設定
x_val = np.arange(0, 11)
# 確率の点の描画
plt.plot(x_val, stats.binom.pmf(k=x_val, n=n, p=p), 'o', ms=8)
# 確率の垂直線の描画
plt.vlines(x_val, 0, stats.binom.pmf(k=x_val, n=n, p=p), lw=2)
# 修飾
plt.xlabel('X', fontsize=12)
plt.ylabel('P(X=x)', fontsize=12);

【実行結果】
点と細い棒で表すと離散型確率分布を感じられると思います！

続いて、p.73 例２）の $${n=10,\ p=0.1}$$ の場合の、表2.2.2「二項分布」（確率分布の表）と図2.2.2「二項分布」（棒グラフ）を作成します。
確率の計算には自作関数を利用します。

### 例2) n=10, p=0.1の場合 p.73 表2.2.2

# 設定
n, p = 10, 0.1     # 試行回数n, 成功確率p
x_val = range(11)  # Xのとりうる値
# 二項分布の確率計算を行い、データフレーム化
data3 = pd.DataFrame(
    {'X': x_val, 'P(X=x)': [binom_prob(x, n, p) for x in x_val]})
# 結果確認
data3.round(5)

【実行結果】

### 例1) n=10, p=0.5の場合 p.73 図2.2.2
sns.barplot(data=data3, x='X', y='P(X=x)', width=0.97, alpha=0.7);

【実行結果】
$${n=10,\ p=0.1}$$ の二項分布の確率のチャートは、左右非対称です！

scipy.stats で確率を算出して可視化してみましょう。

### scipy.statsでBin(10, 0.1)を可視化

# 設定
x_val = np.arange(0, 11)
# 確率の点の描画
plt.plot(x_val, stats.binom.pmf(k=x_val, n=n, p=p), 'o', ms=8)
# 確率の垂直線の描画
plt.vlines(x_val, 0, stats.binom.pmf(k=x_val, n=n, p=p), lw=2)
# 修飾
plt.xlabel('X', fontsize=12)
plt.ylabel('P(X=x)', fontsize=12);

【実行結果】

■ 二項分布の例 p.74
p.74 例）不良品の割合 $${p=0.03}$$ の電子部品について、箱から $${n}$$ 個の標本を復元抽出して $${x}$$ 個の不良品が含まれる確率を可視化します。
確率計算には scipy.stats を利用します。

### 例：不良品の割合p=0.03の電子部品について、
# 箱からn個の標本を復元抽出してx個の不良品が含まれる確率 p.74 ※scipy.stats利用

# 設定
n, p = 100, 0.03
x_val = np.arange(0, 11)
# 確率の点の描画
plt.plot(x_val, stats.binom.pmf(k=x_val, n=n, p=p), 'o', ms=8)
# 確率の垂直線の描画
plt.vlines(x_val, 0, stats.binom.pmf(k=x_val, n=n, p=p), lw=2)
# 修飾
plt.xlabel('X: 不良品数', fontsize=12)
plt.ylabel('P(X=x)', fontsize=12)
plt.title(f'標本サイズ{n}, 不良率{p}');

【実行結果】

p.74 例）コイントスを$${n}$$回行って表の出る回数を $${x}$$ とするときの確率を可視化します。
確率計算には scipy.stats を利用します。

### 例：コイントスをn回行って表の出る回数をxとする p.74 ※scipy.stats利用

# 設定
n, p = 10, 0.5
x_val = np.arange(0, 11)
# 確率の点の描画
plt.plot(x_val, stats.binom.pmf(k=x_val, n=n, p=p), 'o', ms=8)
# 確率の垂直線の描画
plt.vlines(x_val, 0, stats.binom.pmf(k=x_val, n=n, p=p), lw=2)
# 修飾
plt.xlabel('X: 表の出る回数', fontsize=12)
plt.ylabel('P(X=x)', fontsize=12)
plt.title(f'コイントス{n}回, 表の出る確率{p}');

【実行結果】

Section 2.3 超幾何分布 － 非復元抽出のときは

超幾何分布の定義・平均・分散に関するテキストの数式をお借りします。
■ 超幾何分布の定義 p.75
確率変数 $${X}$$  が$${0, 1, 2, \cdots, n}$$ の値を取るとき、

■ 超幾何分布の平均・分散の公式 p.75

超幾何分布の定義の数式を用いて、超幾何分布関数を【Python標準縛り】で作成します。

### 超幾何分布の関数定義 p.75

# 階乗関数
def factorial(x):
    return math.prod([i for i in range(1, x + 1)])

# 超幾何分布の確率算出関数
def hyper_geom_prob(x, N, n, p):
    # N * pを計算して整数化
    Np = int(N * p)
    # 組み合わせNp_C_xの計算
    c1 = factorial(Np) / (factorial(x) * factorial(Np - x))
    # 組み合わせN-Np_C_n-xの計算
    c2 = factorial(N - Np) / (factorial(n - x) * factorial(N - Np - (n - x)))
    # 組み合わせN_C_nの計算
    c3 = factorial(N) / (factorial(n) * factorial(N - n))
    # 戻り値: 確率の算出
    return (c1 * c2) / c3

# テスト
hyper_geom_prob(x=3, N=100, n=8, p=0.3)

【実行結果】
$${N=100,\ p=0.3,\ n=8}$$の超幾何分布に従う確率変数$${X=3}$$の確率を算出しました。

scipy.stats で検算します。

### scipy.statsで答え合わせ
x, N, n, p = 3, 100, 8, 0.3
stats.hypergeom.pmf(k=x, M=N, n=N*p, N=n)

【実行結果】
検算結果は一致しました。

■ 超幾何分布の確率とグラフ p.76
p.76 例１）の $${N=100,\ p=0.1,\ n=10}$$ の場合の、表2.3.1「超幾何分布」（確率分布の表）と図2.3.1「超幾何分布」（棒グラフ）を作成します。
確率の計算には自作関数を利用します。

### 例1) N=100, p=0.1, n=10の場合 p.76 表2.3.1

# 設定
N, p, n = 100, 0.1, 10     # 全数N, 割合p, 標本サイズn
x_val = range(8)           # Xのとりうる値
# 超幾何分布の確率計算を行い、データフレーム化
data4 = pd.DataFrame(
    {'X': x_val, 'P(X=x)': [hyper_geom_prob(x, N, n, p) for x in x_val]})
# 結果確認
data4.round(5)

【実行結果】

### 例1) N=100, p=0.1, n=10の場合 p.76 図2.3.1
sns.barplot(data=data4, x='X', y='P(X=x)', width=0.98, alpha=0.7);

【実行結果】
$${N=100,\ p=0.1,\ n=10}$$ の超幾何分布の確率のチャートは、左右非対称です。

scipy.stats で確率を算出して可視化してみましょう。

### scipy.statsでHyperGeo(N=100, p=0.1, n=10)を可視化

# 設定
x_val = np.arange(0, 11)
# 確率の点の描画
plt.plot(x_val, stats.hypergeom.pmf(k=x_val, M=N, n=N*p, N=n), 'o', ms=8)
# 確率の垂直線の描画
plt.vlines(x_val, 0, stats.hypergeom.pmf(k=x_val, M=N, n=N*p, N=n), lw=2)
# 修飾
plt.xlabel('X', fontsize=12)
plt.ylabel('P(X=x)', fontsize=12);

【実行結果】

続いて、p.76 例２）の $${N=1000,\ n=10,\ p=0.1}$$ の場合の、表2.3.2「超幾何分布」（確率分布の表）と図2.3.2「超幾何分布」（棒グラフ）を作成します。
確率の計算には自作関数を利用します。

### 例2) N=1000, p=0.5, n=10の場合 p.76 表2.3.2

# 設定
N, p, n = 1000, 0.5, 10    # 全数N, 割合p, 標本サイズn
x_val = range(11)          # Xのとりうる値
# 超幾何分布の確率計算を行い、データフレーム化
data5 = pd.DataFrame(
    {'X': x_val, 'P(X=x)': [hyper_geom_prob(x, N, n, p) for x in x_val]})
# 結果確認
data5.round(5)

【実行結果】

### 例2) N=1000, p=0.5, n=10の場合 p.76 図2.3.2
# ※Nの値が大きくなると超幾何分布は二項分布に近づく

sns.barplot(data=data5, x='X', y='P(X=x)', width=0.98, alpha=0.7);

【実行結果】
$${N=1000,\ n=10,\ p=0.1}$$ の超幾何分布の確率のチャートは、ほぼ左右対称です！

scipy.stats で確率を算出して可視化してみましょう。

### scipy.statsでHyperGeo(N=1000, p=0.5, n=10)を可視化

# 設定
x_val = np.arange(0, 11)
# 確率の点の描画
plt.plot(x_val, stats.hypergeom.pmf(k=x_val, M=N, n=N*p, N=n), 'o', ms=8)
# 確率の垂直線の描画
plt.vlines(x_val, 0, stats.hypergeom.pmf(k=x_val, M=N, n=N*p, N=n), lw=2)
# 修飾
plt.xlabel('X', fontsize=12)
plt.ylabel('P(X=x)', fontsize=12);

【実行結果】

■ 超幾何分布の例 p.77
p.77 例）不良品の割合 $${p=0.03}$$ の電子部品について、$${N}$$ 個入っている箱から $${n}$$ 個の標本を非復元抽出して $${x}$$ 個の不良品が含まれる確率を可視化します。
確率計算には scipy.stats を利用します。

### 例：不良品の割合p=0.03の電子部品について、
# N個入っている箱からn個の標本を非復元抽出してx個の不良品が含まれる確率 p.77
# ※scipy.stats利用

# 設定
N, p, n = 100, 0.03, 10
x_val = np.arange(0, 11)
# 確率の点の描画
plt.plot(x_val, stats.hypergeom.pmf(k=x_val, M=N, n=N*p, N=n), 'o', ms=8)
# 確率の垂直線の描画
plt.vlines(x_val, 0, stats.hypergeom.pmf(k=x_val, M=N, n=N*p, N=n), lw=2)
# 修飾
plt.xlabel('X: 不良品数', fontsize=12)
plt.ylabel('P(X=x)', fontsize=12)
plt.title(f'全数{N}, 標本サイズ{n}, 不良率{p}');

【実行結果】

p.77 例）児童数 $${N=715}$$ 名の小学校で、虫歯のある児童の割合 $${p=0.6}$$、$${n=10名}$$ の児童を選んだときに虫歯のある児童が $${x}$$ 人含まれる確率を可視化します。
確率計算には scipy.stats を利用します。

### 例：児童数N=715名の小学校で、虫歯のある児童の割合p=0.6、
# 10名の児童を選んだときに虫歯のある児童がx人含まれる確率 p.77 ※scipy.stats利用

# 設定
N, p, n = 715, 0.6, 10
x_val = np.arange(0, 13)
# 確率の点の描画
plt.plot(x_val, stats.hypergeom.pmf(k=x_val, M=N, n=N*p, N=n), 'o', ms=8)
# 確率の垂直線の描画
plt.vlines(x_val, 0, stats.hypergeom.pmf(k=x_val, M=N, n=N*p, N=n), lw=2)
# 修飾
plt.xlabel('X: 虫歯の児童数', fontsize=12)
plt.ylabel('P(X=x)', fontsize=12)
plt.title(f'全数{N}, 標本サイズ{n}, 虫歯率{p}');

【実行結果】

Section 2.4 ポアソン分布 － めったに起こらないときは

ポアソン分布の定義・平均・分散に関するテキストの数式をお借りします。
■ ポアソン分布 $${P(\lambda)}$$ の定義 p.78
確率変数 $${X}$$ が $${0, 1, 2, \cdots, n}$$ の値を取るとき、

■ ポアソン分布の平均・分散の公式 p.78

ポアソン分布の定義の数式を用いて、ポアソン分布関数を【Python標準縛り】で作成します。

### ポアソン分布の関数定義 p.78

# 階乗関数
def factorial(x):
    return math.prod([i for i in range(1, x + 1)])

# ポアソン分布の確率算出関数
def poisson_prob(x, lam):
    # 戻り値: 確率の算出
    return lam**x / factorial(x) * np.exp(-lam)

# テスト
poisson_prob(x=3, lam=5)

【実行結果】
ポアソン分布 $${P(5)}$$ に従う確率変数 $${X=3}$$ の確率を算出しました。

scipy.stats で検算します。

### scipy.statsで答え合わせ
stats.poisson.pmf(k=3, mu=5)

【実行結果】
検算結果は一致しました。

■ ポアソン分布の確率とグラフ p.78
p.76 例）の $${\lambda=1}$$ の場合の、表2.4.1「ポアソン分布」（確率分布の表）と図2.4.1「ポアソン分布」（棒グラフ）を作成します。
確率の計算には自作関数を利用します。

### 例) λ=1の場合 p.78 表2.4.1

# 設定
lam = 1           # 平均λ
x_val = range(8)  # Xのとりうる値
# ポアソン分布の確率計算を行い、データフレーム化
data6 = pd.DataFrame(
    {'X': x_val, 'P(X=x)': [poisson_prob(x, lam) for x in x_val]})
# 結果確認
data6.round(5)

【実行結果】

### 例) λ=1の場合 p.78 図2.4.1
sns.barplot(data=data6, x='X', y='P(X=x)', width=0.98, alpha=0.7);

【実行結果】
$${\lambda=1}$$ のポアソン分布の確率のチャートは、左右非対称です。

テキストによるとポアソン分布は・・・

• ポアソン分布と二項分布の関係

  • 二項分布において、$${np}$$ を一定値 $${\lambda}$$ に固定したまま、$${n \rightarrow \infty}$$、つまり $${p \rightarrow 0}$$ としたときの極限分布がポアソン分布になる

• ポアソン分布の特徴は、めったに起こらないような場合によく当てはまる分布である

■ ポアソン分布の例 p.79
p.79 例）19世紀のプロシヤで、ある１年間に馬に蹴られて死亡した兵士の数とその軍団の数から、1年間に馬に蹴られて死亡した一軍団の兵士数の平均を計算します。
表2.4.2「馬に蹴られて死亡した兵士の相対度数」を作成します。

### 馬に蹴られて死亡した兵士の相対度数 p.79 表2.4.2

# データフレームの作成の作成
data7 = pd.DataFrame([109, 65, 22, 3, 1, 0],
                     columns=['1年間に兵士がx人死亡した軍団の数'],
                     index=[0, 1, 2, 3, 4, '5以上']
)
# 相対度数列の追加
data7['相対度数x/200'] = data7['1年間に兵士がx人死亡した軍団の数'] / 200
# インデックスの列名の追加
data7.index.name = '死亡した兵士数x'
# 結果表示
data7.T

【実行結果】

1年間に馬に蹴られて死亡した一軍団の兵士数の平均を算出します。

### 1年間に馬に蹴られて死亡した一軍団の兵士数の平均
(data7.index[0:5] * data7.iloc[0:5, 1]).sum()

【実行結果】
1年間に馬に蹴られて死亡した一軍団の兵士数の平均は 0.61 です。

$${\lambda=0.61}$$ のポアソン分布の確率分布の表を作成します。
表2.4.3「$${\lambda=0.61}$$ のときのポアソン分布」に相当します。

### λ=0.61のときのポアソン分布 p.80 表2.4.3

# ポアソン分布の確率計算を行い、データフレーム化
data8 = pd.DataFrame(
    {'P(X=x)': [poisson_prob(x=x, lam=0.61) for x in range(8)]})
# インデックスに列名を追加
data8.index.name = 'X=x'
# 結果確認
data8.round(5).T

【実行結果】
表2.4.2「馬に蹴られて死亡した兵士の相対度数」の相対度数と近似しています。

今回の写経は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。