概要

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)

以上