データサイエンス講座(統計編) 18〜21 確率論（かめ@米国データサイエンティスト）のノート

2026-04-29statistics Python NumPy SciPy matplotlib

「データサイエンス講座(統計編)」（かめ@米国データサイエンティストさん）のノートです。多変数の記述統計学がテーマだった前回の続きになります。第18回「確率変数」から第21回「正規分布」まで、内容は 確率論 となります。

今回のハイライトは、ポアソンの少数の法則 や 二項分布の正規近似 などが、実際に「体験」できた事です。コンピューターを利用すれば、実際に母数を動かしながら確率分布を描くことができます。また、それぞれの確率分布について、乱数を生成して期待値と分散を計算しています。理論が実験で確かめられる、これがコンピューターの魅力ですね！

かめ@米国データサイエンティストさんの講座のお供に、今回も是非よろしくお願いします！JupyterLabのノートブックはこちらにあります。

19 二項分布をPythonを使って理解する(最も基本的な確率分布)

二項分布をPythonを使って理解する(最も基本的な確率分布)【統計学入門19】 - 米国データサイエンティストのブログ

二項分布
- PMF: scipy.stats.binom.pmf(k, n, p, loc=0)
- 乱数生成: scipy.stats.binom.rvs(n, p, loc=0, size=1, random_state=None)
- 期待値: scipy.stats.binom.mean(n, p, loc=0)
- 分散: scipy.stats.binom.var(n, p, loc=0)

scipy.stats.binom — SciPy v1.17.0 Manual

確率質量関数 (PMF)

P(x) = {}_n C_x p^x (1-p)^{n-x}

def my_binom_pmf(x: int, n: int, p: float):
    return math.comb(n, x) * p**x * (1 - p) ** (n - x)

math --- 数学関数 — Python 3.14.4 ドキュメント

実行例 ( $n = 3, p = 1 / 6$ )

>>> n = 3
>>> p = 1 / 6
>>> [my_binom_pmf(x, n, p) for x in range(n + 1)]
[0.5787037037037038,
 0.34722222222222227,
 0.06944444444444445,
 0.0046296296296296285]

scipy.stats.binom.pmf() の実行例 ( $n = 3, p = 1 / 6$ )

>>> [stats.binom.pmf(x, n, p) for x in range(n + 1)]
[np.float64(0.5787037037037036),
 np.float64(0.3472222222222223),
 np.float64(0.06944444444444445),
 np.float64(0.0046296296296296285)]

scipy.stats.binom — SciPy v1.17.0 Manual

乱数生成：期待値と分散

母数 $n = 3, p = 1 / 6$ の二項分布に従う乱数（サンプルサイズ $2160$ ）を生成する。

>>> n = 3
>>> p = 1 / 6
>>> size = 2160
>>> rng = np.random.default_rng(1728)
>>> data_binom = stats.binom.rvs(n, p, size=size, random_state=rng)
>>> data_binom
array([1, 1, 1, ..., 2, 0, 1], shape=(2160,))

scipy.stats.binom — SciPy v1.17.0 Manual

ヒストグラム

>>> sns.set_theme()
>>> fig, axes = plt.subplots()
>>> sns.histplot(data=data_binom, kde=False, ax=axes)
>>> plt.show()

度数を理論値と比較

>>> np.unique(data_binom, return_counts=True)
(array([0, 1, 2, 3]), array([1233,  766,  150,   11]))
>>> [stats.binom.pmf(x, n, p) * size for x in range(n + 1)]
[np.float64(1249.9999999999998),
 np.float64(750.0000000000002),
 np.float64(150.0),
 np.float64(9.999999999999998)]

確かに二項分布に従って乱数が生成できるようだ。

二項分布の期待値は $np$ 、分散は $np(1-p)$ となる。これを生成した乱数で確かめてみる。

平均

>>> np.mean(data_binom)
np.float64(0.5087962962962963)
>>> stats.binom.mean(n, p) # 期待値
np.float64(0.5)

分散

>>> np.var(data_binom)
np.float64(0.4193670696159123)
>>> stats.binom.var(n, p) # 分散
np.float64(0.4166666666666667))

当たり前だが、ほぼ一致している。

20 ポワソン分布を完全解説！正規分布や二項分布との関係は？

ポワソン分布を完全解説！正規分布や二項分布との関係は？【統計学入門20】 - 米国データサイエンティストのブログ

注意: ポアソン分布の母数は $\lambda$ で表すことが多いが、ここでは SciPy の流儀に合わせて $\mu$ で表記する。

ポアソン分布
- PMF: scipy.stats.poisson.pmf(k, mu, loc=0)
- CDF: scipy.stats.poisson.cdf(k, mu, loc=0)
- 期待値: scipy.stats.poisson.mean(mu, loc=0)
- 分散: scipy.stats.poisson.var(mu, loc=0)
- 乱数生成: scipy.stats.poisson.rvs(mu, loc=0, size=1, random_state=None)

参照：scipy.stats.poisson — SciPy v1.17.0 Manual

少数の法則

\lim_{\substack{np = \mu \\ n \to \infty}} {}_n C_x p^x (1-p)^{n-x} = \frac{\mu^x}{x!} e^{-\mu}

$np = \mu$ を一定に保ちつつ $n$ を限りなく大きくしたとき、二項分布の確率は右辺に収束する。証明はポアソン分布【統計検定®準１級のための数学②】 | とけたろうブログなどを参考にして頂きたい。指数関数のマクローリン展開を考えると、この極限も確率分布となる事が分かる。これを ポアソン分布 と呼ぶ。

例：100個に1個の割合で不良品が発生する。500個ずつ箱に梱包したとき、1箱に不良品が5個入っている確率

このとき、不良品の個数 $x$ は $n = 500, p = 1 / 100$ の二項分布に従う。

>>> stats.binom.pmf(5, 500, 1 / 100)
np.float64(0.17635104507334973)

一方、ポアソン分布で計算すると

f(x) = \frac{\mu^x}{x!} e^{-\mu}

def my_poisson_pmf(x: int, mu: float):
    return mu**x * math.exp(-mu) / math.factorial(x)

実行例 ( $\mu = np = 5, x = 5$ )

>>> my_poisson_pmf(5, 5)
0.1754673697678507

scipy.stats.poisson.pmf() の実行例

>>> stats.poisson.pmf(k=5, mu=5)
np.float64(0.17546736976785068)

ポアソン分布で計算した確率は、確かに二項分布の確率に近い。

次に、実際に $n$ を大きくして、二項分布の確率分布をプロットしてみる。

>>> mu = 5
>>> n_list = [10, 20, 100]
>>>
>>> x = np.arange(11)
>>>
>>> fig, axes = plt.subplots()
>>> # 二項分布を描画
>>> for n in n_list:
>>>     axes.plot(x, stats.binom.pmf(k=x, n=n, p=mu / n), "o--", label=f"Binom (n={n})")
>>> # ポアソン分布を描画
>>> axes.vlines(x, 0, stats.poisson.pmf(k=x, mu=mu), label=f"Poisson")
>>> axes.legend()
>>> axes.set_xlabel("x")
>>> axes.set_ylabel("Probability")
>>> plt.show()

$n$ が大きくなるにつれ、二項分布がポアソン分布に近づいている事が確認できる。

確率分布

例: 10分間に平均5回鳴るコールセンターで、1時間に40回電話が鳴る確率は？

$\mu = 5 / 10 \times 60 = 30$ なので、

>>> mu = 30
>>> x = 40
>>> stats.poisson.pmf(k=x, mu=mu)
np.float64(0.013943463479967897)

$x=40$ の確率は約1.4%しかない。

確率分布を描画してみると、

>>> x = np.arange(51)
>>> prob_list = stats.poisson.pmf(x, mu)
>>>
>>> fig, axes = plt.subplots()
>>> axes.plot(x, prob_list, "o")
>>> axes.vlines(x, 0, prob_list)
>>> axes.set_xlabel("Nom. of calls wihtin a hour")
>>> axes.set_ylabel("Probability")
>>> plt.show()

最頻値となる $x=30$ でも、確率は7%程度である。「ちょうどその回数」となる確率は高くない。

>>> 1 - stats.poisson.cdf(k=(35 - 1), mu=mu)
np.float64(0.20269167451688286
>>> 1 - stats.poisson.cdf(k=(40 - 1), mu=mu)
np.float64(0.0462530376458421)

$x \ge 35$ となる確率は20%、 $x \ge 40$ となる確率は4.6%となる。離散確率なので、CDF（累積分布関数）の引数で1引いている事に注意！

期待値と分散

ポアソン分布の期待値と分散は共に $\mu$ となる。

>>> mu = 30
>>> stats.poisson.mean(mu)
np.float64(30.0)
>>> stats.poisson.var(mu)
np.float64(30.0)

乱数生成をして確かめてみる。

>>> size = 1000
>>> rng = np.random.default_rng(1728)
>>> data_poisson = stats.poisson.rvs(mu, size=size, random_state=rng)
>>>
>>> fig, axes = plt.subplots()
>>> sns.histplot(data=data_poisson, kde=True, ax=axes)
>>> plt.show()

>>> np.mean(data_poisson)
np.float64(30.191)
>>> np.var(data_poisson)
np.float64(29.694519)

共に $\mu = 30$ に近いことが確かめられた。

21 幾何分布と指数分布を解説！ポワソン分布との関係は？

幾何分布と指数分布を解説！ポワソン分布との関係は？【統計学入門21】 - 米国データサイエンティストのブログ

幾何分布
- PMF: scipy.stats.geom.pmf(k, p, loc=0)
- 期待値: scipy.stats.geom.mean(p, loc=0)
- 分散: scipy.stats.geom.var(p, loc=0)
- 乱数生成: scipy.stats.geom.rvs(p, loc=0, size=1, random_state=None)

参照：scipy.stats.geom — SciPy v1.18.0 Manual

指数分布（尺度母数 $1/\lambda$ は引数scaleで指定）
- PDF: scipy.stats.expon.pmf(x, loc=0, scale=1)
- CDF: scipy.stats.expon.cdf(x, loc=0, scale=1)
- 期待値: scipy.stats.expon.mean(loc=0, scale=1)
- 分散: scipy.stats.expon.var(loc=0, scale=1)
- 乱数生成: scipy.stats.expon.rvs(loc=0, scale=1, size=1, random_state=None)

参照：scipy.stats.expon — SciPy v1.18.0 Manual

幾何分布 (Geometric distribution)

P(x) = (1-p)^{x-1} p

補足: 試行回数 $x$ は最後の成功を含む。そのため $x$ の範囲は $x \ge 1$ の整数となる。一方で最後の成功を含めない（失敗の回数のみ）場合もある。その時は $x \ge 0$ の整数となる。SciPy は前者。

def my_geom_pmf(x: int, p: float):
    return (1 - p) ** (x - 1) * p

実行例 ( $p = 1 / 6$ )

>>> p = 1 / 6
>>> x_list = range(1, 6)
>>> [my_geom_pmf(x, p) for x in x_list]
[0.16666666666666666,
 0.1388888888888889,
 0.11574074074074076,
 0.09645061728395063,
 0.08037551440329219]

scipy.stats.geom.pmf() の実行例 ( $p = 1 / 6$ )

>>> [stats.geom.pmf(x, p) for x in x_list]
[np.float64(0.16666666666666666),
 np.float64(0.1388888888888889),
 np.float64(0.11574074074074076),
 np.float64(0.09645061728395063),
 np.float64(0.08037551440329219)]

確率分布を描画してみると、

>>> x = range(1, 11)
>>> y = stats.geom.pmf(x, p)
>>>
>>> fig, axes = plt.subplots()
>>> axes.plot(x, y, "o")
>>> axes.vlines(x, 0, y)
>>> axes.set_xlabel("Num. of trial")
>>> axes.set_ylabel("Probability")
>>> plt.show()

単調減少する等比数列（幾何数列）となる。

幾何分布の期待値と分散

幾何分布の期待値は $1/p$ 、分散は $(1-p)/p^2$ となる。証明は幾何分布の期待値(平均)・分散・標準偏差とその導出証明 | 数学の景色などを参照。定義からの導出と特性関数からの導出が紹介されている。

>>> p = 1 / 6
stats.geom.mean(p)
np.float64(6.0)
>>> stats.geom.var(p)
np.float64(30.000000000000007)

乱数生成をして確かめてみる。

>>> size = 1000
>>> rng = np.random.default_rng(1728)
>>> data_geom = stats.geom.rvs(p, size=size, random_state=rng))
>>>
>>> fig, axes = plt.subplots()
>>> sns.histplot(data=data_geom, kde=True, ax=axes)
>>> plt.show()

>>> np.mean(data_geom)
np.float64(5.95)
>>> np.var(data_geom)
np.float64(29.7775)

共に理論値に近いことが確かめられた。

指数分布 (Exponential distribution))

P(x) = \lambda e^{- \lambda x} = \frac{1}{\theta} e^{- x / \theta}

補足: 指数分布の母数 $\lambda$ はポアソン分布と同様に「単位時間あたりの平均発生回数」を表す。一方で SciPy のように 尺度母数 $\theta (= 1 / \lambda)$ とする場合もある。尺度母数は引数 scale で指定する。

def my_expon_pdf(x: float, scale: float):
    return math.exp(-x / scale) / scale

例：1時間 に平均10回電話が鳴るコールセンター。5分後 に次に電話が鳴る確率密度は？

時間から分へ単位を変換する。

\begin{aligned} \lambda &= 10 / 60 = 1/6 \\ \theta &= 1 / \lambda = 6 \end{aligned}

実行例 ( $\theta = 6, x = 5$ )

>>> lam = 10 / 60
>>> scale = 1 / lam
>>> x = 5
>>> my_expon_pdf(x, scale)
0.0724330347511797

scipy.stats.expon.pdf() の実行例 ( $\theta = 6, x = 5$ )

>>> stats.expon.pdf(x, scale=scale)
np.float64(0.0724330347511797))]

確率分布を描画してみると、

>>> x = range(0, 41)
>>> y = stats.expon.pdf(x, scale=scale)
>>>
>>> fig, axes = plt.subplots()
>>> axes.plot(x, y)
>>> axes.set_xlabel("Time[min]")
>>> axes.set_ylabel("Probability Density")
>>> plt.show()

まさに単調減少する指数関数である。

>>> stats.expon.cdf(5, scale=scale)()
np.float64(0.5654017914929218)

5分 以内に 次に電話が鳴る確率は 56.5% となる。

指数分布の期待値と分散

指数分布の期待値は $1/\lambda = \theta$ 、分散は $1/\lambda^2 = \theta^2$ となる。証明は【徹底解説】指数分布とは | Academaid など参考にして頂きたい。確率密度関数の導出（幾何分布との関係など）も詳しい。

>>> lam = 1 / 6
>>> scale = 1 / lam
>>> stats.expon.mean(scale=scale)
np.float64(6.0)
>>> stats.expon.var(scale=scale)
np.float64(36.0)

乱数生成をして確かめてみる。

>>> size = 1000
>>> rng = np.random.default_rng(1728)
>>> data_exp = stats.expon.rvs(scale=scale, size=size, random_state=rng)))
>>>
>>> fig, axes = plt.subplots()
>>> sns.histplot(data=data_exp, kde=True, ax=axes)
>>> plt.show()

>>> np.mean(data_exp)
np.float64(5.944093619768581)
>>> np.var(data_exp)
np.float64(35.69083040632959)

共に理論値に近い。

22 正規分布の正体を暴く

正規分布の正体を暴く【統計学入門22】 - 米国データサイエンティストのブログ

正規分布
- PDF: scipy.stats.norm.pdf(x, loc=0, scale=1)
- CDF: scipy.stats.norm.cdf(x, loc=0, scale=1)
- 期待値: scipy.stats.norm.mean(loc=0, scale=1)
- 分散: scipy.stats.norm.var(loc=0, scale=1)
- 乱数生成: scipy.stats.geom.rvs(loc=0, scale=1, size=1, random_state=None)

参照：scipy.stats.norm — SciPy v1.18.0 Manual

確率密度関数 (PDF)

正規分布 $N(\mu, \sigma^2)$ の確率密度関数 (PDF) は、次の式で表される。

\frac{1}{\sqrt{2 \pi} \sigma} \exp{\bigg(- \frac{(x-\mu)^2}{2 \sigma^2} \bigg)}

補足: SciPyでは、期待値 $\mu$ を loc で、標準偏差 $\sigma$ をscale で指定する。分散 $\sigma^2$ ではないことに注意！

def my_norm_pdf(x: float, loc: float = 0, scale: float = 1):
    return np.exp(-((x - loc) ** 2) / (2 * scale**2)) / (math.sqrt(2 * math.pi) * scale)

具体例：偏差値 ( $\mu = 50, \sigma = 10$ )

実行例

>>> mu = 50
>>> sigma = 10
>>>
>>> x = np.arange(30, 71, 10)
>>> my_norm_pdf(x, loc=mu, scale=sigma)
array([0.0053991 , 0.02419707, 0.03989423, 0.02419707, 0.0053991 ])

scipy.stats.norm.pdf() の実行例

>>> stats.norm.pdf(x, loc=mu, scale=simga)
array([0.0053991 , 0.02419707, 0.03989423, 0.02419707, 0.0053991 ])

確率分布を描画してみると、

>>> x = np.arange(10, 90, 0.1)
>>> y = stats.norm.pdf(x, loc=mu, scale=sigma)
>>>
>>> fig, axes = plt.subplots()
>>> axes.plot(x, y)
>>> axes.set_xlabel("x")
>>> axes.set_ylabel("Probability Density")
>>> plt.show()

偏差値60は上位15.9%以内、偏差値70は上位2.3%以内である。

>>> 1 - stats.norm.cdf(x=60, loc=mu, scale=sigma)
np.float64(0.15865525393145707)
>>> 1 - stats.norm.cdf(x=70, loc=mu, scale=sigma)
np.float64(0.02275013194817921)

期待値と分散

正規分布 $N(\mu, \sigma^2)$ の期待値は $\mu$ 、分散は $\sigma^2$ となる。証明は正規分布の期待値(平均)・分散・標準偏差とその導出証明 | 数学の景色などを参照。定義からの導出と特性関数からの導出が紹介されている。

標準正規分布 ( $\mu = 0, \sigma = 1$ ) を考える。

>>> mu = 0
>>> sigma = 1
>>> stats.norm.mean(loc=mu, scale=sigma)
np.float64(0.0)
>>> stats.norm.var(loc=mu, scale=sigma)
np.float64(1.0)

乱数生成をして確かめてみる。

>>> size = 1000
>>> rng = np.random.default_rng(1728)
>>> data_norm = stats.norm.rvs(loc=mu, scale=sigma, size=size, random_state=rng)
>>>
>>> fig, axes = plt.subplots()
>>> sns.histplot(data=data_norm, kde=True, ax=axes)
>>> plt.show()

$k \sigma$ より大きくなる割合を見てみる。

>>> for k in [1, 2, 3]:
>>>     print(
>>>         f"Probability of greater than {k}: {stats.norm.cdf(x=-k, loc=mu, scale=sigma):.3f}"
>>>     )
>>>     print(f"Samples of greater than {k}: {np.count_nonzero(data_norm > k) / size:.3f}")
Probability of greater than 1: 0.159
Samples of greater than 1: 0.167
Probability of greater than 2: 0.023
Samples of greater than 2: 0.029
Probability of greater than 3: 0.001
Samples of greater than 3: 0.003

期待値と分散を計算する。

>>> np.mean(data_norm)
np.float64(0.012092411878496589)
>>> np.var(data_norm)
np.float64(1.064604085685525)

共に理論値に近い。

二項分布の正規近似

二項分布 $Bin(n, p)$ は $n$ が大きい時に正規分布で近似できる。

$n$ を大きくしながら確率分布を見てみると、

>>> n_list = [5, 10, 20, 50]
>>> p = 1 / 6
>>> x = np.arange(18)
>>>
>>> fig, axes = plt.subplots()
>>> for n in n_list:
>>>     axes.plot(x, stats.binom.pmf(k=x, n=n, p=p), "o--", label=f"Bin({n}, {p:.3f})")
>>> axes.legend()
>>> axes.set_xlabel("x")
>>> axes.set_ylabel("Probability")
>>> plt.show()

徐々に丸みを帯びて正規分布に似てくるように思える。

これをより正確に述べると、二項分布に従う確率変数 $X_n \sim Bin(n, p)$ を標準化した確率変数は標準正規分布に 分布収束 する。

\frac{X_n - np}{\sqrt{n p (1-p)}} \to_d N(0, 1) \quad (n \to \infty)

二項分布を標準化したものを標準正規分布と重ね合わせてみる。

>>> n_list = [3, 6, 50]
>>> p = 1 / 6
>>>
>>> # 描画する範囲（横軸）
>>> z_range = (-3, 4)
>>>
>>> fig, axes = plt.subplots()
>>> # 標準化した二項分布
>>> for n in n_list:
>>>     x = np.arange(n + 1)
>>>
>>>     mu = n * p
>>>     sigma = np.sqrt(n * p * (1 - p))
>>>     # 変数変換（標準化）
>>>     z = (x - mu) / sigma
>>>     y_z = stats.binom.pmf(k=x, n=n, p=p) * sigma
>>>     # 描画範囲を0付近に限定
>>>     z_filter = (z >= z_range[0]) & (z <= z_range[1])
>>>     axes.plot(z[z_filter], y_z[z_filter], "o--", label=f"Bin({n}, {p:.3f})")
>>>
>>> # 標準正規分布
>>> z = np.arange(z_range[0], z_range[1], 0.1)
>>> y_z = stats.norm.pdf(x=z)
>>> axes.plot(z, y_z, label=f"Norm(0, 1)")
>>>
>>> axes.legend()
>>> axes.set_xlabel("z")
>>> axes.set_ylabel("Probability")
>>> plt.show()

実際に標準正規分布に近づく事が確認できる。

ポアソン分布の正規近似

ポアソン分布 $Po(\lambda)$ も $\lambda$ が大きい時に正規分布で近似できる。

二項分布の時と同様に、 $\lambda$ を大きくしながら確率分布を見てみる。

>>> mu_list = [1, 3, 5, 10]
>>> x = np.arange(20)
>>>
>>> fig, axes = plt.subplots()
>>> for mu in mu_list:
>>>     axes.plot(x, stats.poisson.pmf(x, mu), "o--", label=f"Poisson({mu})")
>>> axes.legend()
>>> axes.set_xlabel("x")
>>> axes.set_ylabel("Probability")
>>> plt.show()

正確に述べると、ポアソン分布に従う確率変数 $X_ \sim Po(\lambda)$ を標準化した確率変数は標準正規分布に 分布収束 する。

\frac{X_n - \lambda}{\sqrt{\lambda}} \to_d N(0, 1) \quad (\lambda \to \infty)

標準正規分布と重ね合わせてみる。

>>> mu_list = [1, 3, 10]
>>>
>>> # 描画する範囲（横軸）
>>> z_range = (-3, 4)
>>>
>>> fig, axes = plt.subplots()
>>> # 標準化したポアソン分布
>>> for mu in mu_list:
>>>     x = np.arange(30)
>>>
>>>     # 変数変換（標準化）
>>>     sigma = np.sqrt(mu)
>>>     z = (x - mu) / sigma
>>>     y_z = stats.poisson.pmf(x, mu) * sigma
>>>     # 描画範囲を0付近に限定
>>>     z_filter = (z >= z_range[0]) & (z <= z_range[1])
>>>     axes.plot(z[z_filter], y_z[z_filter], "o--", label=f"Poisson({mu})")
>>>
>>> # 標準正規分布
>>> z = np.arange(z_range[0], z_range[1], 0.1)
>>> y_z = stats.norm.pdf(x=z)
>>> axes.plot(z, y_z, label=f"Norm(0, 1)")
>>>
>>> axes.legend()
>>> axes.set_xlabel("z")
>>> axes.set_ylabel("Probability")
>>> plt.show()

$\lambda = 3$ でも標準正規分布に結構近づいている。