コンテンツにスキップ

補足説明 分位数#

分位数に関連して本篇のなかで言及・引用した公式や定理について補足する記事です。

分位数とは#

確率変数がとりうる値の区間（データセット）を一定の考え方で分割したとき、その分割位置の値を分位数（Quantile）と呼ぶ．よく使われる分割方法には次のものがある．   

  • 中央値：
    データセットの中央値（Median）を分位数と決める。データセットの要素数Nが偶数の場合には中央の値がデータセット内に無いので等分割線の上下2個の値の平均値が中央値になることに注意する．

  • パーセンタイル：
    データセットを昇順に並べて100等分してできる分位数をNパーセンタイル（N-percentile）と呼ぶ．例えば25パーセンタイルはデータセットの 下から25% 以下の境界を指す．

四分位数#

NパーセンタイルのNとして25,50,75をとり全データを4区域に分類したときの各区切り線を四分位数ともいうそれぞれ第1四分位数，第2四分位数，第3四分位数とも呼ぶ．簡単なデータセットを例として以下に図解する．この図解にある四分位数の決め方は高等学校の指導要領に記載されている方法のようだ．

quantile def1

四分位数はデータの外れ値を検出したいときにも利用することができる．例えば以下のようにIQRの1.5倍幅から外れるデータを外れ値とする．

data < Q1 - 1.5 * IQR
data > Q3 + 1.5 * IQR

matplotlibの箱ひげ図（boxplot関数）を使うと，パーセンタイルと外れ値を可視化することができる．
boxplotは，データの第 1 四分位 (Q1) から第 3 四分位 (Q3) まで広がり、中央値に線が引かれた箱（box）を描く．同時にボックスから四分位範囲 (IQR) の 1.5 倍以内にある最も遠いデータポイントまで伸びたひげ（whisker）を描く．さらにヒゲの先端を越えたデータポイントを外れ値（outliner）として描く。以下にサンプルコードと実行した結果を記載する．

# reg_plot_boxplot.py
# docs: https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.boxplot.html

import matplotlib.pyplot as plt
import numpy as np
import japanize_matplotlib

np.random.seed(42)

# サンプルデータセットを３グループ生成する
data = [np.random.normal(loc=i*5, scale=1, size=100) for i in range(3)]

# ２個目のグループに外れ値のデータを追加してみる
outliers = np.random.normal(loc=5, scale=2, size=5)
data[1] = np.concatenate([data[1], outliers])

#箱ひげ図をつくる
# notch Trueにすると中央値の信頼性区間を示す窪みを描画する
# vert Trueにすると縦方向に描画する
plt.boxplot(data, labels=['Group 1', 'Group 2', 'Group 3'], notch=False, vert=False)
plt.title('箱ひげ図')
plt.xlabel('値')
plt.ylabel('グループ')
plt.grid(True)
plt.tight_layout()
plt.show()
quantile def1

numpyとexcelのパーセンタイル計算の差異#

nをデータ配列の要素数、pをパーセント100分率としたときパーセンタイル位置は、前出の高校指導要領の定義はともかくとして

PercentileLocation=(n+1)∗(P/100) \begin{equation} \begin{aligned} PercentileLocation = (n + 1) * (P / 100) \end{aligned} \end{equation}

という定義をどこかの資料で読んで，これが定義かと単純に思い込んでいた． これに従うと8要素のデータ[1,2,3,4,5,6,7,100]に対して25パーセンタイル位置は (8+1)*25/100 = 2.25

実際excelのPERCENTILE関数は(Pが10未満あるいは90超の場合以外は）この値を返す．しかしながらそうはならない定義もある．numpyにpercentile関数があるので計算してみると、np.percentile(data, 25)は「2.75」を返す．

この理由はnumpy.percentile関数の既定値パラメータ（linear）では，計算式が式(1)ではなく、以下を採用しているからである．

PercentileLocation=(n−1)∗(P/100)+1 \begin{equation} \begin{aligned} PercentileLocation = (n - 1) * (P / 100) + 1 \end{aligned} \end{equation}

式(2)を使うとnumpy.percentileが返す値が得られる．

(8-1)*25/100 + 1 = 2.75

式(1)の計算方法はexcelのPERCENTILE.EXC関数と同じで、式(２)の計算方法はexcelのPERCENTILE.INC関数と同じである．ただしPERCENTILE.EXC関数ではパーセンテージの指定が1/(N+1)以下およびにN/(N+1)以上では＃NUM!を返す．
データセットの値が小さなレンジと大きなレンジで式(1)を使ってエラーが返るようなケースでは式(2)を用いることでエラーを回避できる．ちなみにexcelのPERCENTILE.INC関数はnumpy.percentileと同じ値を返す．

numpyやPERCENTILE.INCでのパーセンタイル定義の図解を以下に示す．前出の定義と比べてみれば相違がわかると思う．

quantile def2

パーセンタイルの定義方法がたくさんできた経緯は良く知らないが，データセット数が大きければどの定義でも差異が小さくなるので一貫した使い方をしていればOKということだろう．

正規分布の四分位数#

標本が標準正規分布に従う場合の四分位数の位置を以下の図で示す．

quantile_of_normal_dist

正規分布の四分位位置は例えばexcelのnorm.inv関数を使うことで簡単に得ることができる．

quantile_of_normal_inv_dist

QQプロット#

Q-Q (分位数-分位数) プロットは、指定したデータセットが特定の理論的分布（典型的には正規分布）に従っているかどうかを評価するためのツールである．理論的分布の分位数を横軸に、残差の分布の分位数を縦軸とした散布図を描く． 以下に簡単なイメージを示す．分位数の位置の比較を図示することで観測データと理論的な分布関数とのギャップを視覚的に理解できる。

quantile_qqplot_image

statsmodelsのqqplot関数を例として、プロットの特徴を説明する．
API参照先：statsmodels.api 

statsmodels.graphics.gofplots.qqplot(
    data, 
    dist=<scipy.stats._continuous_distns.norm_gen object>, 
    distargs=(), 
    a=0, 
    loc=0, 
    scale=1, 
    fit=False, 
    line=None, 
    ax=None, 
    **plotkwargs)

dataは標本データを格納した一次元配列．
distは比較の元になる理論的分布である．省略すると標準正規分布が使用される．
distargs,a,loc,scaleはfitがFalseの場合にdistへ渡されるオプションである．
fitがFalseで，さらにlineを"45"にセットすると標本データの分布が理論分布に一致しているかどうかをプロットで見ることができる．
fitがTrueの場合，理論分布のパラメータは推定値がセットされる．
lineパラメータの各オプションは以下 - "45"
45度の参照用直線をプロットする． - "s"
ｘを理論分布の分位数の配列としたとき，標本データの標準偏差m，平均bを算出し，y = x・m + bを参照用直線としてプロットする．標本データが平均b，標準偏差mの理論分布に従っていれば散布図の各点はこの直線上にプロットされることになる． - "r"
標本データから作成した回帰直線を参照用直線としてプロットする。 - "q" 第1四分位数と第3四分位数を通る参照用直線をプロットする。

以下にstatsmodelsのサンプルコード（微修正版）を記載する． 

# reg_check_plot_longley.py
# slitely modified version of sample program by statsmodels docs. 
#
from matplotlib import pyplot as plt
import statsmodels.api as sm

#example from docstring
data = sm.datasets.longley.load()
data.exog = sm.add_constant(data.exog, prepend=True)
mod_fit = sm.OLS(data.endog, data.exog).fit()
res = mod_fit.resid
left = -1.8   #x coordinate for text insert
fig = plt.figure()

# line='45'
ax = fig.add_subplot(2, 2, 1)
sm.graphics.qqplot(res, line='45', fit=True, ax=ax)
ax.set_xlim(-2, 2)
top = ax.get_ylim()[1] * 0.75
txt = ax.text(left, top, "line='45', \nfit=True", verticalalignment='top')
txt.set_bbox(dict(facecolor='k', alpha=0.1))

# line='s'
ax = fig.add_subplot(2, 2, 2)
sm.graphics.qqplot(res, line='s', fit=False,ax=ax)
top = ax.get_ylim()[1] * 0.75
txt = ax.text(left, top, "line='s', \nfit=False", verticalalignment='top')
txt.set_bbox(dict(facecolor='k', alpha=0.1))


# line='q'
ax = fig.add_subplot(2, 2, 3)
sm.graphics.qqplot(res, line='q', fit=True,ax=ax)
top = ax.get_ylim()[1] * 0.75
txt = ax.text(left, top, "line='q', \nfit=True", verticalalignment='top')
txt.set_bbox(dict(facecolor='k', alpha=0.1))

# line='r'
ax = fig.add_subplot(2, 2, 4)
sm.graphics.qqplot(res, line='r', fit=True,ax=ax)
top = ax.get_ylim()[1] * 0.75
txt = ax.text(left, top, "line='r', \nfit=True", verticalalignment='top')
txt.set_bbox(dict(facecolor='k', alpha=0.1))

fig.tight_layout()
plt.show()

結果のプロット

P-value