scikit-learnの機能を拡張/変更したscikit-learn準拠モデルを作る

これはなにか

scikit-learn APIにはない新たな予測モデルを作成したい。
そのときscikit-learnをベースに作成することで、GridSearchやCrossValidationなどscikit-learnで使える関数をそのまま流用したscikit-learn準拠のモデルを作成できる。

言い方を変えると、完全に0からモデルを作る場合、scikit-learnにあるGridSearchなどに類するものは自分で1から作らないといけない。

scikit-learn準拠の自作予測モデルを作成する

scikit-learn.org

以下の流れでscikit-learn準拠の自作予測モデルを作成する。

1. sklearn.base.BaseEstimatorを継承
2. 回帰の場合sklearn.base.RegressorMixin、分類ならsklearn.base.ClassifierMixinを継承
3. モデルのアルゴリズムをfit関数として定義
4. モデルの評価をおこなうpredict関数を定義

3, 4のfit, predictはsklearnモデルの基礎となるため必須かつ、この関数名じゃないといけない。例えば、GridSearchは内部的には渡されたパラメータをもとにfit,predictを用いているので関数名を変えると利用することができなくなる。

また、fitではモデルアルゴリズムを実装しその結果を返すコードを記述し、predictではfitの返り値を使った計算をおこなって予測結果を出すコードを記述する。

実例

試しに、リッジ回帰を実装する。
classを定義するとき、コンストラクタとして__init__は必須となる。また、ここではclassの引数となる値の代入をするが、ここで計算をおこなう*1と一部APIがエラーになるらしい。

なお、コードは以下の記事をおおいに参考にした。

yamaguchiyuto.hatenablog.com

classを定義するときに、ちゃんと挙動するかを確認する関数は以下の記事が詳しい。例えば、ちゃんとscikit-learn準拠になっているか確認するcheck_estimatorや、 predictをする前にfitがちゃんと機能しているか確認するcheck_is_fittedなど。
この関数を用いると、エラーメッセージとして何が間違っているか明示してくれるのでできる限り使用した方が良い*2。

qiita.com

また、以下の記事のように引数の型定義や長さなどをチェックするassertなども仕込んでいるとなおよい。

danielhnyk.cz

import numpy as np
from sklearn.base import BaseEstimator, RegressorMixin

class RidgeRegression(BaseEstimator, RegressorMixin):
    # コンストラクタの設定
　# 初期値は必須 
    # 今回の場合、罰則項λを渡してモデルオブジェクトを作る。初期値は1.0
    def __init__(self,lamb=1.0):
        self.lamb = lamb
    
    # リッジ回帰のアルゴリズムを実装
    # データはここで渡したり、加工をおこなうこと
    # 具体的には係数を計算
    def fit(self, X, y):
        A = np.dot(X.T, X) + self.lamb * np.identity(X.shape[1])
        b = np.dot(X.T, y)
        
        # fit内で定義された変数にはサフィックスで '_' を用いる(慣習なので必須ではない) 
        self.coef_ = np.linalg.solve(A,b)
        
        # fitではselfを返す
        return self
    
    # ここで推定をおこなう。
   # 実装上の問題で、yを渡す必要があるが必要ない場合はNoneとしておく
    def predict(self ,X, y=None):):
        # fitで計算した係数を用いて推定した値を返す
        return np.dot(X, self.coef_)

このモデルを実際に使うと以下のようになる*3。

#ボストン住宅価格データセットの読み込み
from sklearn.datasets import load_boston
boston = load_boston()

#説明変数
X_array = boston.data

#目的変数
y_array = boston.target

# 罰則項を0.5にする
model = RidgeRegression(0.5)

# fitでモデルの学習
model.fit(X_array, y_array)

# predictでyの予測
y_pred = model.predict(X_array)

また、scikit-learn準拠モデルなのではじめに説明したように既存のscikit-learn APIを用いることができる。

以下はGridSearchCVでパラメータを求める例。

scikit-learn.org

なお、ここではscikit-learn準拠モデルで継承したMixin内で定義されているscoreを用いている。Mixin内のscoreから評価の定義を変えたい場合は自作scikit-learn準拠モデル内でscoreをオーバーライドして定義し直す必要がある。

github.com

nykergoto.hatenablog.jp

from sklearn.model_selection import GridSearchCV

# パラメータの探索範囲の指定
parameters = {'lamb':np.exp([i for i in range(-30,1)])}

model = GridSearchCV(RidgeRegression(), parameters, cv=5)
model.fit(X_array, y_array)

best = model.best_estimator_
# => RidgeRegression(lamb=1.0)

なお、RidgeRegressionを定義するときに引数lambの初期値を設定してないとGridSearchCVなどでエラーになる。

実際にリッジ回帰をおこなう場合のsklearn.linear_model.LinearRegression

実際の sklearn.linear_model.LinearRegression コードは以下となる。

github.com

当たり前だが実際のコードを読むと他にも関数を定義しており、同様に関数を定義すると自作scikit-learn準拠classに所属する自作関数を作ることもできる。また、これらコードを読むことで実装がどうなるか確認することもできる。

自作のTransformerを使う

機械学習モデルを作成する際に、多くのモデルで前処理として標準化や正規化をはじめとした前処理をする必要がある。
scikit-learnではそのためのclassとして、sklearn .preprocessingのStandardScalerやRobustScaler、MinMaxScalerなどのTransformer(変換器)が提供されており、fitやtransform、fit_transform関数で処理をおこなう。

https://helve-python.hatenablog.jp/entry/scikitlearn-scale-conversionhelve-python.hatenablog.jp

mathwords.net

このとき、StandardScalerのような前処理を自作でおこないたい場合、TransformerMixinを継承したclassを作成する。
前節のsklearn.base.RegressorMixinやsklearn.base.ClassifierMixinをベースにしたscikit-learn準拠モデルでは、GridSearchなどが使える機械学習モデルを作成できるというメリットがあった。
本節のTransformerMixinを使うメリットはfit_ transformや、一連の前処理を一括しておこなうsklearn.pipeline (後に記事を書く予定) で自作のTransformerを用いることができる。

scikit-learn.org

zerofromlight.com

qiita.com

以下の流れでscikit-learn準拠の自作のTransformerを作成する。

1. sklearn.base.BaseEstimatorを継承
2. データの情報(統計情報など)を取得する fit関数を定義(情報がいらない場合は特になにもしない)
3. データを加工するtransformを定義

実例

外れ値の置き換え

以下では、特徴量の外れ値を処理するためにnumpy.clipを用いて、0.1パーセンタイル点から99.9パーセンタイル点の範囲に収まるようにそれぞれを超える値を置き換えるTransformerを作成する。

class FeatureClipper(BaseEstimator, TransformerMixin):
    def __init__(self, cols_to_clip_lower=None, cols_to_clip_upper=None):
        self.cols_to_clip_lower = cols_to_clip_lower
        self.cols_to_clip_upper = cols_to_clip_upper
    
    def fit(self, X, y=None):
        # 各特徴量の0.001, 0.999パーセンタイル点を取得
        self.lower_bounds = {c: X[c].quantile(0.001) for c in self.cols_to_clip_lower}
        self.upper_bounds = {c: X[c].quantile(0.999) for c in self.cols_to_clip_upper}
        return self
    
    def transform(self, X):
        # 直接の書き換えが起きないようにcopy
        _X = X.copy()
        
        # 各特徴量を0.999パーセンタイル点に収める(0.999を超える値は0.999で置き換え)
        if self.cols_to_clip_lower is not None:
            for c in self.cols_to_clip_lower:
                _X[c] = _X[c].clip(lower=self.lower_bounds[c])
                
        # 各特徴量を0.001パーセンタイル点に収める(0.001より小さい値は0.001で置き換え)
        if self.cols_to_clip_upper is not None:
            for c in self.cols_to_clip_upper:
                _X[c] = _X[c].clip(upper=self.upper_bounds[c]) 
        
        return _X

先程用いたボストン住宅価格に外れ値を作って適用してみる。

#ボストン住宅価格データセットの読み込み
from sklearn.datasets import load_boston
boston = load_boston()

#説明変数

X = pd.DataFrame(boston.data, columns=boston.feature_names)

# CRIM列をテキトーに外れ値に置き換える
X2 = X.copy()
X2.iloc[0,0] = -500
X2.iloc[1,0] = 500

# 最小と置き換わる値
# X2['CRIM'].quantile(0)
# => -500.0
# X2['CRIM'].quantile(0.001)
# => -247.4954247

# 最大と置き換わる値
# X2['CRIM'].quantile(0.999)
# => 292.4329810000252
# X2['CRIM'].quantile(1)
# => 500.0

# Transformerの適用
tranceformer = FeatureClipper(cols_to_clip_lower=['CRIM'], cols_to_clip_upper=['CRIM'])

X_clipped = tranceformer.fit_transform(X2)

実際元のX2['CRIM']とX_clipped['CRIM']の箱ひげ図を見ると置き換わっていることがわかる。

特徴量追加

下記コードでは特徴量を追加するTransformerを作成する。簡易化のため、ボストン住宅価格のCRIM, ZNを2倍するというよくわからん特徴量を追加する。

from sklearn.base import BaseEstimator, TransformerMixin

class FeatureAdder(BaseEstimator, TransformerMixin):
    def __init__(self, cols=['CRIM', 'ZN']):
        self.cols = cols
    
    def fit(self, X, y=None):
        #　統計情報などは使わないのでそのまま
        return self
    
    def transform(self, X):
         # 直接の書き換えが起きないようにcopy
        _X = X.copy()
        
        # 新たな特徴量の作成
        _X['CRIM_by_2'] = _X['CRIM'] * 2
        _X['ZN_by_2'] = _X['ZN'] * 2
        
        # 新しく変数が出来るので、後で取り出せるよう変数名を格納しておく
        self.colnames = _X.columns.tolist()
        
        return _X


tranceformer = FeatureAdder()

tranceformer.transform(X_df) # fitは使わないのでtransform

その他

なお、Mixin系は以下(Mixin( で検索するとわかりやすい

https://github.com/scikit-learn/scikit-learn/blob/14031f6/sklearn/base.py:embed:cite:w600

参考

コードの一部は以下の書籍にあるコードを参考にした

機械学習を解釈する技術〜予測力と説明力を両立する実践テクニック

*1:下記コードだとself.lamb = lamb*2.0 とか

*2:今回はコードをシンプルにしたいので使用してない

*3:本当はデータ分割をした方がいいが例なのでテキトー

「ドメイン知識」という言葉の解像度を上げてインプットに活かす

この記事はなにか

要約

• 「ドメイン知識」っていうけど、「解の質」を上げるための領域と「イシューの質」を上げるための領域があるのでは
• 「ドメイン知識を得よう！」と思っても領域によってインプットの仕方が異なる
• DSでもいわゆるコンサルタント的な部分でのドメイン知識が今後は大切
• ドメイン知識をインプットするときにどの領域なのか意識するとよい

背景

最近役職や会社のフェイズが変わったこともあり、今まで以上に自分で課題を見つけて自分で仕事を見つけてくることが求められるようになった。

これ系は、DSにおいていわゆる「DSはドメイン知識が大事！」みたいな話に繋がるところではあり、じゃあその「ドメイン知識」ってどう得たらいいのかなーと思い、他の人はどう考えているのだろうと調べたが、結構ふわっとした書き方や「それだけじゃなくね？🤔」と思ったことが多かったので書いた。

また、なんとなく「ドメイン知識」っぽいことをインプットしているものの、目的が曖昧で抽象度が高いのでドメイン知識」の解像度を上げる必要性も感じた。

ドメイン知識の要素分解

ドメイン知識という言葉は前述のように、結構大きな枠組みとして語られている気がする。そのため、何のためのドメイン知識かということをはっきりさせるために、典型的な分析案件のフローをもとに分解すると以下のようになる。

問題設定のためのドメイン知識

クライアントのビジネス目的を聞き出し、潜在ニーズを考え、クライアントが本当に求めていることは何かを知るために必要な、ステークホルダー周辺の知識。

具体的には、その業界の構造やビジネスフロー、業界で使われる用語、一般的に使われている分析内容や実施されていること、（サービスであれば）そのサービスはどういう内容なのか、など。

これは、クライアントの本当に解決したいことをヒアリングからあぶり出すための会話に必要。 つまり、クライアントとクライアントの土壌でビジネスのやり取りや議論をするための前提知識となる。
これがないと、クライアントとやり取りをしても議論をすることができないので、クライアントの言う通りのことをするだけになるが、クライアントが認識している課題と本当の課題(潜在的ニーズ）にずれがあったり、より優先順位が高いアプローチがあることが往々にしてある。

また、これはクライアントと会話をするためだけではなく、課題に対して仮説をたてアプローチを考える際にも必要となる。
仮説を立てる際は、一般的にはビジネスの背後でどういうことが起きていそうなのか、整理し、考え、仮説を立てることが必要になり、このビジネスの背後でどういうことが起きていそうなのかがまさにドメイン知識を指す。

このドメイン知識を得るためには、

• その作業/サービスを実際におこなっている(使っている）人にヒアリング
• その作業/サービスを実際に使ってみる
• その業界のビジネスに関する書籍やレポートを読む

などが考えられる。と、いっても何もわからない状態では何から着手すればいいかわからない場合も多いと思うので個人的には自社の既にドメイン知識がある人に色々と聞いたり、おすすめの書籍やサイト、資料を教えてもらうのが手っ取り早い。

「自社の既にドメイン知識がある人」は、例えば同業界の分析をしたことがあるデータサイエンティストが職種が同じなのでベストですが、知識そのものはコンサル部門や営業部門の人が最も詳しいことが多いので積極的に彼らに助力を得るのが近道となります。

なお、過去にドメイン知識について少し書いたときも同様のことを書いています。

knknkn.hatenablog.com

分析のためのドメイン知識

主に

• 分析用データはどういう癖があるか
• モデルを作成するために、どういった背景がありそうか

の2点のために必要な知識。

前者は例えば、データに欠損がある場合はどういう事情で欠損しているか、どうやって取得されたデータか、カテゴリデータのラベルの意味はどうなのか、など使用するデータを理解し、場合によっては適切に補正するために必要となる。

このドメイン知識を得るためには、基本的にはデータの作成者にヒアリングしつつ、作成者がわからないところだけそれを知っていそうな人に聞く*1くらいだと思います。
つまり、事前にインプットは難しい。なお、このあたりを知るための取っ掛かりを得るのがEDAの役割のひとつ。

後者は、そもそも「モデリング」というものは現実で起きている事象をモデルに落とし込むことを指すので、良いモデルを作るためには現実をよく知る必要がある。そのために、例えば効果を測るにはコントロール変数としてどういったものを入れると良さそうか、交互作用はありそうか、など。分類予測であればラベルA,Bで明確に違いが起きそうな要因（特徴量）とはなにか、、、といったことを考えるための「事象の背景」に対する知識となる。

このドメイン知識を得るためには、「問題設定のためのドメイン知識」と重複しているが、

• その作業/サービスを実際におこなっている(使っている）人にヒアリング
• その作業/サービスを実際に使ってみる

の2点が主。

理由としては、「問題設定のためのドメイン知識」では仮説を考える際に「現実でどういうことが起きていそうか」を知る必要があり、「モデリング」では現実で起きていることを高解像度で知り、それをモデルに落とし込む必要があるので共通して「現実を知る」必要があります。よって、現実を知るために「その作業/サービス」を使うか使っている人に話を聞くことが必要となります。

活用のためのドメイン知識

主に、

• 分析結果をどう実行可能な施策に落とすか
• 施策をどうシステムに落とすか

ための知識となる。

例えば、どれだけ優れた分析をおこなっても、実行不可能な施策であれば誰も使えない。また、実行自体はできても現場の人が使いやすいシステムであったり、腹落ちをしていなかったらそれは使われない施策となってしまう。

そのために、分析結果を受けて実施する人が使いやすい施策およびシステムはどのような内容なのか？を考えるための知識が必要となる。

このドメイン知識を得るためには、

• 施策実施者(≠対面で分析プロジェクトとしてやりとりしているクライアント）の現場はどうなっているか
• 施策実施者の普段の業務はどうなっているか
• 施策実施者の「分析」に対しての気持ちはどうか
• （実施先の）システムがどうなっているか

を知る必要がある。

はじめの3つは、施策実施者にヒアリングしてみないとわからない。

例えば、「工場のAの改善のためにXをすればよい！」とわかり、対面で分析プロジェクトとしてやりとりしているクライアントが納得したとしても、それを実施する工場現場で働いている人に聞いてみると「Xをすることはできるけど、その場合作業Bが遅れてします」「Xをするのは面倒で多分忘れてしまう」「Xのためのモニター画面はあるけど、使いづらいシステムなので施策はおこないづらい」「分析の言っている意味がよくわからんからやり方がわからん」のようなことが原因で実質的に分析結果が使われなかったりする。

そのため、「どういう施策内容なら実施できそうか」を考えるためには実施者にヒアリングするしかない。

また、4つ目のシステムに関しても現状どういったシステムになっているかは既存システムがどうなっているか、新たなシステムを入れる場合はそれが可能かといったことはクライアントのシステム担当者の人にヒアリングをするしかない。

余談：クライアントを介さない自己学習が可能なこと

先程まで、ドメイン知識を「問題設定」「分析」「活用」と分けたが、インプットをする際に「自分一人である程度インプットできるもの」「クライアントあるいは実際に使う人へのヒアリングをしないといけないもの」の2種類があることがわかった。

「自分一人である程度インプットできるもの」はどのタイミングであろうと時間が許す限りおこなうことができるが、「クライアントあるいは実際に使う人へのヒアリングをしないといけないもの」はタイミングやそもそもそれが可能なのか、などインプットを突き詰めるのには限界がある。

つまり、自己学習が可能なことは「問題設定領域のドメイン知識」がメインということになる(「分析領域」も一部可能)。そのため、「XX業界の分析をすることになったので事前に色々知っておくぞ！」と思った場合は「問題設定」をするためにインプットをするという意識で情報を拾っていく意識が必要となる。

前提知識

ここまでは、各ドメイン知識領域がどの部分で必要となり、どうやって得ていくかを書いた。以降では前述の「ドメイン知識の要素分解図」をもとに、各ドメイン知識領域のデータサイエンティストとして使っていく際に知識を得ていく際の優先順位の話をおこなっていく。

今回考えるにあたりコンサル系の本やコンサル寄りのDS本を最近改めて読み直したので、本からの知見をもとに「(DSにおける）ドメイン知識」の観点でまとめた内容となっています。
そのため、まずは中心となる書籍で特に参考にさせていただいた箇所をざっくりと書く。

データサイエンティストの要諦

データサイエンティストの要諦

最近出た本。データサイエンティストのスキル定義でよくみる「ビジネス力」「エンジニアリング力」「統計力」のベン図*2でいう「ビジネス力」にフィーチャーした本は色々あるけれど、それ系の中でも特に良かった。

• モデリングしかできない"データサイエンティスト"は今後厳しい
• モデリングに加え、「分析をおこない、その結果からビジネスに意味のあるインサイトを得る」ことが"データサイエンティスト"に求められている
• 「分析をおこない、その結果からビジネスに意味のあるインサイトを得る」はいわゆるコンサルティング力なのでちゃんとコンサルティング力もつけようね

といったことを主題に深堀りがされている（+そのための教育やスキルアップはどうするかとか）。

イシューからはじめよ

有名な問題解決系のコンサル本。つまり、前述の「コンサルティング力」をどう上手くやっていくか、と考えてもよさげ。

イシューからはじめよ ― 知的生産の「シンプルな本質」

バリューのある仕事とは？そのためにどうしたらいいの？という観点でタイトル通り「イシューからはじめる」ことを解いた本。

バリューの本質は下図のように「イシュー度」と「解の質」の2軸に分けられ、右上にいくほど「バリューのある仕事」としている。

「イシュー度」とは「自分の置かれた局面でこの問題に答えを出す必要性の高さ」 = 「問題に対して何を解くか」「課題設定」、
「解の質」とは「そのイシューに対してどこまで明確に答えを出せているかの度合い」 = 「どうやって解いていくか」「アプローチ方法」、
と定義されている。

一般的に、「解の質」がバリューを決めると考えられ、「イシュー度(課題の質)」は重要視されない。しかし、「解の質」が高くてもそもそもの「イシュー度(課題の質)」が良くないとクライアントにとっての価値が良くない。例えば、何かしらの相談に対して、物凄く詳細でわかりやすく調査レポートを出されても、そもそも調査自体が筋違いな方向性であれば何の意味がない。

そのため、仕事の仕方としてはまずは手当り次第思いついた課題を着手するよりも先に、「イシュー度」が高い問題はどれかじっくり検討をする。その後、その「イシュー度」が高い課題に対して、深堀りや分析を繰り返していくことで「解の質」を上げていくのが一番生産性(=成果/投資時間)が高い。

なお、よくやる「根性に逃げたやり方」は、とりあえず思い浮かんだ課題(検討しきってないので「イシュー度」は低いことが多い）に対してまず着手し、深堀りや分析を繰り返していくことで「解の質」が上がる。また、その繰り返しの中で色々なことが見えてくる（または取り組んでいる課題のよくなさに気づき、課題を変える）ので課題が徐々に変質していくことで「イシュー度」も上がっていく。といった進み方となる。いわゆる、走りながら考えるアプローチ。

そのため、右上の象限にたどりつくには、左回りルートは「一心不乱に大量の仕事を根性でおこうなうことで価値のある仕事にする」やり方なので生産性は低く、右回りルートの方が生産性が高い。

データサイエンスで考えると、例えば「施策Aの効果検証をしたい」という問題がある。その際に、「何を効果とするか」「効果検証手法はどうするか」など様々に考えられる。

左回りルートでは、とりあえずよくされているような効果の測り方、手法で分析をおこなう。その際に、よりよいアルゴリズムに改良したり、前処理の工夫などをおこなって「解の質」を上げる。併せて、分析結果のフィードバックを現場などから得つつちょっとずつ「イシュー度」を上げていくといったアプローチが考えられる。

このとき、アルゴリズムの改良であったり前処理をしたりといったことは非常に時間がかかる。また、分析後にフィードバックを得た結果、分析自体のやり直しが起きることもある。つまり、時間の無駄が多い。

一方で、右回りルートでは、分析をする前に「施策Aの効果検証をしたい」という問題について深堀りをおこなっていく。例えば、「何故施策Aの効果検証をしたいか」「効果を測ってその結果をどう使っていくか」「施策Aの効果をどう測るとビジネスに活かせるか」・・・など、問題そのものの設定をする段階で時間を使う。そして、ある程度その中で筋が良さそうなものに着手し分析を進めていき分析の質(「解の質」）を上げていく。

この例からわかるように、「解の質」を上げるのには分析が伴うため「イシュー度」を上げるより時間がかかるため、右回りルートの方が時間がかからず生産性が良い。

「ドメイン知識」がどう役に立つか

各書籍を併せて考えると、

• データサイエンティストにはコンサルティング力が必要
• コンサルティングを上手くやるには、「イシュー度」を上げてから考える

といえる。

各ドメイン領域の性質

ドメイン知識の各領域に対して、それぞれが「イシュー度」と「解の質」どちらに該当するか考えると以下のようになる。

つまり、「バリューのある仕事をする」ために必要なドメイン知識は「問題設定領域」のドメイン知識となる。
一方で、「分析領域」「活用」のドメイン知識を増やしてもそれはあくまで「解の質」を上げるための知識となる。

また、各ドメイン領域が「コンサルティング」と「モデリング」に分けると以下のようになる

つまり、モデリング屋としてやっていくのであれば「分析領域」のドメイン知識を得ればいいし、コンサルティング力を上げたいのであれば「問題設定」「活用」のドメイン知識を得る必要がある。

まとめ

今までの内容をまとめると以下のようになる。

つまり、

• コンサルティング力のための知識で、「イシュー度」を高い問題設定をするために必要な知識を得たい → 「問題設定領域」のドメイン知識を得る
• コンサルティング力のための知識で、問題に対して「解の質」を高めるための知識を得たい → 「活用」のドメイン知識を得る
• モデリング力のための知識で、問題に対して「解の質」を高めるための知識を得たい → 「分析」のドメイン知識を得る

となる。

ドメイン知識をインプットする際に、「どの領域のドメイン知識なのか？」「このドメイン知識は何に活かせるか」を意識しながら学んでいくとそのインプットの使いどころがわかりやすいのではないか。

Appendix

「イシューからはじめよ」のうち、今回取り上げたのは「イシューからはじめよ」のメインメッセージである何故「イシュー度を上げる必要があるか」という部分ですが、書籍内には（今回書いたドメイン知識以外での普遍的な）「イシュー度」をどうやって上げ、それをもとにどうやって「解の質」上げる方法も書かれてます。更には、質の高いイシューを立てたあとそれをもとにどうビジネスを展開していくか（「イシュードリブン」という書かれ方がされてます）という今回挙げた「コンサルティング力」にまつわる話も大いに書かれています。

また、「データサイエンティストの要諦」は本記事のメインとなる「イシュー度」「コンサルティング」を兼ねる「問題設定ドメイン領域」を特に書きましたが、「活用」の部分も「コンサルティング」のひとつとして重要視して深堀りをして、これらのスキルをどうやって習得していくか？ということが書いています。特に、本書のP39 図表2「アナリティクス一連の流れ」は、前述の「問題設定」「モデリング」「活用」を更に詳細に分け、各フェーズが（DS協会のベン図と類似する）「ドメイン知識」「機械学習＆統計知識」「ITツール操作スキル」のどのスキルが必要かを◎◎△✗でわかりやすく記載されているので必見です*3。

前述のように、この「コンサルティング力」は、データサイエンティスト協会のベン図の「ビジネス力」の根底をなす部分なので、データサイエンティストの「ビジネス力」を上げたい方はこの2冊はおすすめです。

データサイエンティストの要諦

イシューからはじめよ ― 知的生産の「シンプルな本質」

おまけ：「DSはドメイン知識が大事！」系で良かった記事

job.newspicks.com

tjo.hatenablog.com

*1:例えば、工場のログ分析でイレギュラーに欠損している日時がある→データの作成者は欠損が多い理由まではわからない→現場に聞くと、その時間は事故があって緊急停止をしていた、みたいな

*2:「データサイエンティスト ベン図」とかでググってください

*3:本当は載せたいけど著作権の関係で載せないです

「機械学習を解釈する技術」のここがすごい

はじめに

弊社の森下が書籍を出版することになりました。

「機械学習を解釈する技術 ~ 予測力と説明力を両立する実践テクニック」

gihyo.jp

本書のレビューに関わらせてもらったのでここが良かったぞ！という部分を書こうかなと思います。

本書の特徴

本書の概要としては、以下のような書籍となっている

・機械学習モデルは予測精度が高い一方で、何故そのような予測になったかがブラックボックス(解釈性が低い）
・しかし昨今は説明責任を果たすためにモデルの振る舞いの把握が求められている
・このような予測精度と解釈性のトレードオフを克服するための手法を紹介

※ 技術書評者のサイトの「この本の概要」をざっくり要約

拙著「機械学習を解釈する技術」が8月4日に技術評論社さんから出版されます！PFI, PD, ICE, SHAPの4つの機械学習の解釈手法を、ライブラリを使わずPythonでゼロから手法を実装することを通じてアルゴリズムを理解し、データ分析を通じて勘所を掴んで頂けるよう工夫しました！https://t.co/3ofl7zopti

— 森下光之助 (@dropout009) 2021年5月29日

機械学習の解釈性については数年前にInterpretable Machine Learningが公開されてからここにあるSHAPなりICEなりの手法をちょくちょく日本でも耳にすることはあったが理屈を知るのにこのInterpretable Machine Learningを読むなり、原論文を読む必要があった。つまり英語。めんどい。

一応ググってみるといくつか日本語ブログ記事が出ていたり、数ヶ月前に上記Interpretable Machine Learningの日本語化がされたりと日本語で解説を読むこともできるようになってきたが、結局数式の理解が不可欠なのであまり数式が得意でない自分は「なるほど？🤔」で止まっていた。

そのような中、もうちょっと平易な日本語解説でわかりやすくまとまった書籍がないかなーと思っていたところそのニーズを満たしてくれる書籍を出してくれた。神やんけ。

で、詳細は後述の「ここがすごいぞ！」にそれぞれ書きますが、「平易な日本語解説」「わかりやすく」という部分がとても丁寧に工夫されているのが本書の一番の特徴かな、と感じました。

なお、章立ては以下のようになっています。

1章:機械学習の解釈性とは
2章:線形回帰モデルを通して「解釈性」を理解する
3章:特徴量の重要度を知る〜Permutation Feature Importance(PFI)〜
4章:特徴量と予測値の関係を知る〜Partial Dependence(PD)〜
5章:インスタンスごとの異質性を捉える〜Individual Conditional Expectation(ICE)〜
6章:予測の理由を考える〜SHapley Additive exPlanations(SHAP)〜
付録A: R による分析例~ tidymodelsとDALEXで機械学習モデルを解釈する~
付録B: 機械学習の解釈手法で線形回帰モデルを解釈する

ここがすごいぞ！① 線形回帰モデルを起点に「解釈」について各手法を説明している

この本の特徴として、まずはいわゆる機械学習モデル*1ではなく、高い解釈性を持ち、皆が馴染み深い線形回帰モデルを通して予測モデルにおける解釈性とはどういうことか？ということを(全体概要である一章を除いた）初っ端の第二章からまず理解させます。

そして、その馴染み深い線形回帰モデルを起点にまずは「線形回帰モデルだとこうすればこういう解釈がわかる」と説明した後に、同じことに対して「機械学習ではこの手法(各章で紹介する手法)を使うことで同じように解釈ができる」と追従する形で各章の機械学習の解釈手法について解説している。

よくある専門書だと、その手法の解説のみに終始する。つまり0ベースで数式をベースに色々と解説がされ、それを土台に話が展開されていくので、ベースの理解(特に数式）をちゃんとしていないとなんとなくわかった気になりがちだが、まずは馴染み深いしわかりやすい「線形回帰モデル」を土台として、そこから各手法の解説に入るのでめっちゃわかりやすい*2。神かな？

ここがすごいぞ！②数式に対応するアルゴリズムをコードでも書いて説明している

論文はもとより、Interpretable Machine Learningなどいわゆる専門書は数式をベースに解説がされていることが多い。つまり、数式が読めることを前提に、その数式をベースに他の手法との差異やデメリットの説明を数学的にすることが多い。

もちろん、数式こそ何かを説明するのに端的で論理的で明示的な言語ではあるのだが、数式に強くない場合「なるほど？🤔」となってしまいがちである。

本書では、数式とともに数式に対応するアルゴリズムを平易なPythonで書くことで数式そのものだけでは理解が難しくても、アルゴリズムさえ見ればその1つ1つの処理の流れが明示されているのでその手法で何をしているかがとてもよくわかる。また、副次的にアルゴリズムから数式の読み方の学習もできる。神。

ここがすごいぞ！③ 実際の値を確認しつつ解説する

その手法の解説の際に、簡単なわかりやすいデータで解説を踏まえた挙動の確認を毎回丁寧にしてくれます。そのため、そのデータの実際の値を見ながらどうなっていったかの確認が非常にわかりやすいです。

イメージとしては、自分で関数とか書いたときにtestとして簡単なデータを渡して処理を確認するのに近い。

また、簡単なデータ例のあとにどの章でも実データのBoston Housing(ボストン住宅価格)を使っての解説もおこなっているので、簡易なデータによる挙動の確認→実データによる解釈の確認といった流れになっていてとても理解がしやすいです。神。

ここがすごいぞ！④ 手法の特徴・メリット・デメリットが腹落ちしやすい

「各手法の特徴はなにか、その特徴を踏まえたメリット・デメリットはなにか」ということを毎回明示的に記載がされてます。また、その際に、他の手法と比較しつつ、今までの解説を踏まえて書かれているのでなんとなくわかった気になるのではなく、ちゃんと「こうなるからこういうデメリットがあるのか！確かに！」と腹落ちしやすい構造となっています。

このあたりは構造として論文の読み方にも近いのでしょうが、この本で書いているような

1. この手法の特徴は？
2. なんでそういうことができるの？→コード書きながら理解
3. 他の手法と比べたときのメリット・デメリットは？ →2をベースにしたコードで再度理解

という流れは、あらゆる手法の理解にも通じるので、つよつよな人はこうやって手法を理解して使いこなしているんだなー、という観点でも非常に勉強になりました。

逆に言えば、こういう流れで理解できるような構成になっているので「つよつよな人の理解の流れで理解できる」からこそ、わかった気になりづらい書籍になっています。神。

まとめ

一般的な、専門書籍や論文を読むのが難しい人にも非常に配慮されている書籍となっていました。また、そのような書籍ではある種「猿でもわかる・・・」的な、ものすごく丁寧になりすぎて逆にわかりづらい・抽象化しすぎて表面しかわからん、みたいなことも起きない良バランスの書籍なように感じました。

自社の人が書いた書籍という点は関係なく、結構まじで自信を持っておすすめできる書籍だと思うので、機械学習に携わる人は是非買ってみてください。

機械学習を解釈する技術〜予測力と説明力を両立する実践テクニック

*1:定義にもよりますが。。。

*2:後述のように数式を省略しているという意味ではない

GBDTのハイパーパラメータの意味を図で理解しつつチューニングを学ぶ

この記事は何か

lightGBMやXGboostといったGBDT(Gradient Boosting Decision Tree)系でのハイパーパラメータを意味ベースで理解する。
その際に図があるとわかりやすいので図示する。

なお、ハイパーパラメータ名はlightGBMの名前で記載する。XGboostとかでも名前の表記ゆれはあるが同じことを指す場合は概念としては同じ。ただし、アルゴリズムの違い（Level-wiseとLeaf-wise）によって重要度は変わるし、片方にのみ存在するハイパーパラメータもあるので注意。

lightgbm.readthedocs.io

また、記事の構成などは以下を大いに参考にさせていただいた。

nykergoto.hatenablog.jp

網羅的には以下の記事もよさげ

qiita.com

そもそもGBDTとは

ハイパーパラメータの話のために、最低限ざっくりGBDTとは。

まず決定木(Decision Tree)を作成し、その決定木で学習しきれなかった部分(誤差)を予測するための新たな決定木を作成し、更にその決定木で学習しきれなかった部分(誤差)を予測するための...ということを繰り返し、それらのシーケンシャルな決定木の結果を統合(Boosting)することで高い精度で予測をおこなう手法です。

copypaste-ds.hatenablog.com

www.acceluniverse.com

qiita.com

つまり、「決定木を」「複数つくり」「結果を統合する」。この3点を意識するとGBDT系のハイパーパラメータの意味は概ね理解できる気がします。

以降では基本的に上図をベースに図示しつつポイントを赤字で記載していきます。

ハイパーパラメータ

ハイパーパラメータは主に

• どのような決定木を作るか
• どういった目的関数にするか
• どうやって学習していくか

の3パターンがあります。

以下では、それぞれ主要なハイパーパラメータを記載します*1。なお、名前はドキュメントのはじめに書いているものを記載しますが、よくみるエイリアスとなるハイパーパラメータ名を()で併記しています*2*3。

どのような決定木を作るか

num_iterations (n_estimators)

作成する木の数(図中のNに相当)。基本的にアーリストッピングを使うと思うので、その場合は実質無限に相当する値を指定しておけば良いっぽい。

default = 100

max_depth

木の深さの最大値。大きすぎるとオーバーフィッティングに繋がるので3~8くらいでやることが多い。7ぐらいが無難っぽい。

default = -1(上限なし)

num_leaves (max_leaves)

葉(ノード)の最大数。大きいほど複雑になるが過学習につながる。複雑さに直結するので最重要かもしれない。

例えば、下図では葉が6つあるのでnum_leavesを5以下にすると木の形が変わる。

1分割される毎に葉は2つ増える*4。つまり、max_depthが大きいほど分岐が深くなり葉は増えがちなので合わせて調整する。また、理屈上2^max_depthより大きくはならないのでそれ以上大きい値にするとエラーになる。
転じて、最大ノード数の7割にしたい場合max_leaves = int(0.7 * 2 ** max_depth)のように指定したりするらしい。

default = 31

min_data_in_leaf(min_child_samples)

葉に所属する(割り振られる）最小データ数。つまり、分岐条件で割り振り数がこれ以下になるような分岐はされなくなる。値が小さいと細かい葉の分割もされるようになるが過学習に繋がるし、逆に大きいと分割が大雑把になる。

そもそも学習に使うデータが少ないと、各葉内のデータ自体が少なくなりがちなので学習データ数によって小さくする必要がある。また、前述のnum_leavesにも影響される。

default = 20

min_gain_to_split (min_split_gain)

決定木ではゲインが最大となるような分割がされる。

その分割するゲインの最小値を指定。それ以下のゲインでしか分割ができない場合は分割がされない。

default = 0.0

min_sum_hessian_in_leaf (min_child_weight)

葉を分割するのに必要な（ロスの）hessianの合計値。小さければ小さいほどロスを小さくしようと葉を分割するが、オーバーフィッティングを引き起こす。

default = 1e-3(=0.001)

正直私もちゃんと理解できてないが、意味の詳細は以下に詳しい（min_sum_hessian_in_leafはXGBoostではmin_child_weight）

stats.stackexchange.com

feature_fraction (colsample_bytree)

各木を作成するときに使用可能な特徴量の割合（何パーセントの特徴量をランダムで利用するか）。ちなみに、ランダムなので各木で選ばれる特徴量は基本的に異なる（たまたま同じなことはある）。

default = 1.0

bagging_fraction(subsample)

使用するデータの割合（何％の訓練データをランダムで利用するか）

subsample_freqと合わせて使う。

default = 1.0

boosting(boosting_type)

boostingのアルゴリズムをどうするか。

基本的にはgbdt(勾配ブースティング)を使う。ただ、dart( Dropouts meet Multiple Additive Regression Trees) や, rf(random foreset )、goss(Gradient-based One-Side Sampling) も指定できる*5。

stacking前提だとgbdt以外のアルゴリズムを使うとアルゴリズムだけ変えたモデルが複数できるので多様性的によさげ。

どういった目的関数にするか

目的関数は 目的関数 = 損失関数 + λ正則化項で表され、この目的関数を最小化(最大化)するようなパラメータを求めるのが機械学習のタスク。以下ではこの各項をどうするかのハイパーパラメータ。

objective

損失関数をなににするか。これはハイパーパラメータではなく、機械学習の目的自体なので調整対象外。

lambda_l1 (reg_alpha)

L1正則化項の影響力を調整する係数λ*6部分。
L1正則化に相当するため、重要じゃない特徴量が落とされる。過学習を抑制する。

default = 0.0のため、デフォルトでは正則化はかからない。

lambda_l2 (reg_lambda)

L2正則化項の影響力を調整する係数λ部分。
L2正則化に相当するため、重要じゃない特徴量の影響が小さくなる。過学習を抑制する。

default = 0.0のため、デフォルトでは正則化はかからない。

なお、lambda_l1とlambda_l2どちらとも指定する場合、いわゆるエラスティックネットの形になる。

どうやって学習していくか

以下はチューニングしても意味がないのでチューニング対象外(参考）。

learning_rate

学習率。各木を足し合わせるときの重み。つまり、大きくするほど誤差に対して各木1つ1つの情報を多く使うのでトータルとして使用する木の数が減る。つまり、学習完了までの時間が短縮できるがその分学習に使用する木の本数が減っているので精度は落ちる。 そのため、時間が無限にあったらものすごく小さくすればいい気もするがいったん0.01あたりに固定するのが無難なようだ。

default = 0.1

eval_metric

validation データを評価する損失関数。デフォルトではobjectiveと同じになる。

early_stopping_round(early_stopping)

validation データに対して損失関数(eval_metric)の改善がされなかったら学習を途中で止める(いわゆるアーリーストッピング)。

default = 0 なので、デフォルトではアーリーストッピングはかからずnum_iterationsの指定数の木で学習する。

verbose

指定した学習回数毎にeval_metricを出力。

seed (random_state)

新たな木を作成する際の初期値として、木の学習に使うインスタンスの選択や特徴量候補の選択はランダムにで選択される(Bagging。なお、それぞれ自体の割合はbagging_fractionとfeature_fractionで制御）。そのときの乱数シード。

固定しておかないと再現性が担保できないので固定しとく。

おまけ：ハイパーパラメータをどういじるか

lightGBMのドキュメントにXをしたい場合はどのハイパーパラメータをどう調整すればいいのか？といったことが書いている。

前述のように、それぞれのハイパーパラメータがなにを調整しているか意味を理解した上で読むとなるほどねーとなると思うのでおまけとして書いておく。ちなみに、図は先程までの図のコピーです。

なお、以下の記事の「5. LightGBMのハイパーパラメータ」に端的に書いている。

www.codexa.net

良い結果のために特に重要

lightgbm.readthedocs.io

以下の3つがTreeの複雑さに関する重要度なハイパーパラメータなので、とりあえずこれらをちゃんと調整すると良さげ。

num_leaves

Treeの複雑さを制御する。理論的には2^max_depthが最大の複雑さだが、オーバーフィッティングに繋がる。

前述のように、最大ノード数の7割にしたい場合max_leaves = int(0.7 * 2 ** max_depth)のように指定したりするらしい。

min_data_in_leaf

num_leavesと同様にTreeの複雑さを制御する。オーバーフィッティングを防ぐのに重要。サンプルサイズとnum_leavesに依存するが大規模データセットだと一般的に数百から数千くらい。

max_depth

木の最大の深さを制御するので、深くなるほど複雑になる。つまり、Treeの複雑さを制御する。

速度を上げる

lightgbm.readthedocs.io

計算リソースを増やす

num_threadsで並列計算をするスレッド数を増やせる。

ツリーのノード数を減らす

LightGBM の総学習時間は、使用するツリーノードの総数に応じて増加するのでツリーのノード数を少なくすると高速化する。具体的にはmax_depth、num_leavesを小さくする。ただし、学習が浅くなるので精度が落ちる可能性があるので注意。

(max_depth=2と4) (max_depth=2と4)

また、min_gain_to_splitを増やす(デフォルトは0.0)ことで、分割をするために必要な最低ゲインを調整できる。例えば、ゲインが0.00001のような分割はほぼ分割されてないのと同じ(情報量が増えない)なのであまり意味のない分割となる。そのため、意味のあまりない分割しかできないくらい分岐がされたらその段階で分割をストップするようにする。

更に、レコードがあまり存在しないようなノードが作成されることがあるがこれはオーバーフィットの要因にもなり捉え方によっては不要な分割となる。前述の'max_depth' 'num_leaves'で結果的に制御されることもあるがmin_data_in_leafと min_sum_hessian_in_leafで明示的に制御することもできる。
具体的には、min_data_in_leafはノードを作るのに最低限必要なレコード数 で、 min_sum_hessian_in_leafはノードを分割するときのロスの(hesianの合計の)小ささ。なのでこれらを大きくすることでノードを作成しづらくできる。

木の数を減らす

作成する木の数(boosting round)を制御するnum_iterationsを小さくする。なお、その際は学習時間には関係ないが精度が落ちるのでlearning_rateを増やす必要がある。

また、アーリストッピング(early_stopping_rounds)を用いることで学習する木の数を小さく(途中で学習をやめる）することができる。

分割数をへらす

木の分割数を減らすことで学習時間を短縮できる。ただし、その分木がシンプルになるので精度が落ちる可能性が高くなる。

lightGBM特有だが、データを読み込んだ際に各特徴量ごとに値のbinningが自動でされる。これはmax_binで binを最大でいくつ作るかを制御できるので値を小さくするほど作るbinの数が減って分割時の計算が簡単になる(max_bin_by_featureで特徴量毎に決めれる）。

また、max_binと似ているがbin内のデータ数の最小値min_data_in_binを大きくすることで結果的にmax_binを小さくするのと同じようなことが起きてbinの粒度を荒くできる（イメージとしては）。

feature_fractionは分割の検討をする特徴量の割合を制御する。デフォルトでは1.0なので全特徴量を検討するが、例えば0.5のときはランダムで半分の特徴量しか分割の検討をしない。つまり、検討数が減るため高速化する。

max_cat_thresholdは特徴量の分割点数を制御する。例えば3であれば分割の仕方が最大で3分割となる。そのため、小さくするほど計算が少なく住むので値を小さくすると高速化する。

データ数を減らす

各木の作成時にデフォルトでは常に全データを使うが、bagging_fractionを用いることで使用するデータの割合(ランダムサンプリング)を制御することができる。なお、セットで使うbagging_freqで何回ごとにランダムサンプリングをするかを制御できる。

精度を上げる

lightgbm.readthedocs.io

以下は力尽きたので詳細は省き、端的に。基本的には複雑な木を作る、という意味合い。そのため、どの場合も学習時間が伸びるとともに、(トレーニングデータへの）過学習に繋がりやすいので注意。

大きなmax_binを使う

max_binを上げるということは 、要は分割点を増やす=複雑な木にする、ということなのでその分精度が上がる。

小さいlearning_rateと大きなnum_iterationsを使う

learning_rateを小さくするほど多くの木を使用することになるので精度を上げることができる。また、この際に作成する木の上限数自体が少ないとあまり意味がないのでnum_iterationsも増やす。

大きなnum_leavesを使う

複雑な木を作成

トレーニングデータを増やす

学習に使うデータは多いほど良いよね的なよくある話。

boostingのアルゴリズムでdartを使う

GBDT(のMART)の正則化としてDrop Outの概念を入れて過学習を防ぐアルゴリズムらしい。

詳細は以下がよくまとまっている。

www.simpletraveler.jp

オーバーフィッティングを避ける

lightgbm.readthedocs.io

こっちも力尽きたので詳細は省き、端的に。こっちはシンプルな木を作る、という意味合いが強い。そのため、どの場合も学習時間が速くなるが、(トレーニングデータへの）精度は落ちやすいので注意。

小さなmax_binを使う

binの数が減るのですなわち分岐点を大雑把にする

小さなnum_leavesを使う

木をシンプルにする

min_data_in_leafやmin_sum_hessian_in_leafを増やす

データの分割を大雑把にする

bagging_fractioとbagging_freqをつかってbaggingの調整をする

木毎に使われるインスタンスが異なりやすくなるため、木の多様性に繋がる。

feature_fractionを使って特徴量のサンプリングを調整する

木毎に使われる特徴量が異なりやすくなるため、木の多様性に繋がる。

トレーニングデータを増やす

学習に使うデータは多いほど良いよね的なよくある話。

lambda_l1、lambda_l2、min_gain_to_splitで正則化の調整をする

正則化項の影響を強める

max_depthを小さくする

木が深くなりすぎないようにすることで複雑さを抑える

extra_treesでExtremely Randomized Treesを作る

理屈上、多様な木ができるため汎用性が上がる

path_smoothを増やしスムージングの調整をする

正則化項の影響を強めるのと同意。

参考

www.codexa.net

neptune.ai

有名ライブラリと比較したLightGBM の現在のP29あたり

*1:一部ハイパーパラメータ以外もあり

*2:エイリアスは1つだけ書いているが、他のエイリアスも存在することもある

*3:あえてエイリアスの方で書いている理由はなんだろ？多分そっちの方が直感的にわかりやすい名前とかなのか？

*4:C4.5とかだと3つ以上のときもあるがlightGBMはCARTっぽいので2つしか増えない

*5:例えば、rfってすると普通のRabdom Forestが作られる気はしますがboosting_typeって名前的にboosting系文脈なので本当に普通のRandomForestが作られるのか少しひっかかる。RFはboostingではないような・・・

*6:エイリアスパラメータ名的にはα表記の方が正しいかもしれない

Optunaを使ってみる

Optunaとは

ざっくり書くと、 良い感じのハイパーパラメーターを見つけてくれる ライブラリ。

ちゃんと書くと、

Optuna はハイパーパラメータの最適化を自動化するためのソフトウェアフレームワークです。ハイパーパラメータの値に関する試行錯誤を自動的に行いながら、優れた性能を発揮するハイパーパラメータの値を自動的に発見します。現在は Python で利用できます。

Optuna は次の試行で試すべきハイパーパラメータの値を決めるために、完了している試行の履歴を用いています。そこまでで完了している試行の履歴に基づき、有望そうな領域を推定し、その領域の値を実際に試すということを繰り返します。そして、新たに得られた結果に基づき、更に有望そうな領域を推定します。具体的には、Tree-structured Parzen Estimator というベイズ最適化アルゴリズムの一種を用いています。

とのこと。

ただ、実際に触ってみると『「ある目的関数に対して最適化をする」に対してハイパーパラメータチューニングで使いやすいようなフレームワークになってる』の方が正確な印象を受けた。つまり、ハイパーパラメータチューニング以外の最適化もできる。

tech.preferred.jp

optuna.readthedocs.io

PyData.Tokyo Meetup #21 講演資料「Optuna ハイパーパラメータ最適化フレームワーク」太田 健 from Preferred Networks

www.slideshare.net

実装1: 簡単な例

まずはハイパーパラメータの最適化ではなく、普通の最適化について例示することで全体でやってることの雰囲気を掴む。

以下のような、 でxについての最適化(yの最小化)を考える。

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-2, 2, 100)
plt.figure(0)
plt.plot(x, (3*x**4 - 2*x**3 - 4*x**2 + 2))
plt.show()

評価関数

まずは、評価関数を定義する。今回はそのままyが評価関数となる。

# 目的関数 (f) 
def f(x):
    return (3*x**4 - 2*x**3 - 4*x**2 + 2)

目的関数

次に、先程の評価関数を用いた目的関数を定義する。

# f(x)をラップするobjective関数を定義
def objective(trial): # 引数 (trial) はTrial型の値
    x = trial.suggest_uniform("x", -5, 5) # 試すxを指定範囲から選ぶ (parameter suggestion) 
    ret = f(x) # 探索 (トライアル) の途中状態を持つ
    return ret

このとき、xの探索範囲も併せて内部で定義をする。上コードではtrial.suggest_uniform("x", -5, 5)が探索範囲に相当する。このtrial.suggest_uniformメソッド部分は探索する内容によってメソッドを変える必要がある。今回は連続値なのでtrial.suggest_uniformだが例えば、対数的な取り方にしたい場合はtrial. suggest_loguniformになるし、カテゴリカル変数であればsuggest_categoricalになる。

optuna.readthedocs.io

最適化

次に先程の目的関数の最適化をおこなう。

study = optuna.create_study(direction="minimize") # 最適化処理を管理するstudyオブジェクト
study.optimize(objective, # 目的関数
               n_trials=100 # トライアル数
              )

まずはStudy型の変数 (study) を生成するためにcreate_study関数でインスタンス初期化をする。

このとき、どのような最適化をするか予め指定をおこなう。

今回は最小化をしたいのでdirectioはminimizeとする。
他にも、storage(これまでのトライアルの結果をどこに保存んするか：デフォルトはInMemoryStorage)、sampler(次のトライアルのパラメータをどう選択するか：デフォルトはTPE)、pruner(トライアルを途中で打ち切るかジャッジ：デフォルトはMedianPruner)といったことを指定することができる。

optuna.readthedocs.io

optuna.readthedocs.io

optuna.readthedocs.io

そして、optimize関数で目的関数と、トライアル数(何回探索するか)n_trialsを指定して最適化をおこなう。

それぞれ実行結果は上キャプチャのように各トライアル毎のxと評価関数の結果、その段階の一番良い結果（評価関数の最小結果）を出力する。

この結果は、studyインスタンスに保存されている。

study.get_trials()で各トライアルの情報がlistで取り出すことができる。また、最小値(目的変数）とそのときのxをそれぞれbest_value、best_paramsで出力できる。

# 探索後の最良値
print(study.best_value)
print(study.best_params)

# 探索の履歴
for trial in study.get_trials():
    print(f"{trial.number}: {trial.value:.3f} ({trial.params['x']})")

実装2: lightGBMでの例

lightGBMのハイパーパラメータチューニングについて、今回はboosting_type num_leaves learning_rateの3つをOptunaで自動チューニングします。コードはPyData.Tokyo Meetup #21 講演資料を参照。

データは以下のように'iris'を使います。

import lightgbm as lgb, numpy as np, optuna, sklearn.datasets, sklearn.metrics
from sklearn.model_selection import train_test_split

iris = datasets.load_iris()

data = iris.data
target = iris.target
train_x, test_x, train_y, test_y = train_test_split(data, target, random_state=0)

通常、lightGBMでは以下のようにハイパーパラメータ(param)を決め打ちだったりGridSearchするなりで決めます。

def main():
    param = {
        'objective': 'multiclass',
        'num_class': 3,
        'boosting_type': 'gbdt',
        'num_leavrs': 3,
        'learning_rate': 0.1
    }
    
    train_xy = lgb.Dataset(train_x, train_y)
    val_xy = lgb.Dataset(test_x, test_y, reference=train_xy)

    gbm = lgb.train(param, train_xy, valid_sets = val_xy)
    
    pred_proba = gbm.predict(test_x)
    pred = np.argmax(pred_proba, axis=1)
    return sklearn.metrics.accuracy_score(test_y, pred)

print('Accuracy:', main())

このコードをもとに、boosting_type num_leaves learning_rateの3つをOptunaで自動チューニングするコードに書き換えます。

def objective(trial):
    param = {
        'objective': 'multiclass',
        'num_class': 3,
        'boosting_type': trial.suggest_categorical('hoge', ['gbdt', 'dart']),
        'num_leavrs': trial.suggest_int('num_leaves', 10, 1000),
        'learning_rate': trial.suggest_loguniform('learning_rate', 1e-8, 1.0)
    }
    
    train_xy = lgb.Dataset(train_x, train_y)
    val_xy = lgb.Dataset(test_x, test_y, reference=train_xy)

    gbm = lgb.train(param, train_xy,valid_sets = val_xy)
    
    pred_proba = gbm.predict(test_x)
    pred = np.argmax(pred_proba, axis=1)
    
    acc = sklearn.metrics.accuracy_score(test_y, pred)
    return acc

study = optuna.create_study(direction='maximize')
study.optimize(objective, n_trials=100)

ほぼ一緒ですが、探索対象であるboosting_type num_leaves learning_rateそれぞれがtrial.suggest_...にして探索範囲を指定しています。また、main()をobjective(trial)のように、引数としてtrialを渡しています（なお、main(trial)のままでも動きはしますが色々ややこしいのでobjective(trial)としています）。

その後、この目的関数に対して、create_studyおよびoptimizeをすることで今回指定した探索対象/範囲でのreturn部分(accuracy_score)の最大となるハイパーパラメータを探索してくれます。

# 探索後の最良値
print(study.best_value)
print(study.best_params)

=>
0.9736842105263158
{'hoge': 'dart', 'num_leaves': 687, 'learning_rate': 0.28862099397009183}

実際の探索範囲などは以下の1つ目のpdfのP31や、2つめの記事の「チューニング」見出しあたりが参考になる。

PyData Tokyo Meetup #21 LightGBM

nykergoto.hatenablog.jp

なお、lightGBMに限ると、optuna側でlightGBMで簡単に利用するためのLightGBM Tunerが用意されている。

これは、通常のlightGBMを使う際に、import lightbgmをfrom optuna.integration import lightgbmに差し替えるだけでそのまま既存のコードで動かすことができる。

ただし、全ハイパーパラメータを調整してくれるわけではなく、以下のみ調整してくれる。

• lambda_l1
• lambda_l2
• num_leaves
• feature_fraction
• bagging_fraction
• bagging_freq
• min_child_samples

そのため、他も調整したい場合はLightGBM Tunerを使わず前述のobjectiveを使う方法でやる必要がある。

optuna.readthedocs.io

# ベースとなるパラメータ
param = {
        'objective': 'multiclass',
        'num_class': 3,
}

train_xy = lgb.Dataset(train_x, train_y)
val_xy = lgb.Dataset(test_x, test_y, reference=train_xy)

gbm = lgb_opt.train(param, train_xy, valid_sets = val_xy
                   )

※'objective': 'multiclass'だと、validationの評価のときにbinary_loglossで評価しようとしてエラー吐くっぽいがバグかやり方が間違ってるかわからん・・・。とりあえず上記だと動かないので注意。

チューニング結果などを取り出したいときは、通常同様にparamsやbest_iteration、best_scoreでアクセスできる。

print(gbm.params)

# =>
{'objective': 'regression', 'metric': 'rmse', 'feature_pre_filter': False, 'lambda_l1': 0.008437118066241178, 'lambda_l2': 6.201313740165131e-06, 'num_leaves': 31, 'feature_fraction': 0.5, 'bagging_fraction': 0.41674201711462905, 'bagging_freq': 2, 'min_child_samples': 20}

tech.preferred.jp

qiita.com

実装3：閾値の最適化

今まではハイパーパラメータを変えて目的関数の最適化だったが、同様のことを分類問題の確率閾値の評価関数最適化にも使ってみる。

まずは、LightGBM Tuner経由で二値分類をおこなう(コードはこの記事参照）。

import optuna.integration.lightgbm as lgb
from sklearn import datasets
from sklearn.model_selection import train_test_split
import numpy as np
import pandas as pd
from sklearn import metrics

# Breast Cancer データセットを読み込む
bc = datasets.load_breast_cancer()
X, y = bc.data, bc.target

# 訓練データとテストデータに分割する
X_train, X_test, y_train, y_test = train_test_split(X, y)

# データセットを生成する
lgb_train = lgb.Dataset(X_train, y_train)
lgb_eval = lgb.Dataset(X_test, y_test, reference=lgb_train)

# LightGBM のハイパーパラメータ
params = {
    # 二値分類問題
    'objective': 'binary',
    # AUC の最大化を目指す
    'metric': 'auc',
    # Fatal の場合出力
    'verbosity': -1,
}

# 上記のパラメータでモデルを学習する
model = lgb.train(params, lgb_train, valid_sets=lgb_eval,
                  verbose_eval=50,  # 50イテレーション毎に学習結果出力
                  num_boost_round=1000,  # 最大イテレーション回数指定
                  early_stopping_rounds=100
                 )

y_prob = model.predict(X_test, num_iteration=model.best_iteration)

このとき、y_predは1となる確率として返ってきます。何も考えずに閾値を0.5として、0.5以上を1、0.5未満を0とするよりこの閾値も調整した方が評価結果がよくなることがあります。

以下の例では評価関数F1の最適化をしています。

# 閾値付きのF1
def f1(y_test, y_prob, threshold):
    y_pred = [1 if prob >= threshold else 0 for prob in y_prob]
    score = f1_score(y_test, y_pred)
    return score

# 目的関数
def objective(trial): 
    threshold = trial.suggest_float('threshold', 0.0, 1.0) # 0~1.0で探索
    ret = f1(y_test, y_prob, threshold)
    return ret

# 最適化
study = optuna.create_study(direction="maximize") 
study.optimize(objective, n_trials=100)


# 探索後の最良値
print(study.best_value)
print(study.best_params)

# =>
0.9710982658959537
{'threshold': 0.019298362201426927}


# 探索の履歴
for trial in study.get_trials():
    print(f"{trial.number}: {trial.value:.3f} ({trial.params['threshold']})")

# =>
0: 0.951 (0.9881809312984907)
1: 0.958 (0.8109994485475401)
2: 0.963 (0.9577886452457406)
3: 0.963 (0.9777420721270973)
4: 0.964 (0.37669967567838136)
5: 0.958 (0.8019522885395326)
6: 0.964 (0.7026656746225185)
7: 0.958 (0.8412677014903346)
8: 0.964 (0.6663596796212387)
9: 0.964 (0.5085085173697037)
10: 0.965 (0.048865011443941064)
11: 0.971 (0.019298362201426927)
12: 0.955 (0.0013199363878491632)
13: 0.966 (0.004203349012100257)
14: 0.970 (0.17029221135052455)
15: 0.970 (0.20333826422743217)
16: 0.970 (0.16478681177393573)
...

以下の記事は評価関数をQWKにして閾値調整をOptunaでおこなっている。

blog.amedama.jp

その他

sample

書き方(このハイパーパラメータのsuggestionは何がいいか？とか)は公式のsampleで使いたいロジックのものを見て真似たら良さげ

github.com

複数アルゴリズムの使用

以下のように、複数アルゴリズム(例ではSVMとRandomForest)で調べることも可能。

def objective(trial):
    classifier_category = trial.suggest_categorical("classifier", ["SVC", "RandomForest"])
    
    if classifier_category == "SVC":
        # ハイパーパラメータ
        svc_c = trial.suggest_loguniform("SVC_C", 0.01, 100) # log一様分布から試すxを選ぶ 
        svc_gamma = trial.suggest_loguniform("SVC_gamma", 0.01, 100) # 同上
        
        # モデル定義
        classifier = svm.SVC(
            C=svc_c, 
            gamma=svc_gamma, 
            random_state=0
        )
        
    elif classifier_category == "RandomForest":
        # ハイパーパラメータ
        randomforest_n_estimators = trial.suggest_int("RandomForest_n_estimators", 1, 3)
        randomforest_max_depth = trial.suggest_int("RandomForest_max_depth", 1, 3)
        
        # モデル定義
        classifier = ensemble.RandomForestClassifier(
            n_estimators=randomforest_n_estimators,
            max_depth=randomforest_max_depth,
            random_state=0
        )
    
    # train
    classifier.fit(tr_x, tr_y)
    
    # pred
    va_pred = classifier.predict(va_x).astype("uint8")
    
    # 評価
    acc = metrics.accuracy_score(va_y, va_pred)
        
    return acc


sampler = optuna.samplers.TPESampler(seed=0)

# Study型の変数 (study) を初期化して生成、sampler(次のトライアルのパラメータを選択)は上記、最大化
study = optuna.create_study(sampler=sampler, direction="maximize")

study.optimize(objective, n_trials=30)

参考

qiita.com

qiita.com

ohke.hateblo.jp

www.nogawanogawa.com

Classを用いて、特徴量作成を仕組み化する@ぐるぐる

atma#10を読んでいて、運営の初心者用講座の「よりたくさんの特徴量の作成」あたりがとてもよかったので自分なりにまとめてみる。

www.guruguru.science

なお、再解釈などはおこなっているが上記リンクをベースに色々と書くだけなので、端的に知りたい場合は上記リンクを読んでください。

これはなにか

特徴量を作成する際に、各処理を1ブロックとして記載することで可読性を上げたり、テストがしやすかったり、train/testで共通のインスタンスを対象とする際にそれぞれで計算をしなくてもすむなど色々と利便性が高かったので使っていきたい。そのため自分用に整理をおこなう。

自作関数での処理との違い

特徴量作成関数の構造 第一回目では関数形式での実装方法をお伝えしました。関数形式でも問題はないのですが、少し不満があります。それは学習時とテスト時で同じ用な計算をしなくてはならない、という点です。

例えば CountEncoding のコードを思い出しましょう。CountEncoding は出現回数の計算を行なう必要がありますが、テストデータやあるいは新しい別のデータに対して CountEncoding する場合には、過去に計算したカウント情報を使って単に変換だけを行なうべきでしょう。(#1でやった方法では学習・推論で二回カウント計算が走っていることに注意してください。) このように機械学習の特徴量変換は

1.学習時に内部状態を記憶しておいて 2.推論時には 1 でつくった情報を用いて変換をする という構造をもっていることが非常に多いです(これは機械学習モデル自体も同様の構造ですね)。ですのでこの論理的な構造をコードにも反映して行きたいと思います。

([講座#2] EDA・モデルの改善から)

上引用の例としてCountEncoding(FrequencyEncoding)が出ています。これは 過去に計算したカウント情報を使って単に変換だけを行なうべきでしょうとあるように「あるカテゴリカル特徴量に対して、TrainとTestで共通の値にする」ということになります。

これは何故かというと、TrainとTestで別の値を持つ場合そのカテゴリカル特徴量がTrainとTestで意味合いが変わる、つまりある意味では別の特徴量なのに同じ特徴料として解釈させて予測させるようなことになるからです。そのため、TrainとTestを合わせた全体でCountをおこなったり、TrainだけでCountした値をそのままTestにも使うのが一般的です。

また、別の理由としてTrain,Test片方にしか存在しないカテゴリがある場合に学習器が機能しないという問題を解決するという意図もあります。
例えばLabelEncoding(OrdinalEncoding)で以下のようなデータを考えます。

Train

動物
犬
猫
犬
犬
猫
犬

Test

動物
猿
豚
猿
豚
犬
犬

このとき、Train/Testがどう変換されるか考えるとそれぞれ以下のようになります

動物	動物(LabelEncoding)
犬	0
猫	1
犬	0
犬	0
猫	1
犬	0

Testでは

動物	動物(LabelEncoding)
猿	0
豚	1
猿	0
豚	1
犬	2
犬	2

みてわかるように、「動物」という特徴量で同じ値、例えば0が指すものがTrainでは犬、Testでは猿となり意味合いが変わってしまいます。また、Trainのみを使ってTestにLabelEncodingしようにもTrainには値「猿」「豚」がないのでエラーとなってしまうのでTrainとTestを合体させたデータでLabelEncodingの変換ロジックを作成する必要があります。

これと同様のことがCountEncodingでもいえます。

話を戻すと、このように全体のデータを使って計算する特徴量は関数で処理している場合、Trainデータに対する処理・Testデータに対する処理
この記事の口座#1のコードではCountEncodingを関数で処理をしているのでコードをみてみると以下のようになってます。

def create_count_encoding_feature(input_df):
    use_columns = [
        'acquisition_method',
        'title',
        'principal_maker',
        # and more
    ]

    out_df = pd.DataFrame()
    for column in use_columns:
        vc = train_df[column].value_counts()
        out_df[column] = input_df[column].map(vc)

    return out_df.add_prefix('CE_')

関数内でtrain_dfのあるカテゴリカル特徴量に対してvalue_counts()をしてCountEncodingで変換したものを引数で渡したDataFrameに渡して上書きしています。

実際にこの関数を使う際は、train_df, test_dfそれぞれにこの関数を走らせます。つまり、vc = train_df[column].value_counts()の部分を計2回走らせることになります。しかし、これは引数にtrain_df, test_dfどちらでもまったく同じことをやるので計算の無駄、1回計算してその値を保持しておけばいいよね、というのが引用の

1.学習時に内部状態を記憶しておいて 2.推論時には 1 でつくった情報を用いて変換をする

で書いていることになります。

参考記事をトレス

ここからは、本論の「Classを用いて特徴量を作成していく」について。
上述のような状態管理のために、Classを用いて特徴量を作成します。この構成としては、AbstractBaseBlockという抽象クラスを作成し、各特徴量作成処理毎にAbstractBaseBlockclassを継承して実装していきます。

各classの中身としては、学習(どう変換などをしていくか）のためのfit関数で状態を更新し保持し、推論（特徴量化)のためのtransform関数でinputしたDFに対して特徴量を作成し、作成された特徴量のDFを返す実装となります。

なお、どの特徴量作成処理classでも共通した名前（中身は異なる）なので、後に説明するように各特徴量作成処理classを一括で処理することができます。

この各特徴量作成処理classをブロックという単位で作成していき、最終的にそのブロックをつなぎ合わせて一括で処理をおこないます。

なお、ここからはしばらく元記事中にあるコードをコピーしつつコメントをしていきます。

ブロックを使った特徴量作成処理(コピペ)

ベースとなるAbstractBaseBlock は以下

class AbstractBaseBlock:
    def fit(self, input_df: pd.DataFrame, y=None):
        return self.transform(input_df)

    def transform(self, input_df: pd.DataFrame) -> pd.DataFrame:
        raise NotImplementedError()

抽象クラスなので、fitとtransformを仮定義。作成したい特徴量によっては内部状態の更新(fit)をしないこともありますが、推論(transform)は必ずおこなうので実装が矯正されます。

内部状態が更新されるブロック例・ CountEncoding

ここでは、CountEncodingのように何かしら状態を保持し、それを用いて特徴量を作成するような処理をブロックで実装します。なお、使用した特徴量の接頭にCE_とつけた新特徴量になる実装となってます。

前述のCountEncodingをおこないます。またこのとき、train,testの全体を使う必要があるので、ついでにtrain,testを結合して読み込むread_whole_df関数も実装します。

まずは前述のようにAbstractBaseBlockclassを継承します。今回はCountEncodingなので、変換したい列名をattributeとして持っておき(__init__)、fitでラベルの出現回数をそれぞれカウントしその値(状態)をattributeに保持し、transformでCountEncodingを適応したいDF(input_df : train_df or test_df)に対して、保持した結果を用いて**CountEncodingをおこないます。

なお、fitでは返り値にtransformの結果が返ってくるのでfit_transformの方がわかりやすい気がしなくもない。

def read_whole_df():
    return pd.concat([
        read_csv('train'), read_csv('test')
    ], ignore_index=True)

class CountEncodingBlock(AbstractBaseBlock):
    """CountEncodingを行なう block"""
    def __init__(self, column: str):
        self.column = column

    def fit(self, input_df, y=None):
#         vc = input_df[self.column].value_counts()
        master_df = read_whole_df()
        vc = master_df[self.column].value_counts()
        self.count_ = vc
        return self.transform(input_df)

    def transform(self, input_df):
        out_df = pd.DataFrame()
        out_df[self.column] = input_df[self.column].map(self.count_)
        return out_df.add_prefix('CE_')

これでCountEncoding処理をおこなうブロックが作成されました。

内部状態更新が行われないブロック例・StringLength

次に、 内部状態の更新がされないブロックを作成していきます。

ここでは入力される元のデータをそのまま加工するような処理をブロックで実装します。 そのまま加工で何かしらの内部状態を更新して保持する必要はないのでfit関数はオーバーロードする必要はありません。

記事中の例では文字数をカウントする処理StringLengthBlockを実装しています。なお、使用した特徴量の接頭にStringLength_とつけた新特徴量になる実装となってます。

class StringLengthBlock(AbstractBaseBlock):
    def __init__(self, column):
        self.column = column

    def transform(self, input_df):
        out_df = pd.DataFrame()
        out_df[self.column] = input_df[self.column].str.len()
        return out_df.add_prefix('StringLength_')

各特徴量処理ブロックをまとめて処理

内部状態を更新するCountEncodingBlock、内部状態を更新しないStringLengthBlockの2ブロックを走らせるためのrun_blocks関数を定義します。

まずはどのブロックをどう使用するかをfeature_blocksとして定義します。

feature_blocks = [
    *[CountEncodingBlock(c) for c in ['art_series_id', 'title', 'description', 'long_title',
       'principal_maker', 'principal_or_first_maker', 'sub_title',
       'copyright_holder', 'more_title', 'acquisition_method',
       'acquisition_date', 'acquisition_credit_line', 'dating_presenting_date',
       'dating_sorting_date', 'dating_period', 'dating_year_early',
       'dating_year_late',]],
    *[StringLengthBlock(c) for c in [
        'title', 'description', 'long_title',
       'principal_maker', 'principal_or_first_maker', 'sub_title',
    ]]
]

ここでは、各ブロックとそのブロックに使用したい特徴量を組み合わせて記載してます。

ちなみに、この[*[counter for counter in iterator]]という表現は

1. リスト内表記
2. *演算子によるunpacking

の組み合わせです。

リスト内表記は [counter for counter in iterator] の形でfor loopが実行されます。コード中のここのみ抜き取ると以下のような結果になります。

また、*演算子はiterator(タプル、リスト、セット)の中身をバラすことができます。つまり、ここでは*なしだと各ブロック数の多次元配列になるので*を用いることで1次元配列にしています。

docs.python.org

qiita.com

そのため、このような書き方をすることで全ブロック処理から生成されたclassオブジェクトを1次元配列に格納しています。

次に、この各classオブジェクトのリストが格納されたfeature_blocksを用いて各classオブジェクト実際に走らせるrun_blocksについて。

def run_blocks(input_df, blocks, y=None, test=False):
    out_df = pd.DataFrame()

    print(decorate('start run blocks...'))

    with Timer(prefix='run test={}'.format(test)):
        for block in feature_blocks:
            with Timer(prefix='\t- {}'.format(str(block))):
                if not test:
                    out_i = block.fit(input_df, y=y)
                else:
                    out_i = block.transform(input_df)

            assert len(input_df) == len(out_i), block
            name = block.__class__.__name__
            out_df = pd.concat([out_df, out_i.add_suffix(f'@{name}')], axis=1)

    return out_df

ここは

1. 各オブジェクト処理を処理時間付きで出力
2. testデータかどうかによって処理を変える
3. おかしなことになってないかassertされる

の3点がおこなわれています。
特徴量処理って結構しくっていることが多いので1,3点目がすごく大事だと思ってます。

2点目のtestはここではtrainデータに対して処理したい場合は各ブロックにある共通関数fitだけをおこない、testデータに対して処理したい場合はtransformのみをおこないます。正確には、fit関数の返り値がtransformで処理した内容なので、trainデータではfitとtransformを、testデータではtransformだけおこなっています。 つまり、trainデータでは各クラスオブジェクトの状態を更新(作成)をする必要があるのでfitを走らせ、testデータは状態の更新をする必要がないのでtransformだけを走らせればよいという仕様になります。 逆に、これがないと冒頭に書いたようにtrainでもtestでも2回同じ計算をすることになって計算の無駄が生じます。

なお、当たり前ですがtestデータの前には状態の更新(作成)の必要があるのでtrainを先に走らせないとエラーになります。

run_blocks(train_df, blocks=feature_blocks)
run_blocks(test_df, blocks=feature_blocks, test=True)

# =>
★★★★★★★★★★★★★★★★★★★★ start run blocks... ★★★★★★★★★★★★★★★★★★★★
    - <__main__.CountEncodingBlock object at 0x7fb981308e80> 0.032[s]
    - <__main__.CountEncodingBlock object at 0x7fb981308978> 0.037[s]
    - <__main__.CountEncodingBlock object at 0x7fb981308f98> 0.035[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132ef28> 0.030[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132e1d0> 0.018[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132e748> 0.020[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132eb70> 0.030[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132e278> 0.019[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132efd0> 0.027[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132e908> 0.015[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132e6a0> 0.019[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132e9b0> 0.017[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132edd8> 0.021[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132e4e0> 0.014[s]
    - <__main__.CountEncodingBlock object at 0x7fb9514c29b0> 0.016[s]
    - <__main__.CountEncodingBlock object at 0x7fb9514c2dd8> 0.014[s]
    - <__main__.CountEncodingBlock object at 0x7fb99184c6a0> 0.013[s]
    - <__main__.StringLengthBlock object at 0x7fb99184c5f8> 0.007[s]
    - <__main__.StringLengthBlock object at 0x7fb9915f4780> 0.007[s]
    - <__main__.StringLengthBlock object at 0x7fb9915f46d8> 0.008[s]
    - <__main__.StringLengthBlock object at 0x7fb9915f44e0> 0.008[s]
    - <__main__.StringLengthBlock object at 0x7fb99180a400> 0.007[s]
    - <__main__.StringLengthBlock object at 0x7fb9812de8d0> 0.006[s]
run test=False 0.502[s]
★★★★★★★★★★★★★★★★★★★★ start run blocks... ★★★★★★★★★★★★★★★★★★★★
    - <__main__.CountEncodingBlock object at 0x7fb981308e80> 0.005[s]
    - <__main__.CountEncodingBlock object at 0x7fb981308978> 0.006[s]
    - <__main__.CountEncodingBlock object at 0x7fb981308f98> 0.005[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132ef28> 0.005[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132e1d0> 0.004[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132e748> 0.004[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132eb70> 0.005[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132e278> 0.005[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132efd0> 0.005[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132e908> 0.004[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132e6a0> 0.004[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132e9b0> 0.004[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132edd8> 0.004[s]
    - <__main__.CountEncodingBlock object at 0x7fb98132e4e0> 0.003[s]
    - <__main__.CountEncodingBlock object at 0x7fb9514c29b0> 0.003[s]
    - <__main__.CountEncodingBlock object at 0x7fb9514c2dd8> 0.003[s]
    - <__main__.CountEncodingBlock object at 0x7fb99184c6a0> 0.003[s]
    - <__main__.StringLengthBlock object at 0x7fb99184c5f8> 0.006[s]
    - <__main__.StringLengthBlock object at 0x7fb9915f4780> 0.007[s]
    - <__main__.StringLengthBlock object at 0x7fb9915f46d8> 0.007[s]
    - <__main__.StringLengthBlock object at 0x7fb9915f44e0> 0.006[s]
    - <__main__.StringLengthBlock object at 0x7fb99180a400> 0.007[s]
    - <__main__.StringLengthBlock object at 0x7fb9812de8d0> 0.007[s]
run test=True 0.183[s]

3点目のassertはここではinputとoutputの行数が一致しない（つまり、処理した結果データ数が増減した）場合assertされるようになります。

で、1点目で処理しているclassオブジェクトが書かれるのでどこでしくっているかとか、処理のボトルネックがどこかがわかる出力となっています。

なお、run_blocksの結果特徴量処理をして作成された新規の列がDFとして返ってきます。

内部状態更新が行われる例2:tf-idf

記事より参照。

tf-idfでテキスト特徴量をベクトル化しSVDで次元圧縮をする。

具体的には、fitでinputされたDFの指定テキスト特徴量に対して前処理(text_normalization関数)およびベクトル化(tf-idf→SVD)をおこなってその状態を記憶する。そして、transformでその記憶された状態を用いて特徴量を作成する。

import texthero as hero
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.decomposition import TruncatedSVD
from sklearn.pipeline import Pipeline

def text_normalization(text):

    # 英語とオランダ語を stopword として指定
    custom_stopwords = nltk.corpus.stopwords.words('dutch') + nltk.corpus.stopwords.words('english')

    # 表記ゆれなどを前処理
    x = hero.clean(text, pipeline=[
        hero.preprocessing.fillna,
        hero.preprocessing.lowercase,
        hero.preprocessing.remove_digits,
        hero.preprocessing.remove_punctuation,
        hero.preprocessing.remove_diacritics,
        lambda x: hero.preprocessing.remove_stopwords(x, stopwords=custom_stopwords)
    ])

    return x

class TfidfBlock(AbstractBaseBlock):
    """tfidf x SVD による圧縮を行なう block"""
    def __init__(self, column: str):
        """
        args:
            column: str
                変換対象のカラム名
        """
        self.column = column

    def preprocess(self, input_df):
        x = text_normalization(input_df[self.column])
        return x

    def get_master(self, input_df):
        """tdidfを計算するための全体集合を返す. 
        デフォルトでは fit でわたされた dataframe を使うが, もっと別のデータを使うのも考えられる."""
        return input_df

    def fit(self, input_df, y=None):
        master_df = self.get_master(input_df)
        text = self.preprocess(input_df)
        self.pileline_ = Pipeline([
            ('tfidf', TfidfVectorizer(max_features=10000)),
            ('svd', TruncatedSVD(n_components=50)),
        ])

        self.pileline_.fit(text) # sklearn.pipelineのfit
        return self.transform(input_df)

    def transform(self, input_df):
        text = self.preprocess(input_df)
        z = self.pileline_.transform(text)

        out_df = pd.DataFrame(z)
        return out_df.add_prefix(f'{self.column}_tfidf_')

内部状態更新が行われる例3:NAを全体の平均値で穴埋め

例用に自作。

NAをを全体の平均値で埋める。

fitで全体の平均値を状態に保存し、tranformではそれを使ってfillnaをおこなう。

class FiillAvgBlock(AbstractBaseBlock):
    """tfidf x SVD による圧縮を行なう block"""
    def __init__(self, column: int):
        """
        args:
            column: str
                変換対象のカラム名
        """
        self.column = column

    def fit(self, input_df, y=None):
        master_df = read_whole_df()
        self._avg = master_df[self.column].mean()
        return self.transform(input_df)

    def transform(self, input_df):
        out_df = pd.DataFrame()
        out_df[self.column] = input_df[self.column].fillna(self._avg, inplace=False)  # 平均値
        return out_df.add_prefix('FillAvg_')

内部状態更新が行われない例1:特徴量の組み合わせ演算

例用に自作。

数値特徴量2つを足し合わせるブロックを作る。

class FeatureAddingBlock(AbstractBaseBlock):
    """数値特徴量2つを足し合わせる block"""
    def __init__(self, column1: int, column2: int):
        self.column1 = column1
        self.column2 = column2

    def transform(self, input_df):
        out_df = pd.DataFrame()
        adding = input_df[self.column1] + input_df[self.column2]
        out_df[self.column1 + '+' + self.column2] = adding
        return out_df.add_prefix('Adding_')

データで例示するのに良い列がないので、特に意味がない特にとなるが、dating_periodとdating_year_lateの和、dating_periodとdating_year_earlyの和の2特徴量を作る。

2列を渡す必要があるのでfeature_blocksではzipを用いて以下のように渡す(ミスが増えそうな渡し方なので他になにかいい渡し方ありそう。いっそ1listの全順列組み合わせをitertools.permutationsで作ったがいいかも？)。

feature_blocks = [
    *[FeatureAddingBlock(i, t) for i, t in zip(['dating_period', 'dating_period'],['dating_year_late', 'dating_year_early'])]
]

run_blocks(train_df, blocks=feature_blocks)

内部状態更新が行われない例2:（log系データに対して）前のレコードとの時系列差分

例用に自作。

（log系データに対して）前のレコードとの時系列差分を取る。

class DiffLagBlock(AbstractBaseBlock):
    """principal_maker視点での前のレコードとの時系列(dating_sorting_date)差分 block"""
    def __init__(self, column: int):
        self.column = column

    def transform(self, input_df):
        out_df = pd.DataFrame()
        copied_input_df = train_df.copy()
        sorted_grouped_input_data = train_df.sort_values(['principal_maker','dating_sorting_date']).groupby('principal_maker')
        lag = sorted_grouped_input_data.shift(1)
        lag = lag.rename(columns=lambda x: 'lag_' + x)
        joined_df = copied_input_df.join(lag['lag_' + self.column])# lagは時系列ソートされてるのでindex結合
        out_df[self.column] = joined_df[self.column] - joined_df['lag_' + self.column]

        return out_df.add_prefix('DiffLag_')

feature_blocks = [
    *[DiffLagBlock(c) for c in ['dating_year_late']]
]
                   
run_blocks(train_df, blocks=feature_blocks)

その他所感

はじめに書いたように、特徴量処理を各ブロックに分けるので可読性も高いし、色々と特徴量用ブロックの引数をいじるなどの試行錯誤がしやすい。また、assertも仕込みやすそう。そのため、この形式で書けそうなら極力書くようにしたい。

一方で、前述の特徴量の和を作るような引数を2つ渡して何かするみたいな処理だと引数の渡し方を工夫しないと可読性が落ちたりとか、ブロックで作った特徴量に更になにかするブロックを作りたかったら、前ブロックで特徴量名が自動生成される関係で次ブロックの引数の指定がちゃんと把握しておかないとやりづらそう（一括でいいならdf.columns.str.contains('XX_')とか使ってするとか？）だったり、処理によってはやりづらそうな場合もありそう。

その場合は無理にこのフレームでやるのではなく、それぞれ別立てで作ればいいのかなーと思う。

自然言語処理を色々楽にするTextheroを使ってみる

Textheroでできること

PythonライブラリTextheroでは、自然言語処理を簡単にできる。機能としては下記が可能。

• 前処理・・・空白やURLの削除、大文字変換など
• 解析・・・文字列が人名か地名かなど(まだあまり機能がない）
• ベクトル変換・・・TF-IDFやPCAなど
• 可視化・・・ワードクラウドや散布図など

なお、まだβ版なようでドキュメントが一部未整備だったり日本語対応をしていない模様。

github.com

DataFrameを用いるサンプルコードのデータはTextheroのサンプル内にあるBBCのレポートを使う。

import texthero as hero
import pandas as pd

df = pd.read_csv(
   "https://github.com/jbesomi/texthero/raw/master/dataset/bbcsport.csv"
)

前処理

自然言語に対して諸々おこなうにあたって、タグの削除や名寄せのための処理（空白の削除や大文字小文字の統一など）を行う必要がある。このライブラリを使わない場合は正規表現でがんばったりする必要があるが（英語であれば）ある程度の処理用メソッドが用意されている。

texthero.org

ただし、現状では日本語未対応なので日本語に対してこのあたりの処理は以下の記事のように相変わらず自分でやる必要がある。

qiita.com

前処理メソッド

様々な前処理用のメソッドがあるので、例示する。なお、例示をしやすいのでpandas.seriesで示しているが通常通りDFに対しても可能。

# データ読み込み
row_text = "    Hi! I've started kàgglé,   and I'm enjoying."
s = pd.Series(row_text)
# =>     Hi! I've started kàgglé,   and I'm enjoying.

# スペースの削除
hero.remove_whitespace(s)
# =>Hi! I've started kàgglé, and I'm enjoying.

# ダイアクリティカルマーク（発音区別符号。àやéなど）の削除
hero.remove_diacritics(s) 
# =>    Hi! I've started kaggle,   and I'm enjoying.

 # 小文字への統一
hero.lowercase(s)
# =>   hi! i've started kàgglé,   and i'm enjoying.

# ストップワードの削除
from texthero import stopwords

default_stopwords = stopwords.DEFAULT # プリセット
custom_stopwords = default_stopwords.union(set(["Hi"]))


hero.remove_stopwords(s, custom_stopwords) 

# =>    I started kàgglé    I enjoying

clean

hero.cleanでは上述含む、一般的によくされるような前処理は一括でおこなってくれる。
具体的には、デフォルトで以下のメソッドがまとめて処理される。

• fillna(s) Replace not assigned values with empty spaces.
• lowercase(s) Lowercase all text.
• remove_digits() Remove all blocks of digits.
• remove_punctuation() Remove all string.punctuation (!"#$%&'()*+,-./:;<=>?@[]^_`{|}~).
• remove_diacritics() Remove all accents from strings.
• remove_stopwords() Remove all stop words.
• remove_whitespace() Remove all white space between words.

# cleanで小文字処理やいらんスペースなどよしなに処理
df['clean_text'] = hero.clean(df['text'])

なお、Textheroではpipe処理をすると可読性よく処理できる。

df['clean_text'] = df['text'].pipe(hero.clean) # pipeでの書き方

# 以下と等価
# df['clean_text'] = hero.clean(df['text'])

また、デフォルトと異なる処理をしたい場合はpipeline引数に処理メソッドをリストで渡す。

from texthero import preprocessing

custom_pipeline = [preprocessing.fillna,
                   preprocessing.lowercase,
                   preprocessing.remove_whitespace]

df['clean_text'] = hero.clean(df['text'], pipeline=custom_pipeline)

# 等価
# df['clean_text'] = df['clean_text'].pipe(hero.clean, custom_pipeline)

解析

named_entitiesで各単語が人物名なのか、地名なのかなど解析してラベリングをしてくれる（正直性能が微妙な気もする)。
こちらも日本語未対応。

具体的なラベルは以下

• PERSON: People, including fictional.
• NORP: Nationalities or religious or political groups.
• FAC: Buildings, airports, highways, bridges, etc.
• ORG : Companies, agencies, institutions, etc.
• GPE: Countries, cities, states.
• LOC: Non-GPE locations, mountain ranges, bodies of water.
• PRODUCT: Objects, vehicles, foods, etc. (Not services.)
• EVENT: Named hurricanes, battles, wars, sports events, etc.
• WORK_OF_ART: Titles of books, songs, etc.
• LAW: Named documents made into laws.
• LANGUAGE: Any named language.
• DATE: Absolute or relative dates or periods.
• TIME: Times smaller than a day.
• PERCENT: Percentage, including ”%“.
• MONEY: Monetary values, including unit.
• QUANTITY: Measurements, as of weight or distance.
• ORDINAL: “first”, “second”, etc.
• CARDINAL: Numerals that do not fall under another type.

s = pd.Series("Yesterday I was in NY with Bill de Blasio")
hero.named_entities(s)[0]

# [('Yesterday', 'DATE', 0, 9),
# ('NY', 'GPE', 19, 21),
#  ('Bill de Blasio', 'PERSON', 27, 41)]

s = pd.Series("Crawford lived in KOBE 10 years ago with Tom")
hero.named_entities(s)[0]

# [('Crawford', 'ORG', 0, 8),
# ('KOBE', 'GPE', 18, 22),
# ('10 years ago', 'DATE', 23, 35),
#  ('Tom', 'PERSON', 41, 44)]

ベクトル変換

おそらく使用のメインになりそうなベクトル変換。現状用意されているのは以下。

• DBSCAN
• k-menas
• MEAN SHIFT
• NMF
• PCA
• tf(TERM FREQUENCY)
• tf-idf
• t-SNE

texthero.org

以下のサンプルはtf-idfとPCAを別カラムとして作成している。なお、次元数は各メソッドのデフォルトを使用。

custom_pipeline = [preprocessing.fillna,
                   preprocessing.lowercase,
                   preprocessing.remove_whitespace]

df['clean_text'] = df['clean_text'].pipe(hero.clean, pipeline=custom_pipeline))) # 前処理
df['tfidf_clean_text'] = hero.tfidf(df['clean_text']) # tf-idfでベクトル化
df['pca_tfidf_clean_text'] = hero.pca(df['tfidf_clean_text']) # PCAで2次元に写像

なお、それぞれ列を作らなくてもいい場合は以下のようにpipeを使った方が可読性が高い。

custom_pipeline = [preprocessing.fillna,
                   preprocessing.lowercase,
                   preprocessing.remove_whitespace]

df['pca'] = (
            df['text']
            .pipe(hero.clean, pipeline=custom_pipeline)) # 前処理
            .pipe(hero.tfidf)# tf-idfでベクトル化
            .pipe(hero.pca) # PCAで2次元に写像
   )

可視化

ワードクラウドと単語の頻度ランキング、散布図が作成できる。散布図はわざわざtexthero経由でやらなくてもいい気もするが、ベクトル化した値は[x, y, z, ...]のように文字列で入るのでseabornなどを使う際は正規表現で抜いてきたり、多次元の場合は列数が増えたりして面倒だが、texthero経由だとそのまま扱うことができる。

texthero.org

散布図

hero.scatterplot(df, 'pca_tfidf_clean_text', color='topic', title='topic')

ワードランキング

hero.top_words(df['clean_text'])

ワードクラウド

日本語の場合などフォント指定したい場合は、fontpathでダウンロードしたフォントを指定する

hero.visualization.wordcloud(df['clean_text'], 
                            colormap='viridis', 
                             width=500,
                             height=400, 
                             background_color='White')

# 日本語の場合はオプションに font_path='fonts/NotoSansCJKJP/NotoSansCJKjp-Regular.otf' などとする

その他

実務で使う際は以下の記事の「テキスト系データの取り扱い」あたり読むとイメージしやすいかも。

www.guruguru.science

参考

qiita.com

qiita.com

www.freecodecamp.org