ヒストグラムを作成して、そこに正規分布曲線も重ねたい。matplotlib編

概要

Python + pandas + matplotlib でヒストグラムを作成して、その母集団が正規分布に従うと仮定して求めた正規分布曲線を重ねてみたいと思います。

また、あわせてシャピロ・ウィルク検定（Shapiro-Wilk検定）による正規性のチェック、正規分布を仮定した母集団の平均値に対する区間推定（95%信頼区間を求める）も行なって行きます。

ここでは、例題として、100点満点の試験得点（100人分）についてヒストグラムを作成してみたいと思います。

■ 平均：76.1点、標準偏差：10.3点

• p=0.77 ( p>=0.05 ) であり母集団には正規性があると言える
• 母平均の95%信頼区間CI = [74.09 , 78.19]

関連エントリ

• 正規性のチェックについては「Pythonでt検定　2クラスの試験成績の比較」のなかで解説しています。

実行環境

Google Colab.（Python 3.6.9）で実行・動作確認をしています。ほぼ Jupyter Notebook と同じです。

!pip list 
matplotlib               3.1.2
numpy                    1.17.4
pandas                   0.25.3
scipy                    1.3.3

matplotlibで日本語を使うための準備

matplotlib の出力図のなかで、日本語が使えるようにします。

!pip install japanize-matplotlib
import japanize_matplotlib

以上により、japanize-matplotlib-1.0.5 がインストール、インポートされて、ラベル等に日本語を使っても文字化け（豆腐化）しなくなります。

コード

%reset -f
import numpy as np
import pandas as pd
import scipy.stats as st
import math                        
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker 
import matplotlib.transforms as ts 

# 試験成績のダミーデータ（0～100の範囲の正規乱数（整数））を生成　
def getTestData(mu,sig,n) : # 平均mu、標準偏差sig、正規乱数の個数n
  f = lambda x: min(100,max(0,round( np.random.normal(mu,sig),0)))
  return (np.frompyfunc(f, 1, 1)(np.zeros(n))).astype(np.float)

df = pd.DataFrame( {'p1':getTestData(75,8,40)} )

# グラフ描画パラメータ
target = 'p1'           # データフレームのなかでプロット対象とする列
x_min, x_max = 40, 100  # プロットする点数範囲（下限と上限）
j = 5                   # Y軸（度数）刻み幅
k = 5                   # 区間の幅
bins = 12               # 区間の数　(x_max-x_min)/k  (100-40)/5->12

# ここからグラフ描画処理
plt.figure(dpi=96)
plt.xlim(x_min,x_max)
d = 0.001

# (1) 統計処理
n   = len(df[target])         # 標本の大きさ 
mu  = df[target].mean()       # 平均
sig = df[target].std(ddof=0)  # 標準偏差：ddof(自由度)=0
print(f'■ 平均：{mu:.1f}点、標準偏差：{sig:.1f}点')
ci1, ci2 = (None, None)

# 正規性の検定（有意水準5%）と母平均の95%信頼区間
_, p = st.shapiro(df[target])
if p >= 0.05 :
  print(f'  - p={p:.2f} ( p>=0.05 ) であり母集団には正規性があると言える')
  U2 = df[target].var(ddof=1)  # 母集団の分散推定値（不偏分散）
  DF = n-1                     # 自由度
  SE = math.sqrt(U2/n)         # 標準誤差
  ci1,ci2 = st.t.interval( alpha=0.95, loc=mu, scale=SE, df=DF )
  print(f'  - 母平均の95%信頼区間CI = [{ci1:.2f} , {ci2:.2f}]')
else:
  print(f'  ※ p={p:.2f} ( p<0.05 ) であり母集団には正規性があるとは言えない')

# (2) ヒストグラムの描画
hist_data = plt.hist(df[target], bins=bins, color='tab:cyan', range=(x_min, x_max), rwidth=0.9)
plt.gca().set_xticks(np.arange(x_min,x_max-k+d, k))

# (3) 正規分布を仮定した近似曲線
sig = df[target].std(ddof=1)  # 不偏標準偏差：ddof(自由度)=1
nx = np.linspace(x_min, x_max+d, 150) # 150分割
ny = st.norm.pdf(nx,mu,sig) * k * len(df[target])
plt.plot( nx , ny, color='tab:blue', linewidth=1.5, linestyle='--')

# (4) X軸 目盛・ラベル設定
plt.xlabel('試験点数区分',fontsize=12)
plt.gca().set_xticks(np.arange(x_min,x_max+d, k))

# (5) Y軸 目盛・ラベル設定
y_max = max(hist_data[0].max(), st.norm.pdf(mu,mu,sig) * k * len(df[target]))
y_max = int(((y_max//j)+1)*j) # 最大度数よりも大きい j の最小倍数
plt.ylim(0,y_max)
plt.gca().set_yticks( range(0,y_max+1,j) )
plt.ylabel('人数',fontsize=12)

# (6) 平均点と標準偏差のテキスト出力
tx = 0.03 # 文字出力位置調整用
ty = 0.91 # 文字出力位置調整用
tt = 0.08 # 文字出力位置調整用
tp = dict( horizontalalignment='left',verticalalignment='bottom',
           transform=plt.gca().transAxes, fontsize=11 )
plt.text( tx, ty, f'平均点 {mu:.2f}', **tp)
plt.text( tx, ty-tt, f'標準偏差 {sig:.2f}', **tp)
plt.vlines( mu, 0, y_max, color='black', linewidth=1 )

コード2

ヒストグラムの場合、区分の境界となる値が左右どちらの度数にカウントされているのか分かりづらい場合があります。以下は、それを分かりやすくしたバージョンです。

%reset -f
import numpy as np
import pandas as pd
import scipy.stats as st
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import matplotlib.transforms as ts
import math

# テスト用ダミーデータ（0～100）生成　
def getTestData(mu,sig,n) : # 平均mu、標準偏差sig、正規乱数の個数n
  f = lambda x: min(100,max(0,round( np.random.normal(mu,sig),0)))
  return (np.frompyfunc(f, 1, 1)(np.zeros(n))).astype(np.float)

df = pd.DataFrame( {'p1':dm} )

# グラフ描画パラメータ
target = 'p1'           # データフレームのなかでプロット対象とする列
x_min, x_max = 40, 100  # プロットする点数範囲（下限と上限）
j = 5                   # Y軸（度数）刻み幅
k = 5                   # 区間の幅
bins = 12               # 区間の数　(x_max-x_min)/k  (100-40)/5->12

# ここからグラフ描画処理
plt.figure(dpi=96)
plt.xlim(x_min-k/2,x_max-k/2)
d = 0.001

# (1) 統計処理
n   = len(df[target])         # 標本の大きさ 
mu  = df[target].mean()       # 平均
sig = df[target].std(ddof=0)  # 標準偏差：ddof(自由度)=0
print(f'■ 平均：{mu:.1f}点、標準偏差：{sig:.1f}点')
ci1, ci2 = (None, None)

# 正規性の検定（有意水準5%）と母平均の95%信頼区間
_, p = st.shapiro(df[target])
if p >= 0.05 :
  print(f'  - p={p:.2f} ( p>=0.05 ) であり母集団には正規性があると言える')
  U2 = df[target].var(ddof=1)  # 母集団の分散推定値（不偏分散）
  DF = n-1                     # 自由度
  SE = math.sqrt(U2/n)         # 標準誤差
  ci1,ci2 = st.t.interval( alpha=0.95, loc=mu, scale=SE, df=DF )
  print(f'  - 母平均の95%信頼区間CI = [{ci1:.2f} , {ci2:.2f}]')
else:
  print(f'  ※ p={p:.2f} ( p<0.05 ) であり母集団には正規性があるとは言えない')

# (2) ヒストグラムの描画
hist_data = plt.hist(df[target], bins=bins, color='tab:cyan', range=(x_min, x_max), rwidth=0.9, align='left')
plt.gca().set_xticks(np.arange(x_min,x_max-k+d, k))

# (3) 正規分布を仮定した近似曲線
sig = df[target].std(ddof=1)  # 不偏標準偏差：ddof(自由度)=1
nx = np.linspace(x_min, x_max+d, 150) # 150分割
ny = st.norm.pdf(nx,mu,sig) * k * len(df[target])
plt.plot( nx - k/2 , ny, color='tab:blue', linewidth=1.5, linestyle='--')

# (4) X軸 目盛・ラベル設定
plt.xticks(rotation=90)
plt.xlabel('試験点数区分',fontsize=12)
tmp = list([ f'${i} \\leq x < {i+k}$' for i in range(x_min,x_max-k+1, k) ]) 
tmp[-1] = f'${x_max-k} \\leq x \\leq {x_max}$'
plt.gca().set_xticklabels( tmp )

# (5) Y軸 目盛・ラベル設定
y_max = max(hist_data[0].max(), st.norm.pdf(mu,mu,sig) * k * len(df[target]))
y_max = int(((y_max//j)+1)*j) # 最大度数よりも大きい j の最小倍数
plt.ylim(0,y_max)
plt.gca().set_yticks( range(0,y_max+1,j) )
plt.ylabel('人数',fontsize=12)