0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@rino_yume(rino yume)

Zindi「A Step Ahead of Drought」の実装を組み直した話 ― TWS予測でベースラインから一歩抜けるために考えたこと

0
Posted at

Zindi「A Step Ahead of Drought」の実装を組み直した話 ― TWS予測でベースラインから一歩抜けるために考えたこと

はじめに

前回の記事で、GRACE/GRACE-FO由来のTWS（陸域総水貯留量）を1か月先まで予測するZindiのコンペを紹介して、Persistence ModelとLightGBMのベースラインまで作りました。

ただ、ベースラインを動かしただけだと当然そこまで順位は伸びません。今回はその先、「どこをどう変えたら上に行けそうか」を考えながら実装を組み直したので、その過程をまとめておきます。完成形は python pipeline/02_train_and_submit.py を一回叩けば、生データから提出ファイルまで通るようにしてあります。

評価指標はRMSEですが、このコンペは最終評価でAI TrustworthinessとInnovation & Practicalityも見られます。ここが地味に効いてくるので、後半でその話もします。

まず、Persistenceの壁をどう越えるか

最初にぶつかるのがこれでした。TWSは前月の値との相関がとにかく強くて、TWS_{t+1} = TWS_t とそのまま置くだけのPersistence Modelが、思った以上に強いんです。下手に凝ったモデルを作っても、これを下回れないことが普通にあります。

そこで発想を変えて、TWS_{t+1}を直接当てるのをやめました。代わりに「前月からどれだけ変化したか」を学習させて、最後にTWS_tを足して戻す。要するに差分を予測する形です。

# delta   : target = TWS_{t+1} - TWS_t        → 予測に TWS_t を足し戻す
# anomaly : target = TWS_{t+1} - clim(cell,月)  → 気候値を足し戻す
def define_target(df, cfg):
    d = cfg.data
    if cfg.target.mode == "delta":
        y = df[d.target_col] - df[d.tws_col]
        base = df[d.tws_col].copy()
    elif cfg.target.mode == "anomaly":
        y = df[d.target_col] - df["tws_clim_mean"]
        base = df["tws_clim_mean"].copy()
    else:  # raw（直接予測）
        y, base = df[d.target_col].copy(), 0.0
    return y, base

水準そのものを毎回学習し直すより、変化量という一番難しい部分にモデルのリソースを集中させられるので、RMSEが素直に下がってくれました。raw（直接予測）も残してあるので、3通り試して比べられます。

季節と地域の構造を特徴量に押し込む

EDAで月別・緯度帯別のTWSを眺めていると、結局のところTWSの変動は「その場所の季節サイクル」でだいたい説明できることが見えてきます。だったらその季節サイクルそのものを特徴量にしてしまえばいい、ということで気候値（climatology）を作りました。

セルごと・月ごとのTWS平均と標準偏差を出して、現在値がそこからどれだけズレているか（偏差・標準化偏差）を入れます。

def fit_climatology(train, cfg):
    d = cfg.data
    cell = [d.lat_col, d.lon_col]
    g = (train.groupby(cell + ["month"])[d.tws_col]
              .agg(["mean", "std", "min", "max"]).reset_index())
    return g  # セル×月の気候値テーブル

# 適用側
df["tws_anom_t"] = df[d.tws_col] - df["tws_clim_mean"]            # 偏差
df["tws_z_t"]    = df["tws_anom_t"] / (df["tws_clim_std"] + 1e-6)  # 標準化偏差

注意点としては、この気候値は必ず学習データだけから計算すること。test側の情報を混ぜると検証スコアだけ良くなって、本番でずっこけます。地味ですがここでリークさせると全部台無しなので、build_features.pyの中でtrainだけからfitしてからtest含む全体にapplyする順番を固定しています。

このほか、過去ラグ・移動平均・トレンド（temporal）、緯度帯やKMeansによる空間クラスタ・近傍セルの平均（spatial）、SPEIの短期と長期の差（hydro）あたりも入れていますが、効き方の感触として一番大きかったのは差分ターゲットと気候値の組み合わせでした。

一番神経を使ったのは検証の作り方

正直、モデルやチューニングよりここが本番だと思っています。提出は1日10回・通算300回までで、最終順位はPrivate Leaderboardで決まる。testは未来の月なので時系列の構造があります。

この状況で普通のRandom KFoldをやると、学習と検証に時間的・空間的に近いデータが混ざって、実力より良いスコアが出ます。これを信じて突っ走ると最終日に泣くやつです。

なので、本番の状況に寄せて「過去で学習して、その直後の月を検証する」前向き分割を基本にしました。

def time_forward_splits(df, cfg):
    """過去で学習 → 直後の月を検証（本番に一番近い）"""
    d = cfg.data
    months = np.sort(df[d.date_col].unique())
    for vm in months[-cfg.validation.n_splits:]:
        tr_idx = np.where(df[d.date_col].values < vm)[0]
        va_idx = np.where(df[d.date_col].values == vm)[0]
        yield tr_idx, va_idx

加えて、未知の地域にどれだけ通用するかを見たいので、空間ブロックで切る分割も用意しました。用途で使い分ける感じです。

分割方法 何を見たいか
time（前向き） 将来の月への予測力。基本これ
spatial（空間ブロック） 知らない地域への汎化
group_month 同じ月が学習と検証に跨らないように
random 最初の感触をつかむ用。過信は禁物

運用としては、timeとspatialの両方で改善しているときだけ「本物の改善」とみなして、Public LBはCVとのズレを確認する答え合わせくらいに留める。LBに合わせにいくと、だいたいPrivateで滑ります。

モデルは寄せ集めて平均を取る

単体で一番強いモデルを磨くより、毛色の違うものを混ぜたほうが安定しました。LightGBM・CatBoost・XGBoostをそれぞれ同じCVで回して、OOFでRMSEが最小になる重みを求めて足し合わせます。

def weighted_blend(results, names, y_true, mask):
    M = np.column_stack([results[n]["oof"] for n in names])[mask]
    loss = lambda w: rmse(y_true[mask], M @ w)
    cons = ({"type": "eq", "fun": lambda w: w.sum() - 1},)
    res = minimize(loss, np.full(len(names), 1/len(names)),
                   method="SLSQP", bounds=[(0, 1)]*len(names), constraints=cons)
    return dict(zip(names, res.x))

さらに後処理で、予測を少しだけPersistence側に寄せる（(1-k)*model + k*persistenceのkをOOFで決める）と、たまに出る大外しが抑えられてRMSEが落ち着きました。

スコア以外でも点が取れることを忘れない

冒頭に書いたとおり、このコンペは最終評価でTrustworthinessやInnovation & Practicalityも評価されます。ここはスコアだけ追っている人とちゃんと差がつくところなので、後回しにしないほうがいいと思っています。

具体的には、SHAPで「何が効いているか」を出して、緯度帯別・月別に「どこで外しているか」のマップを作って、乾燥域や半球で誤差が偏っていないかをチェックする。このあたりを最初からパイプラインに組み込んでおきました。

# 地域別・季節別のRMSEを出して、どこで外しているか可視化する
def grouped_rmse(df, y_true, y_pred, by):
    return (df.groupby(by)
              .apply(lambda g: rmse(g[y_true], g[y_pred]))
              .sort_values(ascending=False))

副産物として、弱い地域や季節が分かると次にどんな特徴量を足せばいいかのヒントにもなるので、純粋に精度改善のためにも役立ちます。

フォルダ構成

設定は全部 config/config.yaml に寄せてあるので、配布データの列名が違っても基本ここだけ直せば動くようにしています。

zindi-drought-tws/
├── README.md
├── requirements.txt              # 版数固定（再現性のため）
├── config/
│   └── config.yaml               # 列名・CV・target・モデルをまとめて管理
├── data/
│   ├── raw/                      # 配布データを置く場所
│   └── interim/  processed/  submissions/
├── notebooks/
│   └── 00_eda.py                 # 分布・欠損・季節×緯度・相関の確認
├── pipeline/
│   ├── 01_persistence_baseline.py   # 超えるべき基準を確認する
│   └── 02_train_and_submit.py       # 生データ→提出を一発で通す
├── src/
│   ├── config.py  utils.py  metrics.py
│   ├── data/load_data.py            # NetCDF(xarray)とCSVを結合
│   ├── features/
│   │   ├── climatology.py           # セル×月の気候値・偏差
│   │   ├── temporal.py              # ラグ・移動統計・トレンド
│   │   ├── spatial.py               # 緯度帯・空間クラスタ・近傍集約
│   │   ├── hydro.py                 # SPEI・土壌水分まわり
│   │   └── build_features.py        # 特徴量の統合とリーク防止、target定義
│   ├── validation/splits.py         # time / spatial / group_month / random
│   ├── models/train.py              # CVループと差分の復元
│   ├── ensemble/blend.py            # 重み最適化・スタッキング
│   ├── postprocess.py               # Persistenceブレンド・クリップ
│   └── explain.py                   # SHAP・誤差分析・バイアス監査
└── outputs/
    └── models/  figures/  oof/  logs/

config側はこんな感じで、ここを切り替えながら回しています。

target:
  mode: "delta"        # delta / anomaly / raw
validation:
  scheme: "time"       # time / spatial / group_month / random
  n_splits: 5
models:
  use: ["lgbm", "catboost", "xgb"]

動かし方

python -m venv .venv && source .venv/bin/activate
pip install -r requirements.txt

# data/raw/ に配布データを置く（列名が違えば config.yaml を直す）
python notebooks/00_eda.py                 # EDA（任意）
python pipeline/01_persistence_baseline.py # まず基準値を確認
python pipeline/02_train_and_submit.py     # 学習→アンサンブル→後処理→提出

進め方の順番

提出回数が限られているので、闇雲に投げないよう自分はこの順でやっています。

まずtime CVでdelta + LightGBM単体を回して、Persistenceをどれだけ下回れたかを見る。ここで一度だけ提出して、ローカルCVとLBがちゃんと連動しているかを確認します。連動が取れていれば、あとは特徴量を足し、CatBoostとXGBoostを加えてアンサンブルし、anomalyモードの提出も作ってdeltaと混ぜる。最後にspatial CVで未知地域の崩れを点検して、後処理のkを微調整する、という流れです。

最終提出に選ぶのは、LBの数字ではなくCVが良かった2本。これだけは崩さないようにしています。

おわりに

整理すると、効いたのは「TWSは直接当てずに差分を学習する」「セル×月の気候値で季節と地域の構造を渡す」「timeとspatialのCVでLBに釣られない」「毛色の違うモデルを混ぜる」の4つで、そこにTrustworthiness向けの説明性・再現性を足した形です。

正直まだ伸ばせる余地はあると思っているので、特徴量と検証はこれからも触っていく予定です。

それでは、Let's try!! で賞金ゲットだぜ。

補足

  • コンペ開催中にソリューションを公開していいかは規約を確認してください（終了前のコード共有が禁止されている場合があります）。
  • 記事にスコアを載せるときは、自分の環境で実際に出た数字に差し替えてください。

参考リンク

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?