nykergoto’s blog

機械学習とpythonをメインに

JWT でのログインができるかどうかを確認するテスト

Django Rest Framework + allauth 使用時に JWT での Login ができているかどうかを念の為確認するテスト。

import pytest
from django.urls import reverse

from atma.shops.models import Shop
from .factories import UserFactory

@pytest.mark.django_db
def test_authorized_by_jwt(client):
    """Json Web Token で login できる"""
    user = UserFactory()  # type: User
    user.set_password('hogehoge')
    user.save()

    url = reverse('accounts_me')
    res = client.get(url)
    assert res.status_code == 401, res.json()

    auth_url = reverse('rest_login')
    res = client.post(auth_url, data={'email': user.email, 'password': 'hogehoge'})
    assert res.status_code == 200, res.json()

    jwt_token = res.json().get('accessToken', None)
    assert jwt_token is not None

    header = {'HTTP_AUTHORIZATION': 'JWT {}'.format(jwt_token)}
    res = client.get(url, **header)

    assert res.status_code == 200, res.json()
  • パスワードの設定は user.set_password で行なう
  • header を kwrgs で渡す
  • "HTTP_AUTHORIZATION" を key にする

の3点がポイント

NGBoostを読んで、実装する。

不確実性を考慮した予測が可能と噂の NGBoost の論文を読んでみたので、全体のながれをまとめて見ました。加えて自分でも NGBoost を実装して、その結果を載せています。

元の論文 NGBoost: Natural Gradient Boosting for Probabilistic Prediction はこちら https://arxiv.org/abs/1910.03225

Introduction

一般的な教師あり学習を考えます。このとき予測モデルは入力データ $X$ に対して予測値 $y$ を出力するように学習していきますが、たいていのモデルではひとつのデータに対しては予測値はひとつしか得ることができません。

例えばウェブ上の行動履歴から、ユーザーの年齢を予測してください、という問題があったとすると、ユーザーの期待される年齢そのものを返すようなモデルがそれに当たります。

例えばテーブルデータでデファクトスタンダードとなっている XGBoost や LightGBM のような勾配ブースティングもそうですし一般的な Deep Learning もそのような設計です。

しかし実際には予測には確信度のようなものが存在しているはずです。ユーザの年齢だけ答えてくださいと言われると数字を出力するしかありませんが、実際には「このユーザーは絶対に若い」とわかるようなデータだったり「年取っていると思うけれど自信ない」というようなデータも存在しているでしょう。 このようなデータによる不確実性も同時にモデリングしましょうというのがモチベーションです。

f:id:dette:20200501012357p:plain
NGBoost Weibiste https://stanfordmlgroup.github.io/projects/ngboost/ より。あるデータに対する予測はふつう一つしか得られない(点推定)。それでは予測の確実性がわからないため確率分布(右)で出力したいよねというモチベーションです。

不確実性を考慮した予測

一般的なモデリングでは確率分布それ自体ではなくて、代表するようなパラメータを予測します。 例えば RMSE を目的関数としている場合では平均値 $\mu$ を予測する関数を作ることが目的です。

不確実性を考慮する予測とは入力 $X \in \mathbb{R}^k$ にたいして対応する予測値 $y$ がどのような分布になっているかを関数として予測する必要があります。すなわち

$$ P(y | X) $$

で表される確率分布 $P$ を推論する作業です。 このままだと分布 $P$ は無限の形状をとれて埒が明かないので特定のパラメータ $\theta \in \mathbb{R}^p$で表現できるような分布 $P_{\theta}$ に絞って考えていくことにします。 例えば $P$ に正規分布を仮定して正規分布は平均値 $\mu$ と分散 $\sigma$ の2つで説明できる分布なので $\theta = (\mu, \sigma)$ とおく、というような操作に相当します。

この仮定によって分布自体を考える必要はなくなり「パラメータ $\theta$ をいかに予測するか」という問題に変えることができます。 以下では $X$ から分布のパラメータ $\theta$ を出力するような関数 $F$ を推論することを考えていきます。

$\theta$ の推論のための評価基準

$\theta$ を推論するためには $\theta$ のクオリティを評価できる必要があります。すなわち $\theta$ とそのときの実現値 $y$ を入力として数値を返すような関数 $S$ を決めなくてはいけません(以下では $S$ が小さいほどよい関数だとします)。 最も良い $\theta$ とは何かを考えると、これは間違いなく $P_{\theta}$ が真の分布 $Q$ に一致するような場合ですね。理想的には。 これを言い換えると $\theta$ を評価する関数 $S$ は $P$ が真の分布 $Q$ になったとき、もっとも小さくなっている必要があります。これを数学的に記述すると

$$ \mathbb{E}_{y \sim Q} [S (Q, y)] \le \mathbb{E}_{y \sim Q} [S(P, y)] \ \ \forall P, Q $$

となります。ここで $y \sim Q$ は確率分布 $Q$ からのサンプル $y$ での期待値を表します。すなわち

$$ \mathbb{E}_{y \sim Q} [S (Q, y)] = \int S(Q, y) Q(y) dy $$

です。このような性質を満たしている指標はいくつかあり最も有名なものが負の対数尤度関数 $\mathcal{L}$ (MLE) です。 $$ L(P, y) = - \ln P(y) $$

今回は $P$ がパラメータ $\theta$ で表わされるものと仮定したので実際には $\theta$ の関数となりますので以下のように記述できます。

$$ \mathcal{L}(\theta, y) = - \ln P_{\theta}(y) $$

この他にも CRPS (以下では $\mathcal{C}$ と表します) という指標もありこちらは以下のように記述される値です。 (一般に MLE よりもロバストであるとのことです)

$$ \mathcal{C} = \int_{- \infty}^{y} F_{\theta} (z)^2 dz + \int_{y}^{\infty} (1 - F_{\theta}(z))^2 dz $$

ここで $F_{\theta}$ は分布 $P_{\theta}$ の累積分布関数です。 (予測関数 $F$ と同じ記号ですけれど違います。念の為。)

$\theta$ の最適化

これらのスコア $S$ があればパラメータが最適化できそうです。ここでは一般的な勾配法を使った場合を考えてみます。 $S$ に対しての $\theta$ に対する最急降下方向 $\nabla S$ は以下のように記述できます。

$$ \nabla S (\theta, y) \sim \lim_{\epsilon \to \infty} \arg \min_{d: | d | = \epsilon} S (\theta + d, y) $$

一見これで問題なさそうに見えますが $S$ はもともと分布 $P$ に対して計算されるものですからパラメータ空間 $|d|$ でのユークリッド距離を用いた範囲内で min をとると $S$ による歪みを考慮することができません。

ここで $S$ が満たす必要がある式を思い出して右辺から左辺を引くと

$$ \mathbb{E}_{y \sim Q} [S (P, y)] - \mathbb{E}_{y \sim Q} [S (P, y)] \ge 0 $$

となります。これは非負の値ですから $Q$ から $P$ がどのぐらい離れているかの度合い (距離) とみなすことができます。 例えば負の対数尤度であれば KL Divergence に相当し, CRPS であれば $L^2$ Divegence となります。

これらの Divergence は Q / P がどのようなパラメータを持っているかに不偏です (Divergenceの式自体には $\theta$ が出現しないことに注意してください)。したがってこの Divergence で表されるような空間で議論することができれば $\theta$ ではなく分布 $P$ の距離の意味で最も小さくなる方向を見つけることができるため、パラメータ $\theta$ のとり方に依存しない方向を見つけることができます。

すなわち最急降下方向を求める際の移動量を測る空間であったユークリッド距離 $| d |$ を、もともとのパラメータ $\theta$ と移動した点 $\theta +d$ での確率分布 $P_{\theta}, P_{\theta + d}$ との距離を測る Divergence $D_S$ で置き換えて

$$ \nabla \hat{S} (\theta, y) \sim \lim_{\epsilon \to \infty} \operatorname{argmin}_{d: D_S (P_{\theta} | P_{\theta + d}) = \epsilon} S (\theta + d, y) $$

とします。こうして求められる勾配はパラメータに不偏です。これを自然勾配 (Natural Gradient) と呼びます。 これを解くとスコア $S$ を用いたときの点 $\theta$ でのリーマン計量行列 $\mathcal{I}_{S}(\theta)$ としたとき

$$ \nabla \hat{S} (\theta, y) \propto \mathcal{I}_{S} (\theta)^{-1} \nabla S(\theta, y) $$

と記述できます。例えば $S = \mathcal{L}$ (負の対数尤度) とするとこれはフィッシャー情報量行列です。

$$ \mathcal{I}_{\mathcal{L}} (\theta) = \mathbb{E}_{y \sim P_{\theta}} [ \nabla_\theta \mathcal{L}(\theta, y) \nabla_\theta \mathcal{L}(\theta, y) ^T] $$

自然勾配を使うとパラメータ変換に堅牢に最適化を行えます。 Figure3 では正規分布に対して平均値 $\mu$ と分散を対数変換した $\log \sigma$ 場合の勾配の比較があります。これを見ると通常の勾配 (右側) よりも自然勾配 (左側) のほうが安定して最適解への方向をしめしていることがわかります。

f:id:dette:20200501013136p:plain
Figure3. 左の列が通常の最急降下方向を矢印で表したもので、右が自然勾配。最急降下方向よりも最適解に近い方向への良い勾配となっていることがわかる。

$F$ の最適化

では次に $\theta$ を求めるような関数 $F$ を作っていくのですが NGBoost ではこれを Natural Gradient Boosting で学習します。 Natural Gradient Boosting は Gradient Boosting をちょっといじったものなので、まずは Gradient Boosting から。

Gradient Boosting

Gradient Boosting は求めるモデル $F$ が複数のモデル $f_i$ の足し算 $F_{m} = \sum_{i=1}^m f_i$ で表される加法的モデルと仮定します。 これを $m$ 個同時に最適化することは難しいので $m$ 番目のステップではすでに作った $i = { 1, \cdots, m-1 }$ のモデルは固定して、次に加えるモデルをスコア $S (\theta, y)$ の $\theta$ に対する勾配と一致するように学習します。すなわち $i$ 番目のデータに対する予測結果が

$$ f_{m} (x_i) = - \frac{S ( F_{m-1}(x_i), y_i)}{\partial F_{m-1}(x_i)} $$

となるよう $f_m$ を求めます。 ($m-1$番目までの状態で得られる $\theta_{i}^{m-1}$ は $F_{m-1}$による予測値 $F_{m-1}(x_i)$ に一致していることに注意してください.) 具体的には $g_m^i = - {S ( F_{m-1}(x_i), y_i)} / {\partial F_{m-1}(x_i)}$ とおいて二乗誤差

$$ \sum_{i=1}^N \left( g_m^i - f_{m}(x_i) \right)^2 $$

を最小化するように $f_m$ を決定します。このとき $f$ に使うモデルを弱学習器と呼び、一般には深さを制限した決定木などが使われます。

Natural Gradient Boosting

Gradient Boosting 目標としていた再急降下方向を自然勾配に置き換えてしまおう、というのが Natural Gradient Boosting です。すなわち $f_m$ を近づける勾配を自然勾配に置き換えます。要するにリーマン計量行列の逆行列で補正した

$$ \hat{g}_m^i = - \mathcal{I}_{S}(F_{m-1}(x_i))^{-1} \frac{S ( F_{m-1}(x_i), y_i)}{\partial F_{m-1}(x_i)} $$

を使って

$$ \sum_{i=1}^N \left( \hat{g}_m^i - f_{m}(x_i) \right)^2 $$

を最小化するように $f_m$ を決定します。

アルゴリズム

それでは実際に NGNatural Gradient Boosting のアルゴリズムがどのように動いていくかをまとめてみます。 上記でのべた仮定をふまえるとアルゴリズムを動かすためには初めに

  • 予測対象とする分布 $P_{\theta}$
  • 各 Boosting で学習するためのモデル $f$
  • $\theta$ の良さを評価するための指標 $S$

の3つを決めておく必要があります。 これらを使って、以下の 1 ~ 4 を $M$ 回繰り返すことで学習器とステップサイズ $\{ \rho^m, f^m \}_{m=1}^M$ の集合を作っていきます。 (以下の $m$ は現在の boosting の iteration 数です)

f:id:dette:20200501013409p:plain
NGBoostのアルゴリズム擬似コード

1. 探索方向 (Natural Gradient) の計算

今のパラメータ $\theta^{(m-1)}$ を用いて $i$ 番目のデータに対する自然勾配 $g_i^m$ を計算します。

2. $f_m$ の学習

次にデータ $x_i$ から勾配 $g_i^m$ を予測するモデルを作成します

$$ f^m = \arg \min_{f} \sum_{i=1}^N \left( f(x_i) - g_i^m \right)^2 $$

このとき求めているのは予測値 $y$ に関係するものでなく、単にパラメータ $\theta$ の勾配です。したがって学習記モデルは基本的に Regression を行なうものなら何でも OK です。

3. ステップサイズ $\rho$ の決定

モデル $f^m$ が決定すれば各データ $x_i$ ごとにパラメータ $\theta_i$ がどちらに移動した良いかを $f(x_i)$ によって推論することができます。 (実際に次のステップ 4 ではこの $f(x_i)$ によって $\theta_i^m$ を更新します。) そのままこの更新方法を使っても良いですが、最適解からとても遠い (or 近い) ことを考慮して、勾配を定数倍して最もスコアが下がる点を探します。

$$ \rho^m = \arg \min_{\rho} \sum_{i=1}^N S\left(\theta^{(m-1)} - \rho f^{m} (x_i) , y_i \right) $$

実際には linear search によって探索し最適なステップサイズを $\rho^m$ とします。

4. パラメータ $\theta$ の更新

ステップサイズ $\rho$ とモデルによって予測された更新方向 $f(x_i)$ を使って実際にパラメータを更新します。このとき値をそのまま使わず定数 $\eta$ をかけて縮めます。

$$ \theta_{i}^{m} = \theta^{(m-1)} - \eta \cdot \rho^m f^m(x_i) $$

(Gradient Boosting などの文脈では learning rate と呼ばれることもあるパラメータです。一気に学習が進んでしまうのを防ぐのが狙いです。)

NGBoost の利点

論文中にあげられている NGBoost の利点は大きく以下の3つです。

1. いろんな分布で使える

予測対象としている分布の制限は単にパラメータで表現できるというものだけですから、いろんなケースで使用することができます。(たとえば負の二項分布だったり、ガンマ分布だったり)

2. 複数のパラメータの同時推定でも最適解を見つけやすい

通常の勾配を用いると初期値の時点でたまたまうまく予測できている lucky なデータがある場合それらのデータが支配的になってしまい、分散推定がそれらのたまたまうまく予測できたデータのみ進んでしまい、反対にたまたま初期値が大きく外れている unlucky なデータで学習が足りないという現象が発生します。 (Figure4) 自然勾配を使うとフィッシャー情報量行列などで勾配を補正するため、初期点にせずにまんべんなく学習をすすめることができます。

f:id:dette:20200501013503p:plain
Figure4. 通常の勾配法(上段)をつかっていると学習が進んでも平均値は変わらずに分散ばかり学習してしまっている。自然勾配法(下段)だと平均値・分散ともにまんべんなく学習が進んでいる。

3. パラメータ変換に不偏

自然勾配の定義からもわかるようにパラメータ変換に依存せずに勾配を選ぶことができます。これによってユーザーは推定する分布のパラメータを好きに変数変換しても問題ありません。例えば正規分布が対象であれば予測するために平均値と分散が必要になりますが、分散は正である必要があるためログ変換したパラメータ $\log \theta$ を推論するほうが数値的にも安定して便利です。
このとき自然勾配を使っていれば問題なく最適化できますが普通のニュートン法だともはやこの問題は凸ではないため解くことが難しくなります。

実験と結果

table1 では不確実性を考慮できるいくつかのモデルと NGBoost の性能比較がされています。 評価指標は負の対数尤度 (NLL) です。 これを見ると負けているデータセットもあるものの概ね良いスコアを出せていることがわかります。

table2 では Natural Gradient を使わない通常の Gradient Boosting (Multiparameter Boosting) と、分散を固定して不確実性を考慮しない Boosting (Homoscedatic Boosting) との比較がされています。 これも評価指標は負の対数尤度です。 こちらを見ると NGBoost は Gradient Boosting の場合よりも性能が高いことがわかります(Natural Gradient にした成果)。 また分散を固定したモデルとも匹敵した性能がでていて分散を個別に推定できるという利点を考えると NGBoost を使って損はない、という結果になっています。

f:id:dette:20200501013832p:plain
Table1 & Table2. 評価指標は負の対数尤度

負の対数尤度では分布全体を考慮する NGBoost が有利でしょう、ということで点推定の結果を出したものが table3 です。こちらは評価指標が RMSE となっています。 NGBoost 以外のモデルは RMSE の点推定を目的関数に持っており NGBoost は違うにもかかわらず、他のモデルに匹敵した性能がでていることがわかります。

f:id:dette:20200501013912p:plain
Table3. 点推定値に対してのRMSE。直接RMSEを最適化していないのにもかかわらず NGBoost はその他のモデルに匹敵する性能を出していることがわかる。

結論

  • NGBoost は不確実性を考慮したモデリングを出来るように Gradient Boosting を拡張したモデルです。
  • 安定した探索方向を作成するため Gradient Boosting の各ステップでの最急降下方向を Natural Gradient に変更しました。
  • 基本的には Gradient Boosting と同じなのでモンテカルロ法ベイズの知識がなくとも簡単に不確実性を考慮した予測をすることが可能です。
  • 点推定モデルとくらべて、点推定タスクに置いても遜色ない性能が見込まれるため「コスト無しで (for free)」不確実性を推論できます。

実装

ここまでは論文と調べた内容のまとめです。 実際に NGBoost を使う場合 python package が公式から提供されています https://github.com/stanfordmlgroup/ngboost 。 ライブラリ自体は pip から install できますのでさっくり試してみたい方はこちらを使うことをおすすめします。 User Guide もあり自分で Score $S$ や分布 $P$ を カスタムしたい場合のガイド もあります。

pip install --upgrade git+https://github.com/stanfordmlgroup/ngboost.git

以下は勉強兼ねて、僕が自分で実装してみたよという内容です。 コード自体はちょっと違いますが script にまとめたものは以下から参照ください。

github.com

自然勾配の計算

アルゴリズムの部分にもあったようにコアになるのは現状のパラメータから Natural Gradient を返す部分です。 対数尤度をスコアとしている場合勾配とフィッシャー情報量行列を計算する必要がありますが、上記の議論でもあったように、もとめるパラメータの変数をどのように定義するか(変数変換の定義方法)によってこれらの形式も変動します。

今回は正規分布を予測分布として平均値 m, 分散をログ変換した log_var を推論対象としました。 このとき勾配とフィッシャー情報量行列は以下のように計算できます。*1

var =  np.exp(log_var) # log -> var に変換
gradients = [
    (m - t) / var,
    1 / 2 - ((t - m) ** 2) / (2 * var)
]

fisher_matrix = [
    [1 / var, 0],
    [0, 1 / 2]
]

弱学習器としては決定木 sklearn.tree.DecisionTreeRegressor を用います。 このときあまり精緻に木を作ってしまうとオーバーフィットのもとになるため max_depth / min_sample_leaf だけパラメータ指定をしています。

from sklearn.tree import DecisionTreeRegressor

def build_weak_learner():
    return DecisionTreeRegressor(max_depth=3, min_samples_leaf=4)

学習する

あとは Boosting 部分を実装すればOKです。 全体を記述すると以下のようになります。

import numpy as np
from scipy.stats import norm
from sklearn.tree import DecisionTreeRegressor

def build_weak_learner():
    return DecisionTreeRegressor(max_depth=3, min_samples_leaf=4)

class LogVarianceNorm:
    def calculate_gradient(self, t, m, log_var):
        var =  np.exp(log_var)
        gradients = [
            (m - t) / var,
            1 - ((t - m) ** 2) / var
        ]

        fisher_matrix = [
            [1 / var, 0],
            [0, 1]
        ]

        return np.array(gradients), np.array(fisher_matrix)
    
    def score(self, t, m, log_var):
        var = np.exp(log_var)
        return - norm.logpdf(t, loc=m, scale=var ** .5)

def true_function(X):
    return np.sin(3 * X)

n_samples = 200
np.random.seed(71)
X = np.random.uniform(-3, 1, n_samples)
y = true_function(X) + np.random.normal(scale=abs(X), size=n_samples)

log_var_norm = LogVarianceNorm()

# hyper parameters
learning_rate = 1e-2
n_iterations = 500

estimators = []
step_sizes = []

X = X.reshape(-1, 1)
N = len(X)

# initialize parameters by whole dataset mean and variance
m, log_var = np.zeros(shape=(N,)), np.zeros(shape=(N,))
m += np.mean(y)
log_var += np.log(np.var(y) ** .5)

for iteration in range(n_iterations):

    # 1. compute natural gradient
    grads = []
    for i in range(N):
        score, fisher = log_var_norm.calculate_gradient(y[i], m[i], log_var[i])
        grad = np.linalg.solve(fisher, score)
        grads.append(grad)
    grads = np.array(grads)

    # 2. fit weak learner
    # mean estimator
    clf_mean = build_weak_learner()
    clf_mean.fit(X, y=grads[:, 0])

    # log_var estimator
    clf_var = build_weak_learner()
    clf_var.fit(X, y=grads[:, 1])

    directions = np.zeros(shape=(N, 2))
    directions[:, 0] = clf_mean.predict(X)
    directions[:, 1] = clf_var.predict(X) 
    estimators.append(
        [clf_mean, clf_var]
    )

    # 3. linear search
    scores = []
    stepsize_choices = np.linspace(.5, 2, 21)
    for stepsize in stepsize_choices:
        d = directions * stepsize
        score_i = log_sigma_norm.score(y, m - d[:, 0], log_var - d[:, 1]).mean()
        scores.append(score_i)

    best_idx = np.argmin(scores)
    rho_i = stepsize_choices[best_idx]
    
    stepsize_i = learning_rate * rho_i
    step_sizes.append(stepsize_i)

    # 4. update parameters
    grad_parameters = directions * stepsize_i
    m -= grad_parameters[:, 0]
    log_var -= grad_parameters[:, 1]
    
    if iteration % 50 == 0:
        print(f'[iter: {iteration}]\tscore: {scores[best_idx]:.4f}\tstepsize: {rho_i:.3f}')

推論する

推論の際にはすべての弱学習器の結果を足しあわせてパラメータを生成します。

xx = np.linspace(-4, 2, 501) # 予測対象のデータ

xx_mean = np.zeros_like(xx) + np.mean(y)
xx_log_var = np.zeros_like(xx) + np.log(np.var(y) ** .5)

# すべての弱学習器の予測結果を重み付けして足しあわせ
for i in range(len(estimators)):
    xx_mean -= step_sizes[i] * estimators[i][0].predict(xx.reshape(-1, 1))
    xx_log_var -= step_sizes[i] * estimators[i][1].predict(xx.reshape(-1, 1))

# 分散に関してはログ変換しているので元に戻す
xx_var = np.exp(xx_log_var)

推論結果をグラフにすると以下のようになります。 ノイズ項を入力 $X$ の絶対値に比例するように設定していますが、それに追随するような予測値が得られていることがわかります。

f:id:dette:20200501014022p:plain
一次元入力データに対して正規分布の平均・分散(var)の対数を推定した結果。上段の点が学習データ、オレンジのラインが予測平均(点推定値)で塗りつぶされているのがそれぞれの分散。学習データのばらつきが大きい領域ほど推定された分散がおおきいことがわかる。下段は真の分散と予測分散をプロットしたもの。概ね真の値を推論できていることがわかる。

Boosting Iteration の進行と予測の変化をアニメーションにすると以下のような感じになります。初め平均値を合わせたあと、段々と分散が正しい値に落ち着いていっていることがわかります。

f:id:dette:20200501020029g:plain
学習の進行 (Boosting Iteration) と予測値の変化をアニメーションにしたもの。平均値をざっくりとあわせたあとだんだんと分散が正しい値に落ち着いていっていることがわかる。

弱学習器に木構造を使っているので当たり前ですが、外挿される値は平均も分散も定数になるのが面白いなと思いました。Dropout MCMC の結果 と比べると大きく違う傾向ですので、タスクによってどのように差がでるのかなども見てみたいです。

上記コードですが、Kaggle Notebook: https://www.kaggle.com/nyk510/simple-ngboost-implementation にも載せてみました。 ご参考までに。

感想

  • XGBoost のような GBDT を拡張するという話かと思っていたらちょっとちがった
  • 点推定の性能評価で Gradient Boosting 比較しかなかったので GBDT と比較したときどの程度性能が落ちるのか見てみたい。
  • ベースは Gradient Boosting なので GBDT のような葉の重みに対する正則化 (reg_alpha / reg_beta) が入っておらず weak learner が overfit しないようにうまく調整する必要がありそうだと感じた。実装していても max_depth を指定しなかったり min_samples_leaf を大きくすると予測結果がかなりギザギザして木構造っぽい予測になりがちなので、データに合わせたチューニングが必要そう。
  • table1~3にあるNavalってデータセットは簡単すぎんか? (全部予測誤差 0)。逆にどんなのか見てみたい気もする。
  • linear search でステップ幅を探索する部分があるのが気になった。
    • これらは各 boosting で選ばれたデータに対して学習したモデルの方向で探索するためたまたまいい方向が選ばれたり・逆があったりで大きく $\rho$ の値が変わる場合がありそう。
    • これによって boosting 回数の合計を制御しても, colsample / subsample の度合いでステップ合計が変わるので例えば CV で n_estimator を決定して再学習させたとき期待した量だけパラメータが更新されない可能性があるような気がしている。 (GBDTの場合は learning_rate のみなので n_estimators さえ決めれば更新ステップ量は固定できる) そんな変なこと考えなくてもよい?

参考文献

  • NGBoost: https://github.com/stanfordmlgroup/ngboost
    • 公式の実装です。普通に使う分にはこちらを使いましょう。
  • 統計的学習の基礎 (第十章)
    • 通称カステラ本。久しぶりに開きましたが Boosting の章は得にわかりやすくてとても好きです。
  • Gradient Boosting と XGBoost · ZABURO
    • Pure な Boosting から Gradient Boosting ~ XGBoost までの流れを丁寧に説明されているブログです。

*1:ヘシアンが対角行列ですので, gradient と事前に割り算をして計算量を減らせますが定式化通りやっています

python の DateTime・Timezone と Django での取扱い

django でアプリケーションを作っていて timezone 周りで困ったことがあったので、初めて django および python での時間の取り扱いについて真面目に調べてみた記録です。

はじめに・ われわれが使っている時計の時間とは何か

まず時計が示している時間とは何かということを考えてみましょう。 世界中に時計がひとつしかなければ(すべての人間がひとつの時計を正として生活していれば)特に問題はなくその時計の時刻自体を YYYY-m-dd HH:mm などで保存したものが常に正しいので何も考えることはありません。 皆同じものを参照しているので不整合はありません。

しかし実際には世界には時差があります。時差があるというのはどういうことかというと時計を文字通り読み取った値が同じでも、物理的には違う時間なことがあるということです。

例えば東京の 4/23 18:00 とニューヨークの 4/23 5:00 は物理的には同じ時間を表していますが、時計は違う時刻を指しています。それは東京・ニューヨークが別の Timezone に属していているためです。

Timezone ごとに、標準時 (イギリスのあたり) の時計に対して何時間ずらすのかが決まっています。例えば東京では標準時から +9H ・ ニューヨークなら -4H ときまっているので、先の 18:00 と 5:00 が同じ物理時間を指していることが確認できます (差分が 9 + 4 = 13 時間であるからです)。

python で時刻を扱う

これらの問題は python で時刻を取り扱う場合でも同様なので、 単に数値としての時刻と Timezone の情報を持った時刻を区別して取り扱うように実装されています。 前者を naive time, 後者を aware time と呼びます。

naive time は単に datetime として作成したときの状態です。標準ライブラリの datetime で現在時刻をとったような状態では naive time です。

from datetime import datetime

naive = datetime.now()  
# datetime.datetime(2020, 4, 23, 10, 15, 27, 886221)

aware な時刻かどうかは datetime object が持っている tz によって判定されています。naive な datetime は tz を持っていません。

naive.tzname() is None  # True

tzを設定するひとつの方法が naive.astimezone です。この関数は tzinfo instance を引数にとってその timezone の aware な時刻のオブジェクトを返します。 tzinfo の設定は pytz を使うのが楽ちんです。例えば UTC 時刻として登録する場合は以下のようになります。

import pytz
naive.astimezone(pytz.UTC) 
# datetime.datetime(2020, 4, 23, 10, 15, 27, 886221, tzinfo=<UTC>)

一度 aware になってしまうと naive な時刻とは引き算をすることができません。

くどいですが区別するのは Timezone 込の情報と Timezone 無しのものは比較することができないからです。 これら2つを同様に扱ってしまってナイーブに引き算などすると例えば東京とニューヨークを同様に扱って引き算できてしまうので、先の例だと時間差がゼロになってしまいます。

実際 naive time と aware time を引き算すると TypeError が発生します。親切ですね。

naive - naive.astimezone(pytz.UTC)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-24-98bec8f1194b> in <module>
----> 1 naive - naive.astimezone(pytz.UTC)

TypeError: can't subtract offset-naive and offset-aware datetimes

astimezone で変換さた時刻は、元の naive な時刻を標準時で持つような tzinfo の時刻となります。 そのためここで tokyo や newyork を設定しても引き算すると時刻の差はゼロです。

tokyo_tz = pytz.timezone('Asia/Tokyo')
nyc_tz = pytz.timezone('America/New_York')

naive.astimezone(tokyo_tz) - naive.astimezone(nyc_tz)
# datetime.timedelta(0)

もし「東京の時計で naive な時刻」を持つような時間を作りたい場合には timezone.locale を使います。 例えば東京とニューヨークそれぞれで同じ時計が表示されている状態の差分を見る場合は以下のようになります。期待通り時差 +9H から -4H を引くので 13時間のズレになっていることがわかります。

# 東京とニューヨークで文字盤上同じ時刻を引き算すると時差と一致する
diff = tokyo_tz.localize(naive) - nyc_tz.localize(naive)
diff.total_seconds() / 60 ** 2  # -13.0

Django での時間の取り扱い

Django では USE_TZ=True のとき時差を考慮した処理を行なうようになります。すなわち djangopython object として扱う時刻はすべて aware な datetime になるということです。反対に USE_TZ=False にしていると tz を設定しない naive な datetime を使います。この場合現在時刻を知りたければ datetime を使うことができます。

from datetime import datetime

now = datetime.now()

一方で USE_TZ=True の場合は aware な方法で now を作成しなくてはいけません。このような場合のために django では django.utils.timezone というモジュールがありその内部に now が実装されています。

from django.utils import timezone

aware_now = tiemzone.now()

この内部では datetime.now を呼び出したあとに UTC を tz に設定しています。 このように django では時差考慮するとなった時には基本的に UTC として取り扱います。これによってデータベースに保存する際にすべての時刻が UTC であることが担保されますので、引き算等の演算での不整合が発生することを防ぐことができます。

# django/utils/timezone.py

# UTC time zone as a tzinfo instance.
utc = pytz.utc

# 中略

def now():
    """
    Return an aware or naive datetime.datetime, depending on settings.USE_TZ.
    """
    if settings.USE_TZ:
        # timeit shows that datetime.now(tz=utc) is 24% slower
        return datetime.utcnow().replace(tzinfo=utc)
    else:
        return datetime.now()

webapplication ではリクエスト元のユーザーが違う timezone に属していることがあります。このような場合にはリクエストごとに Timezone を変更したくなりますが、この用途のために django には activate という関数があり、これを呼び出すとデフォルトの TimeZone を上書きすることができます。

from django.utils.timezone import activate
activate('Your-Favorite-Timezone')

テンプレートエンジンや DRF のシリアライザなどユーザーに情報を返す際にはここで設定された TimeZone での locale datetime への変換がされるので、ユーザーの TimeZone の時計で見た時の時刻をユーザーごとに設定して返すことができます。

Example: time から datetime への変換

「TimeField を今日の時間とみなして、その時間と DateTimeField の値を比較したい」といいう要望があったとしましょう。 Time を date に設定する関数 combine を使うと以下のように書くことができるように思えますが USE_TZ の場合エラーになります。それは to_today_datetime が返している datetime が naive だからです。

from django.utils import timezone

def to_today_datetime(time, now):
    today = now.today()
    dt = timezone.datetime.combine(today, time)
    return dt

time = # TimeField の値. Time 型

now = timezone.now()

# Log という model があったとする
created_at = Log.objects.first().created_at # created_at は DateTimeField
created_at - to_today_datetime(time, now) # TypeError: can't subtract offset-naive and offset-aware datetimes

では utc にすれば OK かというとそうでもありません。 なぜなら今日と認識しているのはリクエストを送ったユーザーですから、ユーザーにとっての時間として表現する必要があるからです。 今設定されている timezone は get_current_timezone で取得することができます。 得られた tzinfo を使って localize することで時間をユーザーのtimezoneでtimeを時計で見たときの時刻に変更することができます。

def to_today_datetime(time, now):
    today = now.today()
    dt = timezone.datetime.combine(today, time)

    # これでは UTC の人にとって `dt` が時計に表示されているときの時刻になる。
    # dt = dt.astimezone(timezone.utc)
    
    # 今のリクエスト context での timezone を取得して、その timezone で `dt` が表示された時の時間に直す
    tz = timezone.get_current_timezone()
    dt = tz.localize(dt)
    return dt

参考

python で logging を止める

はじめに: 基本的なお作法

python の logging の話です。logging そのまま呼び出しもできるのですが若干やんちゃやで、ということが公式ドキュメントに書いています。

ロガーに名前をつけるときの良い習慣は、ロギングを使う各モジュールに、以下のように名付けられた、モジュールレベルロガーを使うことです: https://docs.python.org/ja/3/howto/logging.html#advanced-logging-tutorial

おとなしくガイドに則って, お行儀よく getLogger で名前を指定して logger を取得するようにしましょう。

from logging import getLogger

logger = getLogger(__name__)

# この段階ではコンソールには何も出ない
logger.info('foo')

Handler

  • logger を単に call しても何もおこならないのは logger 単体はログの管理をしているだけで、実際にファイルやコンソールに書き出す仕事は Handler インスタンスが行なっているからです。
  • Handler はいろいろと種類があるのですが普通使うのは StreamHandler(コンソール書き出し) と FileHandler(外部ファイルへの書き出し) が多いと思います。
  • この handler を logger に追加することで、logger → handler に情報伝達されます。たとえば StreamHandler であればログがコンソールに表示されるようになります。
from logging import StreamHandler

handler = StreamHandler()
logger.addHandler(handler)

Level

logging にはその重要度を定める level という概念があります。 logger と handler にはそれぞれ level を設定することができます。 logger は自分に設定されている level よりも高い log が call された時に handler を実行します。

handler にも level があり、呼びだされた level が自分に設定されている level よりも高いときに handler の method が呼び出されます。

# このハンドラが感知するのは INFO よりもレベルが高いもの (i.e. info / warning / error / critical) になる。 
handler.setLevel('INFO')

logger.info('foo') # この段階では logger のレベルが設定されていないので handler も動かない

# logger に level を設定すると動くようになる
logger.setLevel('INFO')
logger.info('foo') # foo

Logging を一時的に止めたい時

本題です。 logger.disabled = True にすると handler の呼び出しが止まります。一時的にログを止めたい時とか便利です。 level=NOTSET でも同じことになりますが, この場合前に設定していたレベルをどこかに覚えておかないと同じ状態に戻せないため disabled を使ったほうが便利です。

logger.disabled = True
logger.info('hoge') # なにもでない

logger.disabled = False
logger.info('piyo') # piyo