【python】データが正規分布に従うかどうかの確認(中編)

はじめに

前回の記事でデータが正規分布に従っているかどうか判断する方法について説明しました。
今回はその中の一つであるQ-Qプロットについて説明します。

はじめに
データの正規性の確認方法
Q-Qプロットで確認
- scipyを使ったQ-Qプロットの作成
- Q-Qプロットの作成方法
終わりに
- pythonと統計学の勉強について

データの正規性の確認方法

以下の方法が主に使われます。

ヒストグラムを目視
歪度、尖度を確認
Q-Qプロットで確認
シャピロウィルク検定
コルモゴロフスミルノフ検定

今回はQ-Qプロットについて説明します。

サンプルデータは前回同様下記の4つのデータを使います。

①正規分布に従うデータ
②ログ正規分布に従うデータ
③正規分布に従う2つの母集団からサンプリングしたデータ
④一様分布に従うデータ

データは以下コードで準備します、ヒストグラムも作っておきます。

import numpy as np
from scipy.stats import norm, lognorm, uniform
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

# ①正規分布に従うサンプル
data1 = norm.rvs(loc=7.5, scale=1, size=100)

# ②ログ正規分布に従うサンプル
data2 = lognorm.rvs(0.95, loc=5, scale=1, size=100)

# ③2つの母集団からのサンプル
data31 = norm.rvs(loc=10, scale=1, size=70)
data32 = norm.rvs(loc=6, scale=1.5, size=30)
data3 = np.concatenate([data31, data32])

# ④一様分布に従うサンプル
data4 = uniform.rvs(loc=5,scale=8,size=100)
fig, ax = plt.subplots(2, 2)
ax[0, 0].hist(data1, density=True)
ax[0, 0].set_xlabel('x')
ax[0, 0].set_ylabel('density')
ax[0, 0].set_title('1. Normal')
ax[0, 0].set_xlim(0, 15)

ax[0, 1].hist(data2, density=True)
ax[0, 1].set_xlabel('x')
ax[0, 1].set_ylabel('density')
ax[0, 1].set_title('2. Log Normal')
ax[0, 1].set_xlim(0, 15)

ax[1, 0].hist(data3, density=True)
ax[1, 0].set_xlabel('x')
ax[1, 0].set_ylabel('density')
ax[1, 0].set_title('3. Sample from two population')
ax[1, 0].set_xlim(0, 15)

ax[1, 1].hist(data4, density=True)
ax[1, 1].set_xlabel('x')
ax[1, 1].set_ylabel('density')
ax[1, 1].set_title('4. Uniform')
ax[1, 1].set_xlim(0, 15)

plt.tight_layout()
plt.show()

結果は以下の通り。

Q-Qプロットで確認

QQプロットとは、2つの確率分布の分位数を互いにプロットして比較する方法のことです。
正規分布かどうか判定するには、正規分布の分位数と、サンプルデータの分位数をプロットします。プロットした点がほぼ直線上に位置していれば正規分布であると言えます。

scipyを使ったQ-Qプロットの作成

scipy.statsに関数が用意されてますのでそれを使用すれば1行でできます。
やってみます。

from scipy.stats import probplot
fig, ax = plt.subplots(2, 2)
probplot(data1, plot=ax[0, 0])
probplot(data2, plot=ax[0, 1])
probplot(data3, plot=ax[1, 0])
probplot(data4, plot=ax[1, 1])
plt.tight_layout()
plt.show()

結果は以下。

正規分布に従うサンプルデータである①は直線上にデータが乗っておりますが、それ以外のデータでは直線から離れたデータが散見されますね。
ヒストグラムと比べて少々分かりやすいというのが利点です。

とりあえず1番簡単な手法で実装しましたが、何が起きているのか、x軸、y軸の値はなんなのか等不明点が多いと思いますので、自分で一から作成してみます。

Q-Qプロットの作成方法

上記の関数で書かれるQ-Qプロットは標準正規分布(平均0, 標準偏差1)とサンプルの分布の比較を行っております。
作成手順は以下の通り。

サンプルデータを昇順にソート
サンプルデータの分位数を求める
サンプルデータの分位数に対応する標準正規分布の分位数を求める

①の正規分布に従うデータを使って上記を順番に行ってみましょう

1. サンプルデータを昇順にソート

これは簡単ですね

sorted_data1 = sorted(data1)

2. サンプルデータの分位数を求める

サンプルデータの累積確率密度関数を求めるのと同義です。
サンプルデータの場合、真の分布がわからないので、昇順にしたサンプルデータの $k$ 番目 $(k = 1, 2, ... ,n)$ のデータを $(k - 0.5) / n$ 分位数(累積確率密度が $(k - 0.5) / n$ の時の期待値)として使用します。
分位数の選択はいろいろ種類がありますが今回は上述の方法を使用します。
上記の手法で、サンプルデータの累積確率密度関数をプロットしてみます。

cdf = [(i + 0.5) / len(sorted_data1) for i in range(len(data1))]
fig, ax = plt.subplots()
ax.scatter(sorted_data1, cdf)
ax.set_xlabel('data')
ax.set_ylabel('cdf')
plt.show()

3. サンプルデータの分位数に対応する標準正規分布の分位数を求める

標準正規分布の累積密度関数の逆関数を使います。
累積密度関数の逆関数はscipy.statsのnorm.ppfで実装されています。
norm.ppfの引数に2で求めたサンプルデータの累積確率密度のリストを入れることでサンプルデータの分位数に対応する標準正規分布の分位数を求めることができます。

norm_x = norm.ppf(cdf)
fig, ax = plt.subplots()
ax.scatter(norm_x, sorted_data1)
ax.set_xlabel('Theoretical quantiles')
ax.set_ylabel('Orderd Values')
plt.show()

大体同じような形になったのではないでしょうか？
最後に、ここに最小二乗法で近似した直線を引いて終わりです。
最小二乗法での近似は以下を参照願います。

# 最小二乗近似
A = np.vstack((norm_x, np.ones(len(norm_x)))).T
v = np.array(sorted_data1).T
u, residuals, rank, s = np.linalg.lstsq(A, v, rcond=None)
x = np.linspace(-3, 3, 100)

fig, ax = plt.subplots()
ax.scatter(norm_x, sorted_data1)
ax.plot(x, x * u[0] + u[1], 'r')
ax.set_xlabel('Theoretical quantiles')
ax.set_ylabel('Orderd Values')
plt.show()

ここでは、標準正規分布(平均0, 標準偏差1)との比較を行いましたが、サンプルデータによっては直線の傾きが大きかったり、小さかったりでお互いの比較が難しいです。
そこで、サンプルデータで平均や標準偏差を求めて、その値を使った正規分布の分位数と比較してみます。
こうすることで、どんなデータにおいても $y=x$ の直線上にあるかどうかで正規分布に従うかどうか判定できます。

norm_x = norm.ppf(cdf, loc=data1.mean(), scale=data1.std())
x = [sorted_data1[0], sorted_data1[-1]]

fig, ax = plt.subplots()
ax.scatter(norm_x, sorted_data1)
ax.plot(x, x, 'r')
ax.set_xlabel('Theoretical quantiles')
ax.set_ylabel('Orderd Values')
ax.set_xlim(4, 11)
ax.set_ylim(4, 11)
ax.set_aspect('equal')
plt.show()