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

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

まず正規分布で、データの形を見てみます。

正規分布から確率を求めるための例題は、

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

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

まず正規分布の確率密度関数を、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

正規分布の確率密度関数は以下の記事に載せてますので、興味あれば見てください。

404 NOT FOUND | 数学とか統計とか

正規分布を描いてみる

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

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% になりました。

練習問題

今回も、演習で身につける統計学入門の問題を解いていきますね。

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

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

mu = 1200 # 平均
sigma = 11 # 標準偏差
import matplotlib.pyplot as plt # 描画ライブラリのインポート
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
x = np.arange(1150, 1250)
y = normal_dist_pdf(x, mu, sigma)
ax.plot(x, y, color = 'b', linewidth = 5.0)
plt.grid()

書いてみると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

この時点で、平均と標準偏差が求められたので確率も求められるのですが、もう少しプログラミングらしく書いてみましょう。

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

100.0
50.0
7.0710678118654755

めでたく平均と分散、そして標準偏差が求められたので、正規分布を描いてみます。

import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
x = np.arange(0, 200)
y = normal_dist_pdf(x, mu, sigma)
ax.plot(x, y, color = 'b', linewidth = 5.0)
plt.grid()

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

from scipy.stats import norm

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

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

プログラミング無料体験はこちら↓↓↓

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

正規分布を描いてみる

norm.pdf と norm.cdf

練習問題

関連記事:

コメント