Box-Cox Transormationは次の式による変換。λ=0のときはlog(x)。
λが1より大きい場合は小さな値の間隔が圧縮され、小さい場合は大きな値の間隔が圧縮されるように変換される。
import numpy as np
from scipy.special import boxcox1p
import matplotlib.pyplot as plt
from bokeh.plotting import figure
from bokeh.io import output_notebook, show
output_notebook()
p = figure(
title="Box-Cox Transformations",
x_axis_label='x',
y_axis_label='λ',
)
for lam in [-1, 0, 1, 2]:
v = np.array([boxcox1p(i, lam) for i in range(10)])
v = v / v.max()
p.line(v, lam)
p.circle(v, lam, size=8)
show(p)
これによって左右非対称な分布を対称(skew=0)な正規分布に近づけることができる。 以前正規分布に近づけるのに対数を取ったが、これはBox-Cox transformationの1ケースだといえる。
KaggleのHome Prices CompetitionのKernelからデータの探り方を学ぶ - sambaiz-net
試しに適当な左右非対称な分布のデータを変換してみる。
import pandas as pd
import seaborn as sns
p = np.array([500, 2000, 3000, 2500, 1000, 500, 125, 125, 0.0125, 0.0125])
p = p / p.sum()
data = np.random.choice([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], size=1000, p=p) * np.random.rand(1000)
df = pd.DataFrame({
'x': data
})
print('skew={:.3}'.format(df['x'].skew())) # skew=1.08
sns.distplot(df['x'], hist=False, label=('data'))
sns.distplot(np.random.normal(df['x'].mean(), np.sqrt(df['x'].var()), len(data)), hist=False, label=('normal distribution'))
評価値としてソートしたデータと正規分布のCDF(cumulative distribution function; 累積分布関数)の逆関数の相関係数を用いる。 逆関数はppf()で、(0,1)の定義域で(-∞,∞)の値を取る。
from scipy.stats import norm
import matplotlib.pyplot as plt
x = [x * 0.001 for x in range(0, 1001)]
y = norm.ppf(x) # inverse CDF
plt.plot(x, y)
-2~2まで0.1刻みでλを動かして相関係数をプロットする。
from scipy.special import boxcox1p
lams = [ i * 0.1 for i in range(-20, 20) ]
data.sort()
xs = [norm.ppf(((i+1) - 0.5) / len(data)) for i in range(len(data))]
yss = [np.array([boxcox1p(d, lam) for d in data]) for lam in lams] # boxcox1p(x) = ((x + 1) ** lam - 1) / lam
cors = [np.corrcoef(xs, ys)[0, 1] for ys in yss]
best_cor = -1.0
for i, cor in enumerate(cors):
if best_cor < cor:
best_cor = cor
best_lam = lams[i]
best_ys = yss[i]
plt.plot(lams, cors)
最良のλで変換したところskewが大きく減少し正規分布に近づいた。
from scipy.stats import skew
print('lambda={} skew={:.3}'.format(best_lam, skew(best_ys))) # lambda=-0.1 skew=0.11
sns.distplot(best_ys, hist=False, label=('transformed data'))
sns.distplot(np.random.normal(best_ys.mean(), np.sqrt(best_ys.var()), len(best_ys)), hist=False, label=('normal distribution'))
inv_boxcox1pで元に戻せる。
from scipy.special import inv_boxcox1p
inv_ys = inv_boxcox1p(best_ys, best_lam)