Pythonで連続値の確率分布、正規分布から確率を計算してみる

正規分布は、実生活の多くの現象をモデル化するために使われる統計分布の一つです。

本記事では、正規分布を用いて、ある養鶏場で鶏卵の重さに関する確率を求める方法を紹介します。

具体的には、卵の平均重さと標準偏差から、その重さが70gを超える確率を計算します。

この過程でPythonを用いて確率密度関数（PDF）や累積分布関数（CDF）を活用し、実際のデータに基づいた計算を行います。

正規分布から確率を求める

以下の、演習で身につける統計学入門からの例題用いて、正規分布から確率を求めてみましょう。

ある養鶏場の鶏卵の平均の重さは65g、標準偏差6gの正規分布を示し、さらにランダムに１個取り出した時その重さが70gを超える確率を求める。

Pythonを使って正規分布から確率を求める

ではさっそくPythonを使って実装していきましょう。

まず正規分布の確率密度関数を、Pythonに落とし込みます。

import numpy as np
import pandas as pd
def normal_dist_pdf(x, mu, sigma):
    # ガウスの公式
    y = (1/np.sqrt(2* np.pi*sigma**2)) * np.exp(-1 * (x-mu)**2/(2*sigma**2))
    return y

正規分布を描いてみる

描画ライブラリを使って、正規分布を描いてみます。

mu = 65 # 卵の平均
sigma = 6 # 標準偏差
import matplotlib.pyplot as plt # 描画ライブラリをインポート
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
x = np.arange(40, 90) # x の範囲の設定
y = normal_dist_pdf(x, mu, sigma) # y の確率密度関数
ax.plot(x, y, color = 'b', linewidth = 5.0)
plt.grid()

やはり図で表すと、データがどのような形をとっているのか分かりやすいですね。

ここから、統計計算ライブラリscipyを使って70gを超える確率を求めてみます。

scipyからインポートするのはnormモジュール。

その、normモジュールから累積密度関数 (CDF) を使って確率を求めます。

ちなみにnorm.pdf と norm.cdfの違いを簡単に説明します。

norm.pdf と norm.cdf

まずnorm.pdf と norm.cdfの違いを簡単に説明します。

pdf = Probability Density Function (確率密度関数)

loc = 平均
scale = 標準偏差
x = -∞から∞の間の値
返り値：x の値が発生する確率(%)

cdf = Cumjulative Distribution Function (累積分布関数）

loc = 平均
scale = 標準偏差
x = -∞から∞の間の値
返り値：x 以下の値が発生する確率(%)

ということを踏まえて、卵の重さが70gを超える確率を実装してみましょう。

累積分布関数と確率密度関数のどちらを使うかという問題ですが、累積分布関数はある一点までの確率の和を知ることができ、確率密度関数はある一点の確率を知ることができます。

この場合、70gというピンポイントではなく70g以上という範囲の確率を求めるので、累積分布関数を使います。

from scipy.stats import norm # 統計ライブラリscipyからnormモジュールをインポート
p = norm.cdf(x = 70, loc = 65, scale = 6) # 累積分布関数 cdf (Cumulative Distribution Function)を使う
print(1-p)
0.20232838096364314

計算してみた結果、卵が70gを超える確率は、20%になります。

ではnormを使わずに、実装してみます。

ちなみに累積密度関数の公式は以下。

$f(x) = \frac{1}{2}\left [ 1 + erf\left ( \frac{x - \mu }{\sigma \sqrt{2}} \right ) \right ]$

この式を実装すると、以下の通り。

from scipy.special import erf
import math
def normal_dist_cdf(x, mu, sigma):
    return 1/2 * (1 + erf((x-mu)/(sigma * math.sqrt(2))))

これを、実際に数値を当てはめて計算すると、結果は以下。

import pandas as pd
print(1-(pd.DataFrame([[normal_dist_cdf(70, mu, sigma)]], index = ['卵'], columns=['70g 以上'])))

     70g 以上
卵  0.202328

結果はscipy.norm を使った時と同じ20% になりました。

以下は確率の範囲を描画したプロットになります。

以下コードになります。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['font.family'] = 'MS Gothic' #豆腐化防止
from scipy.stats import norm

# 平均と標準偏差
mu = 65 # 平均
sigma = 6 # 標準偏差

# x軸の範囲を設定
x = np.linspace(mu - 4 * sigma, mu + 4 * sigma, 1000) # 全体の範囲
y = norm.pdf(x, mu, sigma) # 確率密度関数

# 70g以上の領域を計算
x_fill = np.linspace(70, mu + 4 * sigma, 500) # 70g以上の範囲
y_fill = norm.pdf(x_fill, mu, sigma) # その範囲の確率密度

# プロット
plt.figure(figsize=(10, 6))
plt.plot(x, y, label='正規分布', color='blue') # 全体の正規分布
plt.fill_between(x_fill, y_fill, color='red', alpha=0.5, label='70g以上の確率') # 塗りつぶし
plt.axvline(70, color='red', linestyle='--', label='70g') # 70gのライン

# グラフの設定
plt.title('正規分布: 平均=65g, 標準偏差=6g', fontsize=14)
plt.xlabel('重さ (g)', fontsize=12)
plt.ylabel('確率密度', fontsize=12)
plt.legend()
plt.grid()

# グラフ表示
plt.show()

Excelを使って正規分布から確率を求める

Excelを使っても確率を求めることができます。

卵の例題の場合、NORM.DIST()関数を使い累積分布関数を求める際にはTRUE、確率密度関数を求める際にはFALSEを用いて計算します。

以下Excelを用いた計算式です。

Excelを使っても、確率は20.23%になりますね。

練習問題

では練習問題を少しやってみましょう。

問題は、演習で身につける統計学入門からの出題です。

製品Gの重さの平均が1200gで、標準偏差が11gの正規分布の場合、ランダムに1個取り出した時の重さが1180g以下である確率を求める。

ではこの正規分布を描いてみましょう。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# 平均と標準偏差
mu = 1200 # 平均
sigma = 11 # 標準偏差

# x軸の範囲を設定
x = np.linspace(1150, 1250, 500) # 滑らかな範囲を生成
y = norm.pdf(x, loc=mu, scale=sigma) # 正規分布のPDFを計算

# グラフの描画
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(x, y, color='b', linewidth=2.0) # 正規分布の描画
plt.grid()
plt.title('正規分布', fontsize=14)
plt.xlabel('値', fontsize=12)
plt.ylabel('確率密度', fontsize=12)

# グラフ表示
plt.show()

書いてみると1200gを中心にした正規分布が書けます。

ではscipy.normを使って、1180g以下になる確率を求めてみます。

from scipy.stats import norm

p = norm.cdf(x = 1180, loc = 1200, scale = 11)
print(p)
0.03451817399720763

結果の確率は、3%ともめられました。

もう一つ練習問題を解いてみましょう。

コインを200回トスし、表が出る回数の期待値（平均）と標準偏差を求め、表が出る回数が120回以上になる確率を求める。

この問題は二項分布を正規分布の両方を使います。

つまり確率変数Xが各試行で起こるpの二項分布に従うとき、施行をn回行うとXの期待値（平均）と分散は以下になります。

$E\left [ X \right ] = np$ 平均を求める

$V\left [ X \right ] = np(1-p)$ 分散を求める

では問いを四則演算だけで計算してみましょう。

200/2 # np = 200/2: 平均の計算
100.0

(200/2)*(1-(1/2)) # np x (1-1/2): 分散の計算
50.0 

50**0.5 # 50のルートをとる
7.0710678118654755

この時点で、平均と標準偏差が求められたので確率も求められるので、まずは正規分布を描いてみましょう。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# 最初の平均と分散の計算
mu = 200 / 2 # 平均を代入する
var = mu * (1 - (1 / 2)) # 分散を求める
sigma = var**0.5 # 標準偏差を求める

# 正規分布の確率密度関数を計算するためのxの範囲を指定
x = np.arange(mu - 4 * sigma, mu + 4 * sigma, 0.1)

# 正規分布の確率密度関数を計算
y = norm.pdf(x, mu, sigma)

# プロットを作成
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(x, y, color='b')

# グリッドを表示
plt.grid()

# プロットを表示
plt.show()

ではここから、scipyを使って表が出る回数が120回以上になる確率を求めます。

from scipy.stats import norm

p = norm.cdf(x = 120, loc = 100, scale = 7)
print(1-p)
0.002137366980086264

計算によると、表が出る回数が120回以上になる確率は、0.2% という結果になりました。

まとめ

本記事では、正規分布を用いて特定の確率を求める方法を学びました。

Pythonを使用して、確率密度関数や累積分布関数を活用し、鶏卵の重さが70gを超える確率を20%と計算しました。

また、正規分布の視覚的な理解を深めるためにグラフを描画し、計算結果を確認しました。

このような手法は、日常的なデータ分析や統計学の応用に役立ちます。