1. code0327

    Posted

    code0327
Changes in title
+ヒストグラムを作成して、そこに正規分布曲線も重ねたい。matplotlib編
Changes in tags
Changes in body
Source | HTML | Preview
@@ -0,0 +1,204 @@
+# 概要
+
+Python + pandas + matplotlib で**ヒストグラム**を作成して、その母集団が正規分布に従うと仮定して求めた**正規分布曲線**を重ねてみたいと思います。
+
+また、あわせて**シャピロ・ウィルク検定**(Shapiro-Wilk検定)による**正規性のチェック**、正規分布を仮定した**母集団の平均値に対する区間推定**(95%信頼区間を求める)も行なって行きます。
+
+ここでは、例題として、100点満点の試験得点(100人分)についてヒストグラムを作成してみたいと思います。
+
+![fig1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152805/3b44c596-e767-0b74-785a-5b179f923517.png)
+
+
+>■ 平均:76.1点、標準偏差:10.3点
+> - p=0.77 ( p>=0.05 ) であり母集団には正規性があると言える
+> - 母平均の95%信頼区間CI = [74.09 , 78.19]
+
+## 関連エントリ
+
+- 正規性のチェックについては「[Pythonでt検定 2クラスの試験成績の比較](https://qiita.com/code0327/items/a96dd2fbd8a491d2eeaa)」のなかで解説しています。
+
+
+# 実行環境
+
+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` がインストール、インポートされて、ラベル等に日本語を使っても文字化け(豆腐化)しなくなります。
+
+# コード
+
+```python:
+%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
+
+ヒストグラムの場合、区分の境界となる値が左右どちらの度数にカウントされているのか分かりづらい場合があります。以下は、それを分かりやすくしたバージョンです。
+
+![fig2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152805/b5ad6124-2f4d-64d4-e608-ab306031541e.png)
+
+
+```python:
+%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)
+```