概要
imbalancedなclassification問題のsample_weight計算方法を調べてて、regressionならどうなるのだろう、と疑問に感じた。classificationの場合、classの出現頻度の逆数をweightとする。つまり稀なクラスほどweightは大きくなる。
これをregressionに当てはめるなら稀な値、つまりanomalyなyほどweightを大きくすればいいんじゃないかと。そうすると特異検知の話になるのでsample_weightをマハラノビス距離にしてみたらどうだろう、というのが今回の実験備忘録。
classification用のsample_weight計算方法は良い記事を見つけたので、それを試した備忘録を後日投稿することとする。
- 実施期間: 2023年2月
- 環境：Colab
- パケージ：scikit-learn
1. Dataset
sklearnのcalifornia_housingを使用する。
まず、最近お世話になることが多いfeature-engineをinstallしておく。
!pip install feature-engine
importとdatasetのloadまで行う。
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn import datasets, linear_model
from sklearn.pipeline import Pipeline
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error, r2_score
from feature_engine.wrappers import SklearnTransformerWrapper
# DataFrameでloadする。
df_org = datasets.fetch_california_housing(return_X_y=False, as_frame=True).frame
Kaggleでもほぼ同じ内容のcalifornia housing問題があって、そのときに検討した処理をする。sample_weightとは関係しないが、この処理も備忘録として残す。
このdatasetの扱いに困るのがLatitudeとLongitudeだと思う。下図のように完全に２つの山があり、これを直感的に36度線で切ると山を分離することもできるがLatitudeの配慮が全く無くないのでイヤ。
fig, axes = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True) # , sharex="all"
df_org['Latitude'].hist(bins=20, ax=axes[0])
axes[0].set_title('Latitude')
df_org['Longitude'].hist(bins=20, ax=axes[1])
axes[1].set_title('Latitude')
plt.show()
やはりML的に分離したいのでGMMで２つに分離する。ヒューリスティックな36度線分離でヒストグラムはきれいになったので、それぞれの山はガウス分布であると仮定した。
便宜的にcordinate_labelをDataFrameに追加し２つを識別する。
from sklearn.mixture import GaussianMixture
# train, validation, testに分ける。
df = df_org.sample(frac=1, random_state=42, ignore_index=True)
train_df = df.iloc[:10000]
val_df = df.iloc[10000:20000]
test_df = df.iloc[20000:]
# train_dfでfitさせ、
gm = GaussianMixture(n_components=2, n_init=10, random_state=42)
gm.fit(train_df[['Longitude', 'Latitude']])
# fit()で収束したか確認する。
# print(f'Convergence: {gm.converged_}')
# train, validation, testのそれぞれで分離する。
y_pred = gm.predict(train_df[['Longitude', 'Latitude']])
y_pred_prob = gm.predict_proba(train_df[['Longitude', 'Latitude']])
train_df = train_df.assign(cordinate_label = y_pred)
y_pred = gm.predict(val_df[['Longitude', 'Latitude']])
val_df = val_df.assign(cordinate_label = y_pred)
y_pred = gm.predict(test_df[['Longitude', 'Latitude']])
test_df = test_df.assign(cordinate_label = y_pred)
train_df.head()
fig, axis = plt.subplots(1, 2, figsize=(16, 8), tight_layout=True)
train_df[train_df['cordinate_label'] == 0].plot.scatter(x='Longitude', y='Latitude', c='r', ax = axis[0])
train_df[train_df['cordinate_label'] == 1].plot.scatter(x='Longitude', y='Latitude', c='g', ax = axis[0])
val_df[val_df['cordinate_label'] == 0].plot.scatter(x='Longitude', y='Latitude', c='r', ax = axis[1])
val_df[val_df['cordinate_label'] == 1].plot.scatter(x='Longitude', y='Latitude', c='g', ax = axis[1])
LAを中心にする赤点軍と、SFを中心とする緑点群に別れた。
直感的には赤点になってもいいような緑点も内陸部にあるけどsample_weightを検討するには逆にいいのかもしれない。まぁ、境界線があることは明らかに見えるし、ここにはあまり手を入れないほうがいい気がする。
Longitude, Latitudeの分布を確認する。
fig, axes = plt.subplots(2, 2, figsize=(10, 8), tight_layout=True) # , sharex="all"
train_df[train_df['cordinate_label']==0]['Latitude'].hist(bins=20, ax=axes[0, 0])
axes[0, 0].set_title('Latitude_r')
train_df[train_df['cordinate_label']==0]['Longitude'].hist(bins=20, ax=axes[0, 1])
axes[0, 1].set_title('Longitude_r')
train_df[train_df['cordinate_label']==1]['Latitude'].hist(bins=20, ax=axes[1, 0])
axes[1, 0].set_title('Latitude_g')
train_df[train_df['cordinate_label']==1]['Longitude'].hist(bins=20, ax=axes[1, 1])
axes[1, 1].set_title('Longitude_g')
座標系を２つに分け、
def divide_into_2regions(df):
# red: LA
df['Latitude_r'] = df[df['cordinate_label']==0]['Latitude']
df['Longitude_r'] = df[df['cordinate_label']==0]['Longitude']
df['cordinate_label_r'] = (df['cordinate_label']==0)
# green: San Francisco
df['Latitude_g'] = df[df['cordinate_label']==1]['Latitude']
df['Longitude_g'] = df[df['cordinate_label']==1]['Longitude']
df['cordinate_label_g'] = (df['cordinate_label']==1)
df = df.drop(['Latitude', 'Longitude', 'cordinate_label'], axis=1)
return df
train_df = divide_into_2regions(train_df)
val_df = divide_into_2regions(val_df)
test_df = divide_into_2regions(test_df)
train_df.head()
さらにredのみ、greenのみのdatasetをこしらえ、とりあえずgreenのdatasetだけを使ってsample_weightの評価を行う。
train_df_r = train_df[train_df['cordinate_label_r']==True]
train_df_r = train_df_r.drop(['Latitude_g', 'Longitude_g', 'cordinate_label_g'], axis=1)
train_df_g = train_df[train_df['cordinate_label_g']==True]
train_df_g = train_df_g.drop(['Latitude_r', 'Longitude_r', 'cordinate_label_r'], axis=1)
val_df_r = val_df[val_df['cordinate_label_r']==True]
val_df_r = val_df_r.drop(['Latitude_g', 'Longitude_g', 'cordinate_label_g'], axis=1)
val_df_g = val_df[val_df['cordinate_label_g']==True]
val_df_g = val_df_g.drop(['Latitude_r', 'Longitude_r', 'cordinate_label_r'], axis=1)
test_df_r = test_df[test_df['cordinate_label_r']==True]
test_df_r = test_df_r.drop(['Latitude_g', 'Longitude_g', 'cordinate_label_g'], axis=1)
test_df_g = test_df[test_df['cordinate_label_g']==True]
test_df_g = test_df_g.drop(['Latitude_r', 'Longitude_r', 'cordinate_label_r'], axis=1)
分布を確認する。
def draw_all_hist(df, region):
fig, axes = plt.subplots(2, 4, figsize=(20, 8), tight_layout=True) # , sharex="all"
df['MedInc'].hist(bins=20, ax=axes[0, 0])
axes[0, 0].set_title('MedInc')
df['HouseAge'].hist(bins=20, ax=axes[0, 1])
axes[0, 1].set_title('HouseAge')
df['AveRooms'].hist(bins=20, ax=axes[0, 2])
axes[0, 2].set_title('AveRooms')
df['AveBedrms'].hist(bins=20, ax=axes[0, 3])
axes[0, 3].set_title('AveBedrms')
df['Population'].hist(bins=20, ax=axes[1, 0])
axes[1, 0].set_title('Population')
df['AveOccup'].hist(bins=20, ax=axes[1, 1])
axes[1, 1].set_title('AveOccup')
if region == 'r':
df['Latitude_r'].hist(bins=20, ax=axes[1, 2])
axes[1, 2].set_title('Latitude_r')
df['Longitude_r'].hist(bins=20, ax=axes[1, 3])
axes[1, 3].set_title('Longitude_r')
else:
df['Latitude_g'].hist(bins=20, ax=axes[1, 2])
axes[1, 2].set_title('Latitude_g')
df['Longitude_g'].hist(bins=20, ax=axes[1, 3])
axes[1, 3].set_title('Longitude_g')
plt.show()
draw_all_hist(train_df_g, 'g')
Outlierを除去する場合
def find_normal_boundaries(df, variable):
# Calculate the boundaries
# for a Gaussian distribution.
upper_boundary = df[variable].mean() + 3 * df[variable].std()
lower_boundary = df[variable].mean() - 3 * df[variable].std()
return upper_boundary, lower_boundary
def find_skewed_boundaries(df, variable, distance):
# Let's calculate the boundaries
# for skewed distributions
# The parameter "distance" gives us the option to
# estimate 1.5 times or 3 times the IQR when defining
# the boundaries.
IQR = df[variable].quantile(0.75) - df[variable].quantile(0.25)
lower_boundary = df[variable].quantile(0.25) - (IQR * distance)
upper_boundary = df[variable].quantile(0.75) + (IQR * distance)
return upper_boundary, lower_boundary
'AveRooms', 'AveBedrms', 'Population', 'AveOccup', 'Longitude_g'は外れ値があるので、IQRの３倍以上外れた点を除去する。
train_df_g_no_out = train_df_r.copy()
high_bnd, low_bnd = find_skewed_boundaries(train_df_g_no_out, 'AveRooms', 3.)
train_df_g_no_out = train_df_g_no_out[(train_df_g_no_out['AveRooms']>low_bnd) & (train_df_g_no_out['AveRooms']<high_bnd)]
high_bnd, low_bnd = find_skewed_boundaries(train_df_g_no_out, 'AveBedrms', 3.)
train_df_g_no_out = train_df_g_no_out[(train_df_g_no_out['AveBedrms']>low_bnd) & (train_df_g_no_out['AveBedrms']<high_bnd)]
high_bnd, low_bnd = find_skewed_boundaries(train_df_g_no_out, 'Population', 3.)
train_df_g_no_out = train_df_g_no_out[(train_df_g_no_out['Population']>low_bnd) & (train_df_g_no_out['Population']<high_bnd)]
high_bnd, low_bnd = find_skewed_boundaries(train_df_g_no_out, 'AveOccup', 3.)
train_df_g_no_out = train_df_g_no_out[(train_df_g_no_out['AveOccup']>low_bnd) & (train_df_g_no_out['AveOccup']<high_bnd)]
high_bnd, low_bnd = find_skewed_boundaries(train_df_g_no_out, 'Longitude_g', 3.)
train_df_g_no_out = train_df_g_no_out[(train_df_g_no_out['Longitude_g']>low_bnd) & (train_df_g_no_out['Longitude_g']<high_bnd)]
train_df_g_no_out.describe()
print(f'{len(train_df_g) - len(train_df_g_no_out)} samples were omitted.')
除去後の分布を確認する。
draw_all_hist(train_df_g_no_out, 'g')
不要な'cordinate_label_g'列を落とし、train, val, test datasetについてそれぞれX, yに分ける。
train_df_g_no_out = train_df_g_no_out.drop(['cordinate_label_g'], axis=1)
train_X = train_df_g_no_out.drop(['MedHouseVal'], axis=1)
train_y = train_df_g_no_out['MedHouseVal']
val_df_g_no_out = val_df_g.copy().drop(['cordinate_label_g'], axis=1)
val_X = val_df_g_no_out.drop(['MedHouseVal'], axis=1)
val_y = val_df_g_no_out['MedHouseVal']
test_df_g_no_out = test_df_g.copy().drop(['cordinate_label_g'], axis=1)
test_X = test_df_g_no_out.drop(['MedHouseVal'], axis=1)
test_y = test_df_g_no_out['MedHouseVal']
不要な'cordinate_label_g'列を落とし、train, val, test datasetについてそれぞれX, yに分ける。
# greenの場合は下記を実行
train_df_g_no_out = train_df_g.copy().drop(['cordinate_label_g'], axis=1)
train_X = train_df_g_no_out.drop(['MedHouseVal'], axis=1)
train_y = train_df_g_no_out['MedHouseVal']
val_df_g_no_out = val_df_g.copy().drop(['cordinate_label_g'], axis=1)
val_X = val_df_g_no_out.drop(['MedHouseVal'], axis=1)
val_y = val_df_g_no_out['MedHouseVal']
test_df_g_no_out = test_df_g.copy().drop(['cordinate_label_g'], axis=1)
test_X = test_df_g_no_out.drop(['MedHouseVal'], axis=1)
test_y = test_df_g_no_out['MedHouseVal']
y('MedHouseVal')を確認する。
train_y.sort_values(ignore_index=True).plot()
2. Scale用のPipe line
Feature-engineを積極的に使用する。
from feature_engine.transformation import LogTransformer, BoxCoxTransformer
from feature_engine.outliers import Winsorizer
from sklearn.preprocessing import StandardScaler, Normalizer
# from sklearn.preprocessing import QuantileTransformer
# quantile_transformer = QuantileTransformer(output_distribution='normal', random_state=0)
各featureの分布から、pipe lineは下記とする。
Feature-engineで用意されていない変換はsklearnを使用するが、これもSklearnTransformerWrapper()でラップすることによりDataFrameをシームレスに使用できる。
all_clmns = list(train_X.columns)
preprocessor = Pipeline([
('iqr_cap_right', Winsorizer(capping_method='iqr', tail='right', fold=3, variables=['MedInc', 'AveRooms', 'AveBedrms', 'Population', 'AveOccup'])),
('iqr_cap_both', Winsorizer(capping_method='iqr', tail='both', fold=3, variables=['Longitude_g'])), # , add_indicators=True
('boxcox', BoxCoxTransformer(variables=['MedInc', 'Population', 'AveOccup'])),
('std_scaler', SklearnTransformerWrapper(transformer = StandardScaler(), variables = all_clmns)),
# ('norm_scaler', SklearnTransformerWrapper(transformer = Normalizer(), variables = all_clmns)),
])
preprocessor.fit(train_X)
Pipeline(steps=[('iqr_cap_right',
Winsorizer(capping_method='iqr',
variables=['MedInc', 'AveRooms', 'AveBedrms',
'Population', 'AveOccup'])),
('iqr_cap_both',
Winsorizer(capping_method='iqr', tail='both',
variables=['Longitude_g'])),
('boxcox',
BoxCoxTransformer(variables=['MedInc', 'Population',
'AveOccup'])),
('std_scaler',
SklearnTransformerWrapper(transformer=StandardScaler(),
variables=['MedInc', 'HouseAge',
'AveRooms', 'AveBedrms',
'Population', 'AveOccup',
'Latitude_g',
'Longitude_g']))])
train_xf_X = preprocessor.transform(train_X)
val_xf_X = preprocessor.transform(val_X)
test_xf_X = preprocessor.transform(test_X)
draw_all_hist(train_xf_X, 'g')
3. sample_weight計算
transform後のtrain_xf_Xについてマハラノビス距離を求め、そのn乗をsample_weightとする。nを振って精度の変化を確認する作戦。
まず全体について確認してみる。
from sklearn.covariance import MinCovDet
robust_cov = MinCovDet().fit(train_xf_X)
print('Estimated covariance matrix:\n', robust_cov.covariance_)
# mahal_robust_cov = robust_cov.mahalanobis(train_cov_X)
mahal_robust_cov = robust_cov.mahalanobis(train_xf_X)
print(mahal_robust_cov.shape)
print(mahal_robust_cov.max(), mahal_robust_cov.min())
counts, bins = np.histogram(mahal_robust_cov)
plt.hist(bins[:-1], bins, weights=counts)
plt.show()
上図に現れる少量のマハラノビス距離が離れたsampleを特異としてweightできるか？
4. RandomForestRegressor modelで評価
ここでは無調整のRandomForestRegressorを使用する。sample数とn乗を振りmseで精度比較する。各条件で得る精度はdfに入れる。mseの平均で評価するため各条件にて10回fitさせている。
尚、pow=.0 の時はsample_weightは全部1.になるのでsample_weihgt指定しない時と同じ結果になる（確認済み）。
def predict_with_sample_weight(train_xf_X, train_y, val_xf_X, val_y, pow=0.):
robust_cov = MinCovDet().fit(train_xf_X)
mahal_robust_cov = robust_cov.mahalanobis(train_xf_X)
sample_weight = mahal_robust_cov
sample_weight = np.power(sample_weight, pow)
fit_params={'sample_weight': sample_weight}
mse_val_lst = []
r2_val_lst = []
np.random.seed(seed=23)
for i in range(10):
reg = RandomForestRegressor()
reg.fit(train_xf_X, train_y, sample_weight)
# reg.fit(train_xf_X, train_y, sample_weight)
y_pred = reg.predict(val_xf_X)
mse_val_lst.append(mean_squared_error(val_y, y_pred))
r2_val_lst.append(r2_score(val_y, y_pred))
return mse_val_lst, r2_val_lst
data_num_lst = [500, 1000, 1500, 2000, 3000]
pow_lst = [0., 0.1, 0.2, 0.3, 0.4, 0.5, 1., 1.5, 2.]
df = pd.DataFrame()
for data_num in data_num_lst:
for pow in pow_lst:
mse_val_lst, r2_val_lst = predict_with_sample_weight(train_xf_X[:data_num], train_y[:data_num], val_xf_X, val_y, pow)
clmn_name = str(data_num) + '_' + str(pow)
s = pd.Series(mse_val_lst, name=clmn_name)
df = pd.concat([df, s], axis=1)
mse結果表示用に泥臭いdf処理
df_1 = df.copy() # df計算に時間がかかったので、念の為copyする
s_mean = df_1.mean(axis = 0)
s_std = df_1.std(axis = 0)
df_mean = pd.DataFrame()
df_std = pd.DataFrame()
for i, data_num in enumerate(data_num_lst):
s_wk = s_mean[i*9:(i+1)*9].reset_index(drop=True)
s_wk = s_wk.rename(str(data_num))
df_mean = pd.concat([df_mean, s_wk], axis=1)
df_mean['pow'] = pow_lst
df_mean = df_mean.set_index('pow')
fig, ax = plt.subplots(figsize=(8, 6), tight_layout=True)
ax.set_ylabel('mse')
df_mean.plot(ax=ax)
ax.grid(axis = 'y', c = 'black', alpha = 0.3, linestyle = '--', linewidth = 1)
plt.show()
先にも書いた通り、pow=.0はsample_weightが無いときと等しくなる。
どのsample数でも、距離の0.4乗前後でsample_weight無しよりも若干mseが小さくなっていることがわかる。もちろんデータによるので今回はこういう結果であって、違うデータならどうなるかやってみなければわからない。
マハラノビス距離以外でも特異性を各sampleで予め計算できるのであれば、それをsample_weightにすることは可能だろう。
次回は違う非線形modelをoutlierを取り除いたデータでTrainingし（つまりoutlierの存在を知らないmodelを作成）、それにoutlierを含む別データでPredictした結果のmseかmaeを距離に見立ててsample_weightとするケースを実験してみたい。
あと、classificationも。
余談
当初10回fitさせる部分はcross_validate(cv=10)で計算しようとしたんだけれど、fit_paramsに入れたsample_weightの長さがXの長さと一致しないエラーが発生して断念した。Xの数が10で割り切れないと最後の１回分は他の9回と比べて長さが異なって、それにsample_weightが実装上対応できていないような気がした。
時間があるときに原因深堀りしてbugとして報告すること。
sample_weight = np.power(sample_weight, n)
fit_params={'sample_weight': sample_weight}
scores = cross_validate(reg, val_xf_X, val_y, cv=10,
scoring=['neg_mean_squared_error', 'r2'],
fit_params=fit_params)
以上