非正規分布をrobust Z-scoreで標準化する

正規化・標準化には、min-max normalizationとZ-score Normalization(Standardization)がよく使われています。今回、robust Z-scoreを試して、上記の正規化と比較してみました。

min-max normalization

min-max normalizationとは、データが最小値0・最大値1になるようにする方法で、以下の式で正規化します。

x' = \frac{x-min(x)}{max(x)-min(x)}

pythonでは、sklearn.preprocessingのminmax_scaleかMinMaxScalerで計算できます。
この正規化は、データの分布が一様分布であることを前提としています。

Z-score Normalization(Standardization)

Z-score Normalizationとは、データが平均0・分散1になるようにする方法で、以下の式で正規化します。この値をZ-scoreといいます。μ は平均、σ は標準偏差を表しています。

x' = \frac{x-\mu}{\sigma}

pythonでは、sklearn.preprocessingのscaleかStandardScalerで計算できます。
この正規化は、データの分布が正規分布であることを前提としています。

一様分布でも正規分布でもないときはどうするのか？

実際のデータでは一様分布でも正規分布でもないことが多かったので、どうしたらよいのか調べていたところ、以下の記事でrobust Z-scoreを見つけました。

ロバストzスコア：中央値と四分位数で，非正規分布，外れ値を含む標準化
（メモ）ロバストｚスコアを用いた異常値の除外

以下で、Pythonで試してみました。

robust Z-scoreの実装

robust Z-scoreについての詳細は上記の記事を読んでいただければと思います。以下では、簡略に説明して実装します。

Z-scoreは正規分布を前提としていますが、これを非正規分布にあてはめるために、まず、平均 μを中央値 、標準偏差 σ を四分位範囲（IQR）に置き換えます。

x' = \frac{x-median(x)}{IQR}

この式はsklearn.preprocessingのrobust_scaleかRobustScalerで計算できます。

さらに、これを標準正規分布に対応できるものにします。IQRを標準正規分布に対応させたものを正規化四分位範囲（NIQR）といい、IQRをF(0.75) - F(0.25) = 1.3489で割ったものとなります。(F(x)は累積分布関数の逆関数)

NIQR = \frac{IQR}{1.3489}

先ほどの式の分母をIQRからNIQRに置き換えたものが、robust Z-scoreになります。

robust Z score = \frac{x-median(x)}{NIQR}

以上のことを踏まえて実装すると、以下のような関数となります。

def robust_z(x):
    from sklearn.preprocessing import robust_scale
    from scipy.stats import norm

    coefficient = norm.ppf(0.75)-norm.ppf(0.25)
    robust_z_score = robust_scale(x)*coefficient

    return robust_z_score

3つの正規化の比較

これまでに出てきた3つの正規化を比較したいと思います。
まず、データを用意します。一様分布でも正規分布でもないデータがほしいので、一様分布と正規分布を結合したデータを用意しました。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.stats import chisquare, shapiro, norm
from sklearn.preprocessing import minmax_scale, scale, robust_scale

np.random.seed(2020)

# 一様分布と正規分布を結合したデータ。
data = np.concatenate((np.random.uniform(low=5.0, high=10.0, size=100),
                       np.random.normal(loc=5.0, scale=1.0, size=100)))

# ヒストグラムを作画。
fig, axes = plt.subplots()
axes.hist(data)
axes.set_title("Histogram of data")
fig.show()

このデータが一様分布でないことと正規分布でないことを検定で確かめます。
一様性はカイ２乗検定、正規性は（シャピロ–ウィルク検定）で確かめました。

# 度数分布を計算する。
hist_data, _ = np.histogram(data, bins="auto")

# 一様性検定（カイ２乗検定）
_, chisquare_p = chisquare(hist_data)
print("一様性検定（カイ２乗検定）のp値 : {}".format(chisquare_p))

# 正規性検定（シャピロ–ウィルク検定）
_, shapiro_p = shapiro(data)
print("正規性検定（シャピロ–ウィルク検定）のp値 : {}".format(shapiro_p))

結果は以下の通りです。どちらもP値が0.05より小さいので、一様分布でも正規分布でもないといえそうです。

一様性検定（カイ２乗検定）のp値 : 3.8086163670115985e-09
正規性検定（シャピロ–ウィルク検定）のp値 : 8.850588528730441e-06

このデータでmin-max normalization、Z-score、robust Z-scoreを計算し、比較してみます。

# 各方法で正規化を行い、データフレームに入れる。
score_df = pd.DataFrame(data=np.array([minmax_scale(data), scale(data), robust_z(data)]).T,
                        columns=["min-max", "Z-score", "robust Z-score"])


# グラフを作成
fig, axs = plt.subplots(ncols=3, constrained_layout=True)

# x軸の幅の設定
xrange = {"min-max":(0,1),
          "Z-score":(-2.5,2.5),
          "robust Z-score":(-2.5,2.5)}

# 各ヒストグラムの作画
for i, score_name in enumerate(score_df.columns):

    axs[i].hist(score_df[score_name])
    axs[i].set_title(score_name)
    axs[i].set_xlim(xrange[score_name])

fig.show()

結果が以下の図です。あまり違いがありません。データの分布次第では違いが出てくるのかもしれません。

外れ値があるデータで比較する

そもそもrobust Z-scoreの"robust"は、外れ値に対してrobust（頑健）であるという意味です。robust Z-scoreは外れ値検出にも利用されています。
そこで、データに外れ値を入れて比較したいと思います。比較しやすいように、極端で多数の外れ値を入れてみます。

# 外れ値（一様分布）をデータに結合する。
outier = np.concatenate((data,
                         np.random.uniform(low=19.0, high=20.0, size=15)))

# 各方法で正規化を行い、データフレームに入れる。
outlier_df = pd.DataFrame(data=np.array([minmax_scale(outier), scale(outier), robust_z(outier)]).T,
                          columns=["min-max", "Z-score", "robust Z-score"])

# 外れ値なしと外れ値ありのデータフレームを結合。
concat_df = pd.concat([score_df, outlier_df],
               axis=1,
               keys=['without outlier', 'with outlier'])


# グラフを作成
fig, axs = plt.subplots(nrows=2, ncols=3, constrained_layout=True)

# x軸の幅の設定
xrange = {"min-max":(0, 1),
          "Z-score":(-6.5, 6.5),
          "robust Z-score":(-6.5, 6.5)}

# ヒストグラムの作画
for i, (data_name, score_name) in enumerate(concat_df.columns):
    row, col = divmod(i, 3)
    axs[row, col].hist(concat_df[(data_name, score_name)])
    axs[row, col].set_xlim(xrange[score_name])

    title = "\n".join([data_name, score_name])    
    axs[row, col].set_title(title)       

plt.show()

結果が以下の図です。上が外れ値なし、下が外れ値ありの場合です。
min-max normalizationは、外れ値に非常に影響を受けています。Z-scoreも外れ値の影響を受け、外れ値なしの場合と大きく異なっています。robust Z-scoreが最も外れ値の影響が少なく、外れ値なしの場合と比較的似ています。

最後に

robust Z-scoreは、データが正規分布のときはZ-scoreと同じ結果となるので、迷ったらrobust Z-scoreを使おうと考えています。特に、外れ値も使いたい場合にはrobust Z-scoreが有効であると感じました。