【LightGBM】Tweedie回帰でゼロ過剰データを扱う

はじめに

保険金の請求額、店舗ごとの需要、広告のコンバージョン金額——こうしたデータは「多くがゼロ、ときどき大きな正の値」という形をしています。ゼロに大きな山があり、正の側は右に長い裾を引きます。

この形は正規分布を仮定する通常の二乗誤差回帰(objective='regression')とは相性が良くありません。LightGBMには、まさにこの分布のためのobjective='tweedie'が用意されています。この記事で使い方と、通常回帰と比べたときの挙動を実測します。

Tweedie分布

Tweedie分布はtweedie_variance_power(以下 p)というパラメータを1つ持ち、1 < p < 2の範囲ではゼロに確率の塊があり、正の側は連続という分布になります。境界のp=1はポアソン(カウント)、p=2はガンマ(正の連続)に対応し、その間を連続的につなぎます。

1 < p < 2のTweedieは、確率的には「複合ポアソン–ガンマ」——件数がポアソン、1件あたりの大きさがガンマ——として解釈できます。件数が0ならyは0、1件以上なら正の連続値になります。

データを作る

この解釈をそのままコードにして、保険金のようなデータを作ります。事故件数Nをポアソン、1件あたりの金額をガンマ(少額が多く、まれに高額)とし、その合計をyとします。

import numpy as np
import pandas as pd
import lightgbm as lgb

def make(n=60000, seed=0):
    rng = np.random.RandomState(seed)
    x1, x2, x3 = rng.uniform(0, 1, n), rng.uniform(0, 1, n), rng.uniform(0, 1, n)
    lam = np.exp(-1.5 + 1.2*x1 + 0.8*x2 - 0.5*x3)   # 事故頻度(対数線形)
    N = rng.poisson(lam)                            # 事故件数
    # 1件あたり Gamma(0.5, 4.0): 平均2.0 だが裾が重い(まれに高額)
    y = np.array([rng.gamma(0.5, 4.0, k).sum() if k > 0 else 0.0 for k in N])
    return pd.DataFrame({"x1": x1, "x2": x2, "x3": x3}), y

X, y = make()
print(f"ゼロの割合: {(y == 0).mean():.1%}")
print(f"平均: {y.mean():.3f}  最大: {y.max():.2f}")
ゼロの割合: 61.0%
平均: 1.037  最大: 40.24

6割がゼロで、残りが正の裾を引きます。平均は約1ですが、まれに40近い値が出ます。分布を描くと次のようになります。ゼロに大きな山があり、正の側が右に長く伸びています。

lightgbm-tweedie-dist

Tweedieで学習する

objective='tweedie'を指定し、tweedie_variance_powerで p を渡します。今回のデータは複合ポアソン–ガンマなので、まずは中間の1.5にします。

base = dict(num_leaves=31, learning_rate=0.1, min_child_samples=50, verbose=-1)

n = len(X); tr = np.arange(n) < 45000
Xtr, Xte, ytr, yte = X[tr], X[~tr], y[tr], y[~tr]

model = lgb.train({**base, "objective": "tweedie", "tweedie_variance_power": 1.5},
                  lgb.Dataset(Xtr, ytr), num_boost_round=300)
pred = model.predict(Xte)
print(f"予測の最小値: {pred.min():.3f}")   # 負にならない
予測の最小値: 0.009

Tweedie目的関数は内部で対数リンクを使うため、予測は常に非負になります。金額や件数のように負を取らない量では、この性質がそのまま効いてきます。

L2回帰との比較

同じデータを通常の二乗誤差回帰(objective='regression')でも学習し、比べます。評価にはRMSE・MAEに加えて、順序性(大きい値をどれだけ上位に並べられるか)を測る正規化Gini、Tweedie分布のあてはまりを測るTweedieデビアンスを使います。

from sklearn.metrics import mean_tweedie_deviance

rmse = lambda p, t: float(np.sqrt(np.mean((p - t) ** 2)))
mae  = lambda p, t: float(np.mean(np.abs(p - t)))
def gini(y, pred):                       # 正規化Gini(大きいほど順序性が良い)
    o = np.argsort(pred); cum = np.cumsum(y[o]) / y.sum(); g = cum.sum()/len(y) - 0.5
    o2 = np.argsort(y);   cum2 = np.cumsum(y[o2]) / y.sum(); gmax = cum2.sum()/len(y) - 0.5
    return g / gmax

m_tw = lgb.train({**base, "objective": "tweedie", "tweedie_variance_power": 1.5},
                 lgb.Dataset(Xtr, ytr), num_boost_round=300)
m_l2 = lgb.train({**base, "objective": "regression"},
                 lgb.Dataset(Xtr, ytr), num_boost_round=300)
p_tw, p_l2 = m_tw.predict(Xte), m_l2.predict(Xte)

for name, p in [("Tweedie", p_tw), ("L2", p_l2)]:
    dev = mean_tweedie_deviance(yte, np.clip(p, 1e-8, None), power=1.5)
    print(f"{name:8s} RMSE={rmse(p,yte):.4f}  MAE={mae(p,yte):.4f}  "
          f"Gini={gini(yte,p):.4f}  deviance={dev:.4f}  neg={(p<0).mean():.2%}")

5シードの平均は次の通りです。

Tweedie  RMSE=2.5071  MAE=1.3895  Gini=0.2412  deviance= 3.9751  neg=0.00%
L2       RMSE=2.5167  MAE=1.4517  Gini=0.2259  deviance=58.5443  neg=0.22%

RMSEはほとんど変わりません。二乗誤差は裾の大きな値に支配されるため、そこでは差が出にくいためです。一方、MAE・Gini・デビアンスではTweedieが一貫して優り、L2はごく一部で負の予測(金額としてはあり得ない値)を出しています。この差は5シードすべてで同じ向きでした。

図で見ると違いがはっきりします。

lightgbm-tweedie-cmp

左は正解 y と各モデルの予測値の分布です(対数密度)。予測は条件付き平均E[y|x]なので、正解のようなゼロの塊や長い裾は持たず、平均付近に集中した滑らかな分布になります。ここで見えるのは、L2の予測が一部で0を下回っている(左の赤い網掛け)のに対し、Tweedieは対数リンクのおかげで常に0以上に収まっている点です。右は予測十分位ごとの較正です。L2は最上位の十分位で予測が実測を大きく上回っています(少数の高額データに引っ張られて過大に予測している)。Tweedieはこの帯でも実測に近く保たれています。

二乗誤差だけを見ると通常回帰と互角に見えますが、裾の重いゼロ過剰データでは、予測の非負性・順序性・分布のあてはまりでTweedieが優位になります。保険料率の算定や需要予測のように、こうした性質が重要になる場面で効いてきます。

tweedie_variance_power の効果

p を1.1〜1.9で振って、同じデータでのTweedieデビアンス(真のp=1.5で評価、3シード平均)を見ます。

power=1.1: deviance=3.9397
power=1.3: deviance=3.9608
power=1.5: deviance=3.9947
power=1.7: deviance=4.0748
power=1.9: deviance=4.2464

この範囲では差は穏やかで、点予測はpにそれほど敏感ではありませんでした(今回はやや低めのpが良い結果)。pはデータの平均–分散関係(分散 ∝ 平均^p)を決めるパラメータで、真の分布が未知の実データでは交差検証で選ぶのが基本です。ゼロが多く裾が重いほど1に近く、ゼロが少なく連続に近いほど2寄りが目安になります。

注意点

  • 目的変数 y は非負である必要があります(負の値があるとエラーになります)。
  • tweedie_variance_power1 < p < 2で指定します(1はポアソン、2はガンマに対応)。カウントデータならobjective='poisson'、正の連続値だけならobjective='gamma'という選択肢もあります。
  • Tweedie目的関数は対数リンクのため予測は常に非負です。裏を返すと、y が本質的に負を取りうるデータには向きません。
  • RMSEだけで評価すると通常回帰との差が見えにくいので、MAE・Gini・Tweedieデビアンス・予測の非負性など、目的に合った指標で確認します。

まとめ

  • objective='tweedie'は、ゼロ過剰かつ正の裾を持つ連続データ(保険金・需要など)のための目的関数
  • 内部で対数リンクを使うため予測は常に非負
  • 通常の二乗誤差回帰とRMSEはほぼ同等でも、裾が重いデータではMAE・順序性(Gini)・分布のあてはまり(デビアンス)で優位
  • tweedie_variance_power(1〜2)は平均–分散関係を決めるパラメータで、交差検証で選ぶ

参考リンク