はじめに

回帰分析でRMSPEやMAPEといった目的変数の正解と残差の比を評価する場合、目的変数をlog変換すると精度が上がる場合があります。これについて考察して見たいと思います。

環境

macOS Majave 10.14.6
Python 3.6.8 :: Anaconda, Inc.
lightgbm 2.2.3
sklearn 0.19.1

手順

観測値を$y_i$、予測値を$\hat y_i$とします。添字$i=1,2,...,N$はデータの添字です。RMSPEやMAPEでは残差$y_i-\hat y_i=\varepsilon_i$を観測値$y_i$で割ります。良く訓練された予測器では、残差は観測値と相関がなく、残差の分散は観測値によらず一定です。したがってRMSPEやMAPEの和のうち、観測値が小さい項は値が大きくなり、評価に大きな影響を及ぼします。

\begin{align}
\mathrm{RMSPE} &= \sqrt{\frac{1}{N}\sum_{i=1}^N\Big(\frac{y_i - \hat y_i}{y_i}\Big)^2} \\
\mathrm{MAPE} &= \frac{1}{N}\sum_{i=1}^N\Big|\frac{y_i - \hat y_i}{y_i}\Big|
\end{align}

これをモデリングで解決しようとすると、観測値$y_i$が大きいときは精度が悪く、観測値$y_i$が小さいときは精度が良い、というモデルを作らなければなりません。
ポアソン回帰などはそういった効果を取り込むことができますが、今回は目的変数の変換だけでどれだけの効果があるか調べてみたいと思います。

MAPE

目的変数が0より大きい値を取るデータを仮定します。対数変換したあとの目的変数を$\tilde y_i = \log (y_i/C)$とします。$C$は$y_i$の次元と同じ次元を持つ定数です。回帰分析として$\tilde y_i$を予測する予測器を訓練します。この予測器で予測した値を$\hat{\tilde y_i}$とします。
残差を$\tilde{\epsilon_i} = \tilde{y_i} - \hat{\tilde y_i}$とします。予測器は十分な精度をもっている、すなわち$\epsilon_i \sim \mathcal N(0, \tilde \sigma^2)$を満たすとします。予測値を元のスケールに戻すには指数変換します。

\hat{y_i} = C e^{\hat{\tilde y_i}}

元の$y_i$でのMAPEを計算します。

\begin{align}
\mathrm{MAPE} & = \frac{1}{N}\sum_{i=1}^N\Big| \frac{y_i - \hat y_i}{y_i} \Big| \\
& = \frac{1}{N}\sum_{i=1}^N\Big| \frac{e^{\tilde y_i} - e^{\hat {\tilde y_i}}}{e^{\tilde y_i}} \Big| \\
& = \frac{1}{N}\sum_{i=1}^N\Big| 1 - e^{\hat {\tilde y_i} - \tilde y_i} \Big| \\
& = \frac{1}{N}\sum_{i=1}^N\Big| 1 - e^{-\tilde\epsilon_i} \Big|
\end{align}

MAPEの式が$\tilde\epsilon_i$のみで書かれ、$y_i$がなくなっています。したがってMAPEをスケールによらない無次元な量で表す事ができました。

$\tilde\epsilon_i$を確率変数と考え、このMAPEの期待値を計算します。計算は付録に書きました。

E[\mathrm{MAPE}] = \sqrt{\frac{2}{\pi}} e^{\frac{\tilde\sigma^2}{2}}\int_0^\tilde\sigma dy \exp \Big[-\frac{y^2}{2}\Big]

もう少し意味をわかりやすくするため、近似してみましょう。$\tilde\sigma \ll 1$を仮定します。

E[\mathrm{MAPE}] \simeq \sqrt{\frac{2}{\pi}}\tilde\sigma + \mathcal O(\tilde\sigma^3)

よってMAPEは対数変換したあとの予測における残差の分散程度になることが分かります。

MAPEの二乗の期待値は次のようになります。

\begin{align}
E[(\mathrm{MAPE})^2] &\simeq \frac{1}{N}(1-2e^{\frac{1}{2}\tilde\sigma^2}+e^{2\tilde\sigma^2}) + \sqrt{\frac{2}{\pi}}\tilde\sigma^2 \\
&\simeq \frac{1}{N}\tilde\sigma^2 + \sqrt{\frac{2}{\pi}}\tilde\sigma^2
\end{align}

したがって分散は$\tilde\sigma^2/N$となります。

RMSPE

\begin{align}
\mathrm{RMSPE} & = \sqrt{\frac{1}{N}\sum_{i=1}^N\Big(\frac{y_i - \hat y_i}{y_i}\Big)^2} \\
& = \sqrt{\frac{1}{N}\sum_{i=1}^N(1-2e^{-\tilde\epsilon_i} + e^{-2\tilde\epsilon_i})}
\end{align}

これ以上は厳密に計算できないので、荒っぽいですが$\epsilon_i \ll 1$として一番寄与する項のみ取り出して評価してみます。

1 - 2 e^{-x} + e^{-2 x} \simeq x^2 + \mathcal O(x^3)

テイラー展開の結果から次のようになります。

\begin{align}
\mathrm{RMSPE} & = \sqrt{\frac{1}{N}\sum_{i=1}^N \epsilon_i^2}
\end{align}

これの期待値を取ります。球座標$r^2=\sum_i \epsilon_i^2$にしてやれば計算できます。

\begin{align}
\mathrm{RMSPE} & = \sqrt{\frac{1}{N}}\int_{-\infty}^{\infty}\prod_i(d\epsilon_i 
 \frac{1}{\sqrt{2\pi\tilde \sigma^2}}e^{-\frac{\epsilon_i^2}{2\tilde\sigma^2}})\sqrt{\sum_{i=1}^N  \epsilon_i^2}\\
&= \sqrt{\frac{1}{N}}\frac{1}{(2\pi\tilde \sigma^2)^{\frac{N}{2}}}\frac{2\pi^{\frac{N}{2}}}{\Gamma(\frac{N}{2})}2^{\frac{N-1}{2}}(\frac{1}{\tilde\sigma})^{-\frac{N+1}{2}}\Gamma(\frac{N+1}{2})\\
&= \sqrt{\frac{2}{N}}\frac{\Gamma(\frac{N+1}{2})}{\Gamma(\frac{N}{2})}\tilde\sigma
\end{align}

ここで$N$が十分大きいとして、ガンマ関数の比を$\infty$周りで級数展開します。

\frac{\Gamma(\frac{N+1}{2})}{\Gamma(\frac{N}{2})} = \sqrt{\frac{N}{2}}+\mathcal O(\sqrt{\frac{1}{N}})

以上をまとめると、$\tilde\sigma^2\ll 1$かつ$N\sim\infty$であればRMSPEは次のようになります。

\mathrm{RMSPE} \simeq \tilde\sigma

RMSPEも対数変換したあとの予測における残差の分散程度になることが分かりました。

ちなみにRMSPEの二乗は厳密に計算できます。

\begin{align}
E[(\mathrm{RMSPE})^2] &= E[\frac{1}{N}\sum_{i=1}^N(1-2e^{-\tilde\epsilon_i} + e^{-2\tilde\epsilon_i} )] \\
&= \frac{1}{N}\sum_{i=1}^N(1-2e^{\frac{1}{2}\tilde\sigma^2}+e^{2\tilde\sigma^2})\\
&\simeq\tilde\sigma^2
\end{align}

最後のイコールは$\tilde\sigma \ll 1$を仮定して近似しました。RMSPEの分散を計算するには次のオーダーが必要ですが、およそ$\tilde\sigma^4$のオーダーなので、平均値の10%前後となると考えられます。

実験

実験、といっても値が同じくらいになるか確認だけしましょう。本当はデータをシャッフルしたりモデルを変えたりしたいですが。

import os
import random as rn
import numpy as np
import lightgbm as lgb
import matplotlib.pyplot as plt
from sklearn.datasets import load_boston
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import KFold, train_test_split


os.environ['PYTHONHASHSEED'] = '0'
np.random.seed(0)
rn.seed(0)


def train_predict(X_train, y_train, X_test, y_test):
    # create dataset for lightgbm
    lgb_train = lgb.Dataset(X_train, y_train)
    lgb_eval = lgb.Dataset(X_test, y_test, reference=lgb_train)
    # LightGBM parameters
    params = {
            'task' : 'train',
            'boosting_type' : 'gbdt',
            'objective' : 'regression',
            'metric' : {'l2'},
            'num_leaves' : 31,
            'learning_rate' : 0.1,
            'feature_fraction' : 0.9,
            'bagging_fraction' : 0.8,
            'bagging_freq': 5,
            'verbose' : 0
    }

    # train
    gbm = lgb.train(params,
                lgb_train,
                num_boost_round=100,
                valid_sets=lgb_eval,
                early_stopping_rounds=10)
    y_pred = gbm.predict(X_test, num_iteration=gbm.best_iteration)
    return y_pred


def RMSPE(y_test, y_pred):
    return np.sqrt(np.mean((y_test - y_pred)**2/y_test**2))


def MAPE(y_test, y_pred):
    return np.mean(np.abs((y_test - y_pred)/y_test))


boston = load_boston()
y = boston['target']
X = boston['data']

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=999
)
y_train_log, y_test_log = np.log(y_train), np.log(y_test)

# log変換なし
y_pred = train_predict(X_train, y_train, X_test, y_test)

# log変換あり
y_pred_log = train_predict(X_train, y_train_log, X_test, y_test_log)
y_pred_log_exp = np.exp(y_pred_log)


print('RMSPE:', RMSPE(y_test, y_pred))
print('RMSPE_log:', RMSPE(y_test, y_pred_log_exp))

print('MAPE:', MAPE(y_test, y_pred))
print('MAPE_log:', MAPE(y_test, y_pred_log_exp))

sigma_log = np.sqrt(np.var(epsilon_log))
print('sigma_log:', sigma_log)

print('sqrt(2/pi)*MAPE', np.sqrt(2/np.pi)*sigma_log)

結果を見ると、$\tilde \sigma$とRMSPEが、$\sqrt(2/\pi)\tilde\sigma$とMAPEがほとんど一致しています。また、今回のデータではlog変換したほうが精度が上がっています。厳密には、これはただの期待値なので、分散を調べてやる必要がありますが、今回はここまでで良いでしょう。

内容	値
RMSPE	0.1452407772644214
RMSPE_log	0.12980324752007738
MAPE	0.10628665495403554
MAPE_log	0.09833765994684592
sigma_log	0.12943135929115654
sqrt(2/pi)*sigma_log	0.10327128326214231

おわりに

目的変数の変換は、一般的にはバイアスを生じさせたり誤差を増幅させるのでしない方がよいですが、評価関数によっては適切な変換が効果をもたらすようです。

こちらも参考になります。

付録

$x\sim \mathcal N(0, \sigma)$のとき、$|1-e^{-x}|$の期待値を計算します。

\begin{align}
E_x[|1-e^{-x}|] &= \frac{1}{\sqrt{2\pi\sigma^2}}\int_{-\infty}^{\infty}dx |1-e^{-x}| \exp[-\frac{x^2}{2\sigma^2}]\\
&= \frac{1}{\sqrt{2\pi\sigma^2}}\int_{0}^{\infty}dx (1-e^{-x}) \exp[-\frac{x^2}{2\sigma^2}] - \frac{1}{\sqrt{2\pi\sigma^2}}\int_{-\infty}^{0}dx (1-e^{-x}) \exp[-\frac{x^2}{2\sigma^2}]\\
&= \frac{1}{\sqrt{2\pi\sigma^2}}\int_{0}^{\infty}dx  (e^x-e^{-x})\exp[-\frac{x^2}{2\sigma^2}]\\
&= \frac{1}{\sqrt{2\pi\sigma^2}}\int_{0}^{\infty}dx \Big( \exp \Big[-\frac{(x-\sigma^2)^2}{2\sigma^2}+\frac{\sigma^2}{2}\Big] - \exp \Big[-\frac{(x+\sigma^2)^2}{2\sigma^2}+\frac{\sigma^2}{2}\Big] \Big) \\
&= \frac{1}{\sqrt{2\pi\sigma^2}}\int_{\sigma^2}^{\infty}dx  \exp \Big[-\frac{(x-\sigma^2)^2}{2\sigma^2}+\frac{\sigma^2}{2}\Big] 
+ \frac{1}{\sqrt{2\pi\sigma^2}}\int_{0}^{\sigma^2}dx  \exp \Big[-\frac{(x-\sigma^2)^2}{2\sigma^2}+\frac{\sigma^2}{2}\Big]\\
& - \frac{1}{\sqrt{2\pi\sigma^2}}\int_{-\sigma^2}^{\infty}dx\exp \Big[-\frac{(x+\sigma^2)^2}{2\sigma^2}+\frac{\sigma^2}{2}\Big]
- \frac{1}{\sqrt{2\pi\sigma^2}}\int_{-\sigma^2}^{0}dx\exp \Big[-\frac{(x+\sigma^2)^2}{2\sigma^2}+\frac{\sigma^2}{2}\Big]\\
&= \frac{1}{\sigma}\sqrt{\frac{2}{\pi}}\int_{0}^{\sigma^2}dx  \exp \Big[-\frac{(x-\sigma^2)^2}{2\sigma^2}+\frac{\sigma^2}{2}\Big] \\
&=\sqrt{\frac{2}{\pi}} e^{\frac{\sigma^2}{2}}\int_0^\sigma dy \exp \Big[-\frac{y^2}{2}\Big] 
\end{align}

最後のイコールは$y = x/\sigma - \sigma$と変数変換しました。