RMSE を Fold ごとに取ると全体の値より小さくなる証明
この記事を書く前に twitter でお話をしている流れで、まますさんに的確な証明を頂くことができました! 証明にはこちら RMSE.pdf - Google ドライブ からアクセスできます。(まますさんありがとうございましたmm)
そもそも
この記事のお題は RMSE を Fold ごとに取ると全体の値より小さくなる証明をやります ということです。
これをやろうと思ったきっかけは #かぐるーど での kaggle本の本読みです。 前回は第5章だったのですが、その5.2.2で次のような記述があります。
クロスバリデーションでモデルの汎化性能を評価する際は、通常は各foldにおけるスコアを平均して行いますが、それぞれのfoldの目的変数と予測値を集めてデータ全体で計算する方法もあります。なお、評価指標によっては各foldのスコアの平均と、データ全体で目的変数と予測値から計算したスコアが一致しません。例えば、MAEやloglossではそれらが一致しますが、RMSEでは各foldのスコアの平均はデータ全体で計 算するより低くなります。
要するに K-Fold それぞれの rmse 平均値と、データセット全体での rmse の値だとデータセット全体のほうが大きい (K Fold のほうが良いように見積もられてしまう) という話です。たしかに Root をとる操作を毎回やるのと、全体で合わせた後やるのだと前者のほうが小さい値になような感じはしますよね。
これ一般的に示せるかなーという議論があり、僕が「関数の凸性とイェンゼンの不等式でいけますよ」と言ったところじゃあやってほしい!と言われたのが当エントリの経緯になります。
せっかくなので簡単にですが凸関数とイェンゼンの不等式にも触れつつ、お話できればと思っています。
NOTE: 若干細かい定義域についてやイェンゼンの不等式の導出についてなどは省略していますので、それらは別途文献など見ていただけば幸いです。
凸関数とは
下準備として、凸関数というのを定義します。凸関数というのは色々な定義がありますが、以下を満たすような関数 $f: \mathbb{R} \to (-\infty,+\infty]$ のことです
$$ f(t x_1 + (1 - t)x_2) \leq t f(x_1) + (1-t) f(x_2) $$
ただし $x_1, x_2$ は任意の実数 $\mathbb{R}$ の点で $t$ には 0以上1以下の制約がついています。
要するに $x_1$ と $x_2$ の内分点での $f$ の値と最初に $f$ で計算してしまってから $f(x_1)$ と $f(x_2)$ の内分を取るのとだと、後者のほうが大きいような関数、って言うことです。
また $f$ が微分可能な場合 $f$ の二階微分 $f'' \geq 0$ であることと上記の等式は同値になります。
イェンゼンの不等式
これをちょっと発展させて内分点の部分を2つ以上の点に拡張したのがイェンゼンの不等式です。
イェンゼンの不等式は上記の式と同じく特定の関数 $f$ が凸関数である必要十分条件を表した式で, $f$が凸ならば任意の自然数 $n$ と$\sum_{i=1}^n p_i = 1, p_i \geq 0$ を満たすような $p_i$ に対して次の式
$$ f(p_1 x_1 + p_2 x_2 \cdots + p_n x_n) \leq p_1 f(x_1) + p_2 f(x_2) \cdots + p_n f(x_n) $$
がなりたつ、という定理です。
たとえば$n=2$の時を考えてもらうと先ほどの凸関数の定義そのままであることはすぐわかると思いますので、凸関数の定義を変数 $n$ 個の場合に拡張したようなイメージです。
RMSE を考える
RMSE とは入力とラベルの誤差の2乗和を $M$ とした時に
$$ {\rm RMSE}(M) = M^{\frac{1}{2}} $$
で計算される値です。これの二階微分を考えると
$$ {\rm RMSE}''(x) = - \frac{1}{4} M^{- \frac{3}{2}} < 0 $$
です。すなわちRMSEの二階微分は常に負の値となります。これは凸関数と全く正反対の性質で一般に凹関数 (concave) と呼ばれ先のイェンゼンの不等式とちょうど不等号が反対の不等式が成立します。
K-Foldしたときの RMSE
今Fold を$ K$ 個に分割して、それぞれが $n_k$ 個のデータを持っているとします。(データセット全体では $N$ 個とします。) この時各 Fold での MSE (Mean Squared Error) を $M_k$ とすると Fold ごとのデータの数で重みづけた ${\rm RMSE}_{\rm fold}$ は
$$ {\rm RMSE}_{\rm fold} = \sum_{k=1}^K \frac{n_k}{N} \sqrt{M_k} $$
となります。一方で通常の RMSE に関しては
$$ {\rm RMSE} = \sqrt{\frac{1}{N} \sum_{k=1}^K n_k M_k} = \sqrt{\sum_{k=1}^K \frac{n_k}{N} M_k} $$
となります。ここで $M_k$ に $n_k$ をかけているのは $M_k$ が既に Mean Squared Error なので要素の数を掛けて和にになおして全体の $N$ で割算をするためです。
ここで $p_k = n_k / N$, $f(x) = \sqrt{x}$ と考えると $f$ は凹関数でかつ $\sum p_k = 1$ ですのでイェンゼンの不等式が用いることが出来て
$$ {\rm RMSE} = \sqrt{\sum_{k=1}^K \frac{n_k}{N} M_k} \geq \sum_{k=1}^K \frac{n_k}{N} \sqrt{M_k} = {\rm RMSE}_{\rm fold} $$
が成立します。即ち fold ごとで RMSE を計算して重み付きの平均を取った値のほうが、データセット全体での RMSE の値より小さくなることがわかりました。
RMSE 以外でも…
上記の証明を追っていただくと分かるようにこの証明はロス関数の値がデータごとに計算できること、及びそれをデータセット全体の平均したあとに凹関数に代入する、という構造が保たれている限り同様の議論をすることが可能です。
ですので Log を取ってから Root を取る RMSLE (Root Mean Squared Log Error) なども同様の議論が可能です。
参考文献
以下は本記事を書くにあたって使用した凸関数に関する話題や凸最適化に関する日本語の参考文献です。
- wikipedia: 凸関数: 凸関数 - Wikipedia
- 非線形計画法: 山下 信雄
- 非線形最適化の基礎: 福島 雅夫
- 工学基礎 最適化とその応用: 矢部 博
この記事は kaggle その2 advent calendar 2019 の記事です。
分析コンペLT会でLTをさせてもらいました!!
2019/11/30 に行われた分析コンペLT会にLT枠として参加させていただきました。😄
僕は普段大阪で仕事をしているのでもともと発表はおろか参加する予定はなかったのですが(行きたいとはめっちゃおもっていた)、いつもかぐるーどでお世話になっているカレーさんに「発表枠あけて待ってます！！」(原文ママ)と嬉しいオファーを頂いき、ちょっと東京旅行兼ねて発表枠として参加しました。
発表者枠あけて待ってます！！
— カレー🍛専業kaggler (@currypurin) 2019年11月8日
内容に関しては俵さんが素敵にまとめて頂いているのでそちらを参考にしていただければと思います。
発表内容
僕の発表は「初手が爆速になるフレームワークを作ってコンペ設計した話」というタイトルで、自分が作っている https://gitlab.com/nyker510/vivid という機械学習用のフレームワークとそれを使ってコンペ設計が楽になったよ、ということを話しました。
Vivid について
Kaggle でもそうですがお仕事でも当然特徴量とモデルのバージョンを管理するのはとても大切です。なんですが僕はとても大雑把な人間なので、良く「このファイルどのスクリプト or Notebook から出てきたんだっけ…」ということに陥っていました。
そこで「動かすだけで勝手にログもモデルもバージョンも保存してくれるやつがあったらいいなーというかないと僕は無理だな」と思い、自分で色々と試行錯誤をして出来上がったフレームワークです。
大きなコンセプトというか特徴は以下のような感じ。
- 必要なことだけを書くので良い
- 基本は勝手にやってくれる
- k-fold の split をして oof を計算するコード、みたいな定形処理は全部 vivid にお願いして、プロジェクト固有のコードに集中できる用に
- ログの出力やモデルの保存なんかも勝手にやってくれる用に
- テンプレート的な特徴作成の提供
- 毎回 count encoding のコードを書くのは良くないのでそれもやってもらう
- 特徴量をある粒度 (atom とよんでいます) とその集合体 molecule で管理して versioning する機能とかもあります
- スタッキング・アンサンブル対応
- 対応というよりは、一気通貫に出来るというのが売りです。
(全部一気につながって作るので、前の run で作った特徴で学習していて予測するとバージョン違いで精度が出ない・カラムの数が違う、といった悲しいミスを防ぐ)
- 対応というよりは、一気通貫に出来るというのが売りです。
立ち位置としては sickit-learn よりも更に盛ったフレームワークという感じでしょうか。(webのDjango的な)
オレオレフレームワークに過ぎないのでコードもアレだしドキュメントもないしで正直めちゃ恥ずかしいんですが「こうやったらいいんじゃないか」とか「僕はこうやってますー」みたいな知見がでてお話できればとてもうれしいです。;)
ソースコードはこちら https://gitlab.com/nyker510/vivid からアクセスできます。 pip なら
pip install git+https://gitlab.com/nyker510/vivid
でインストールできますのでちょっと触っていただけるとめちゃ喜びます
atmaCupについて
スライドで atmaCup について触れたのですが結構な方に知っていただいてとても嬉しかったです! 次回第3回も鋭意企画中ですのでぜひご参加いただけるとこちらもとてもとても喜びます😆
感想
やはり発表すると思っていることがまとまるので、発表者が得るものが大きいなととても感じました。今回もこの発表のために vivid をちょっと直して(そしてバグを見つけ😌)、スライドにするために自分の考えを一度整理して、という作業をおこなったのでいま自分がやっていることの意味合いがクリアになって非常に良い経験でした。
また発表のことについて質問してディスカッションできたり、僕も他の人の発表について聞けて(そしてどれもめちゃくちゃクオリティが高い!!) 質問できたりというのは実際にあって喋る良い点だなあと思いました。
今後も機会があれば是非参加したいと思います! 企画頂いたかぐるーどの皆様、会場提供頂いた日経BP様、関係者の方々本当にありがとうございましたmm
PRML の本読みをしています @section3
最近(というか今日から) 会社でPRML勉強会をやっています。ふつう第1章からやるのが普通ですがPRMLはちょっと重たいので息切れすると良くないよねということでいきなり第3章から始めるという方針をとっています。*1
今回は僕が担当で、主に
- 線形モデルの導入
- 正則化 (なんで Lasso はスパースになるのかの話)
- バイアス・バリアンス
の話をしていました。たぶんこの話の中で一番一般的に使われるのは「バイアスバリアンス」の話だと思っていて、改めてまとめてみて個人的に良かったです。ざっと今日話した内容含めてまとめると以下のような感じかな?と思ってます。
前提
- $D$: データの集合
- $E_D$ データ集合での期待値を取ったもの
- $y$: 予測関数
- $E_D[ y(x;D)]$ データの集合で期待値を取った予測関数。いろんなデータでモデルを作ってアンサンブルしてるイメージ
- $h$: 理想的な予測値
バイアス (モデルの傾向)
- データを沢山取ってアンサンブル平均した $E_D[ y(x;D)]$ が理想的な関数 $h$ にどれだけ近いか
- ロジック自体の傾向が強い(biasが強い)と大きくなる
- [低] 表現力が大きく多様性があるモデル
- [高] 表現力が小さいモデル (いくら平均しても理想に近づけない)
- たとえば常に 1 を返すようなモデル $y(x) = 1$ などはハイバイアスです。
- 癖が強い(データに合わせる気がない)というイメージ
バリアンス (Variance: モデルの分散)
- 今回のデータでの予測値と平均的なデータでの予測値の差分
- 毎回のデータでどれぐらいばらつくか
おまけ
バイアスバリアンスの説明に使われている Figure 3.5 にあたる絵を書く python script を書いてみました。
自分で書くと色々と考えるので楽しいですね。
import numpy as np import matplotlib.pyplot as plt def gaussian_kernel(x, basis=None): if basis is None: basis = np.linspace(-1.2, 1.2, 101) # parameter is my choice >_< phi = np.exp(- (x.reshape(-1, 1) - basis) ** 2 * 250) # add bias basis phi = np.hstack([phi, np.ones_like(phi[:, 0]).reshape(-1, 1)]) return phi def estimate_ml_weight(x, t, lam, xx): basis = np.linspace(0, 1, 24) phi = gaussian_kernel(x, basis=basis) w_ml = np.linalg.inv(phi.T.dot(phi) + lam * np.eye(len(basis) + 1)).dot(phi.T).dot(t) # bias があるので +1 しています xx_phi = gaussian_kernel(xx, basis=basis) pred = xx_phi.dot(w_ml) return pred n_samples = 100 fig, axes = plt.subplots(ncols=2, nrows=3, figsize=(10, 12), sharey=True, sharex=True) for i, l in enumerate([2.6, -.31, -2.4]): ax = axes[i] preds = [] for n in range(n_samples): x = np.random.uniform(0, 1, 40) xx = np.linspace(0, 1, 101) t = np.sin(x * 2 * np.pi) + .2 * np.random.normal(size=len(x)) pred = estimate_ml_weight(x, t, lam=np.exp(l), xx=xx) if n < 20: ax[0].plot(xx, pred, c='black', alpha=.8, linewidth=1) preds.append(pred) ax[1].plot(xx, np.sin(2 * xx * np.pi), c='black', label=f'Lambda = {l}') ax[1].plot(xx, np.mean(preds, axis=0), '--', c='black') ax[1].legend() fig.tight_layout() fig.savefig('bias_variance.png', dpi=120)
*1:若干荒業感がありますが、一人以外は一度は読んだことがあるというのとやはり息切れが怖いのでちょっと先に応用が出てきそうなところから取り組んでいます
Nelder Mead Method を python で実装
最近会社で Kaggle 本 Kaggleで勝つデータ分析の技術 の本読みをやっています。昨日は第二章でしきい値の最適化の文脈で Nelder Mead Method を使った最適化が記載されていました。ちょっと気になったので実装してみたという話です。
Nelder Mead Method とは
最適化手法の一種です。いくつかのデータ点(単体)での評価値をベースに更新していく方法で、勾配の情報などを使わないため目的関数さえ評価できれば実行できるお手軽さがメリットです。
アルゴリズムの概要は英語版 wikipedia が詳しいです。(自分が見た限りでは日本語版では一部の更新で符号がまちがっているように見受けられたので英語のほうがよさ気)
前提として $N$ 次元のベクトルの最適化をしたいとすると
- 適当に選んだ $N$ 次元ベクトル + 各次元に $+1$ した合計で $N+1$ の $N$ 次元ベクトルを計算
- 1で計算した各点での目的関数値を計算して、その順序を保存。
- 一番悪かった点以外の重心を計算
- 一番悪かった点と重心の対称な点(反射点)を計算
この反射点を使って場合分けを計算していきます。
- 反射点が二番目に悪い点より良い場合は採用として一番悪い点を置き換え
- 反射点が一番良い場合、もっとソッチの方に移動しても良いという発想で、最悪点から反対側に更に移動した点(拡大点)を計算。
拡大点が反射点より良い場合には拡大点を採用として一番悪い点を置き換えます。(拡大点が反射点より悪い時はすなおに反射点で妥協) - それ以外の時(2番めのものより悪い時)は重心と最悪点の内点(縮小点)を計算
- 縮小点が最悪よりは良い時は最悪点を縮小点で置き換え。
- 縮小点も悪い時は最良点以外の点を、最良点との内点に置き換え
python のコードにすると以下のような感じです。Notebook はこちら https://nyker510.gitlab.io/notebooks/posts/nelde-mead/ にあげています。
import numpy as np init_point = np.random.uniform(-.1, .1, 2) - np.array([1.5, 1.5]) alpha = 1. gamma = 2. rho = .5 sigma = .5 scaling_factor = 1. eval_points = np.vstack([init_point + np.eye(2) * scaling_factor, init_point]) best_points = [init_point] for i in range(100): # 目的関数で評価して並び替え objective_values_by_point = objective(eval_points[:, 0], eval_points[:, 1]) orders = np.argsort(objective_values_by_point) f_values = objective_values_by_point[orders] eval_points = eval_points[orders] # 最悪点以外で重心を計算 centroid = np.mean(eval_points[:-1], axis=0) p_worst = eval_points[-1] # 反射点 (最悪点から重心の対象な点) を計算 + そのときの目的関数を評価 p_reflect = centroid + alpha * (centroid - p_worst) f_ref = objective(*p_reflect) f_best = f_values[0] f_second_worse = f_values[-2] f_worst = f_values[-1] # 反射点がそこそこ良い場合 (second worse よりは良い) → 最悪点を反射点で置き換え if f_best <= f_ref < f_second_worse: eval_points[-1] = p_reflect # 反射点が一番良い場合 elif f_ref < f_best: # もうちょっと良い点がないか探りに行くイメージ (拡大点) # 最悪点が一番悪い → 反射点の方にもっと移動するイメージ (反射点 + centroid - p_worst [i.e. 最悪点からみた重心方向]) p_expand = centroid + gamma * (centroid - p_worst) f_expand = objective(*p_expand) # 拡大点が良かったら採用. だめだったら反射点で妥協 if f_expand < f_ref: eval_points[-1] = p_expand else: eval_points[-1] = p_reflect # 反射点が二番目に悪い点より良くない時 elif f_second_worse <= f_ref: # 重心と最悪点の間を計算 (大体すべての点のまんなかあたりに対応する) p_cont = centroid + rho * (p_worst - centroid) f_cont = objective(*p_cont) # 縮小点が最悪よりは良い時 # 重心方向への移動は良いと判断して最悪点を置き換え if f_cont < f_worst: eval_points[-1] = p_cont else: # そうでない時最良点めがけて縮小させる (最良点との内点に移動して三角形ちっちゃくするイメージ) shrink_points = eval_points[0] + sigma * (eval_points - eval_points[0]) eval_points[1:] = shrink_points[1:] best_points = np.vstack([best_points, eval_points[0]])
実験
はじめに自分で適当に作った関数でやります。たぶん (0, -.4) あたりが最適解なはず(あやしい)。
def objective(x, y): return (x - .1) ** 2 + (y - .3) ** 2 + np.sin((x + .5) * 2) * np.sin(y * 3) + sigmoid(x * 4)
これでやると以下のような感じ。わりとよさ気。
ついでに Rosenbrock function でもやってみましょう。
こちらは最適な座標は (1, 1) になるように設定しています(見えにくいですが黒い丸があるところ)。
実行すると以下のような感じ。たどり着いてるっぽいですね。
NOTE
これは Nelder Mead に限った話ではないですが、問題のスケールをちゃんと -1, 1 ぐらいに合わせておく必要があります。というのは一番初めに生成する点列の評価値に依存してその後の点列が決定するため、一番初めの点の性質が悪いと上手く更新が進まないためです。
例えば今回で言うと init_point に対して np.eye(1) を足すことで初期の点を生成していますが、これを np.eye(1) * .00001 などにすると上手く収束してくれません。
反対に大きすぎても更新幅が大きすぎて発散気味になります。
これと同様にそもそもの初期点をどう選ぶかもとても重要です。コード中では np.random.uniform(-.1, .1, 2) - np.array([1.5, 1.5]) としていますが、たとえばこの周りだけ凹んでいる関数だった場合そこから抜け出せずに局所解に陥って終わりになることもあります。ちょっとうまくいかないなーという時はスケールがあっているかを確認+初期点を変えてやってみるのも大切そうです。
