【第6章】言語処理100本ノックでscikit-learnに入門

この記事では[言語処理100本ノック第6章](https://nlp100.github.io/ja/ch06.html)を用いてscikit-learnを解説します。

まずは`pip install scikit-learn`しておきましょう。

# 50. データの入手・整形
>News Aggregator Data Setをダウンロードし、以下の要領で学習データ（train.txt），検証データ（valid.txt），評価データ（test.txt）を作成せよ．
>
>    ダウンロードしたzipファイルを解凍し，readme.txtの説明を読む．
>    情報源（publisher）が”Reuters”, “Huffington Post”, “Businessweek”, “Contactmusic.com”, “Daily Mail”の事例（記事）のみを抽出する．
>    抽出された事例をランダムに並び替える．
>    抽出された事例の80%を学習データ，残りの10%ずつを検証データと評価データに分割し，それぞれtrain.txt，valid.txt，test.txtというファイル名で保存する．ファイルには，１行に１事例を書き出すこととし，カテゴリ名と記事見出しのタブ区切り形式とせよ（このファイルは後に問題70で再利用する）．
>
> 学習データと評価データを作成したら，各カテゴリの事例数を確認せよ．

この問題はscikit-learnはあまり関係無いので好きな方法で解けば大丈夫です。
とにかくまずはファイルをダウンロードしてreadme.txtを読みましょう。

```python
!wget https://archive.ics.uci.edu/ml/machine-learning-databases/00359/NewsAggregatorDataset.zip
!unzip -c NewsAggregatorDataset.zip readme.txt
```

readmeくらいなら良いですが、圧縮ファイルはなるべく展開せずに扱いたいです。データ本体のzipファイルは[zipfileモジュール](https://docs.python.org/ja/3/library/zipfile.html#zipfile.ZipFile.open)に扱わせると良いです。読み込む手法は何でも良いですが、この場合[pandasを使う](https://pandas.pydata.org/pandas-docs/stable/getting_started/intro_tutorials/02_read_write.html)と楽かなと思います。分割は[sklearn.model_selection.train_test_split()](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.train_test_split.html)にやらせましょう。シャッフルもやってくれます。

初歩的な話ですが、ライブラリ名はscikit-learnですがimportするときのモジュール名はsklearnです。


```python
import csv
import zipfile

import pandas as pd
from sklearn.model_selection import train_test_split


with zipfile.ZipFile("NewsAggregatorDataset.zip") as z:
    with z.open("newsCorpora.csv") as f:
        names = ('ID','TITLE','URL','PUBLISHER','CATEGORY','STORY','HOSTNAME','TIMESTAMP')
        df = pd.read_table(f, names=names, quoting=csv.QUOTE_NONE)

publisher_set = {"Reuters", "Huffington Post", "Businessweek", "Contactmusic.com", "Daily Mail"}
df = df[df['PUBLISHER'].isin(publisher_set)]
df, valid_test_df = train_test_split(df, train_size=0.8, random_state=0)
df.to_csv('train.txt', columns=('CATEGORY','TITLE'), sep='\t', header=False, index=False)
valid_df, test_df = train_test_split(valid_test_df, test_size=0.5, random_state=0)
valid_df.to_csv('valid.txt', columns=('CATEGORY','TITLE'), sep='\t', header=False, index=False)
test_df.to_csv('test.txt', columns=('CATEGORY','TITLE'), sep='\t', header=False, index=False)
```

`pandas.read_table()`はTSVファイルを読み込んで`DataFrame`型オブジェクトを生成します。namesは列名を設定してます。`quoting=csv.QUOTE_NONE`はクォーテーションマークを文字列扱いさせる設定です。`csv.QUOTE_NONE`は`3`と書いても同じです。

（昔`read_table()`はdeprecatedだから`read_csv(sep='\t')`を使えと聞いたんですが、警告が出ないのでundeprecatedになったようです）

`df['PUBLISHER']`の部分は列を抽出する操作で、戻り値は`Series`型になります。pandasの`DataFrame`型がテーブル全体の構造を表現しており、各列が`Series`型で表現されていたわけです。そのメソッド`isin()`は各要素に対する`in`演算の真理値の`Series`を返します。そしてそれをdfのキーであるかのように渡すとTrueの行だけ抽出した`DataFrame`が返ってくるといった感じです。



```python
names = ('CATEGORY','TITLE')
df = pd.read_table('train.txt', names=names, quoting=csv.QUOTE_NONE)
df['CATEGORY'].value_counts()
```




    b    4503
    e    4254
    t    1210
    m     717
    Name: CATEGORY, dtype: int64




```python
df = pd.read_table('test.txt', names=names, quoting=csv.QUOTE_NONE)
df['CATEGORY'].value_counts()
```




    b    565
    e    518
    t    163
    m     90
    Name: CATEGORY, dtype: int64



# 51. 特徴量抽出
> 学習データ，検証データ，評価データから特徴量を抽出し，それぞれ`train.feature.txt`，`valid.feature.txt`，`test.feature.txt`というファイル名で保存せよ． なお，カテゴリ分類に有用そうな特徴量は各自で自由に設計せよ．記事の見出しを単語列に変換したものが最低限のベースラインとなるであろう．

この問題では、抽出した特徴量をベクトル（行列）化しろとは言われていません。のちのちエラー分析などに役立てるために、特徴量を人間に読める形式で保存することが求められていると思われます。

（`scikit-learn`の`Countvectorizer`を使うと特徴量抽出とベクトル化がセットで行われてしまい、この問題と馴染みません）

そこで特徴量を自分で抽出して辞書オブジェクトを作って保存し、次の問題で[Dictvectorizer](https://scikit-learn.org/stable/modules/generated/sklearn.feature_extraction.DictVectorizer.html)を使ってベクトル化するという方針で解きます。辞書のkeyは特徴量の名前、valueはすべて1.0にしておきます。バイナリ素性です。特徴量から辞書を作るこの処理は推論時にも必要になってくるので関数化しておきます。

特徴量の保存形式は特に指定されていませんが、可読性の観点でjsonl形式が良いと思います。

一応コンマやクォートを単語から分離しておきたいです。方法は何でもよいです。tokenizerとしては[spaCy](https://spacy.io/)が有名ですが、この問題では`Countvectorizer`のtokenizerも有力に思います。



```q51.py
import argparse
import json


from sklearn.feature_extraction.text import CountVectorizer


def ngram_gen(seq, n):
    return zip(*(seq[i:] for i in range(n)))


nlp = CountVectorizer().build_tokenizer()

def make_feats_dict(title):
    words = nlp(title)
    
    feats = {}
    for token in words:
        feats[token] = 1.0
    for bigram in ngram_gen(words, 2):
        feats[' '.join(bigram)] = 1.0
    for trigram in ngram_gen(words, 3):
        feats[' '.join(trigram)] = 1.0
    return feats


def dump_features(input_file, output_file):
    with open(input_file) as fi, open(output_file, 'w') as fo:
        for line in fi:
            vals = line.rstrip().split('\t')
            label, title = vals
            feats = {'**LABEL**': label}
            feats.update(make_feats_dict(title))
            print(json.dumps(feats), file=fo)

            
def main():
    parser = argparse.ArgumentParser()
    parser.add_argument('input_file')
    parser.add_argument('output_file')
    args = vars(parser.parse_args())
    dump_features(**args)
    
            
if __name__ == '__main__':
    main()

```

    Overwriting q51.py



```python
!python q51.py test.txt test.feature.txt
```


```python
!python q51.py valid.txt valid.feature.txt
```


```python
!python q51.py train.txt train.feature.txt
```

`ngram_gen()`はngramを作っています。bigramであれば`[[I, am, an, NLPer], [am, an, NLPer]]`を（短い方に合わせて）転置する！というエレガントな方法でやっています。

(ラベルは特徴量では無いですが次の問題が楽になるので書き出しています)

main関数の方で変なことをしていますが、これは4章でも使った[引数リストのアンパック](https://docs.python.org/ja/3/tutorial/controlflow.html#unpacking-argument-lists)によって、辞書でキーワード引数を渡そうとしています。`parse_args()`の戻り値はnamespaceオブジェクトなので、（5章で出てきた）`vars()`で辞書オブジェクトに変換しています。

# 52. 学習
> 51で構築した学習データを用いて，ロジスティック回帰モデルを学習せよ．

まずは51で作ったファイルから特徴量を表す辞書から成るリスト`X`を作ります。それを機械学習モデルに入力するには全特徴量の値を並べたベクトルが必要です。そこで[DictVectorizer()](https://scikit-learn.org/stable/modules/generated/sklearn.feature_extraction.DictVectorizer.html)を使います。`DictVectorizer`のメソッド`fit(X)`で特徴量名とインデックスのマッピングを`X`から獲得し、インスタンス内部の変数に格納します。そして`transform(X)`で`X`を`numpy`行列に変換します。これを一気にやるのが`fit_transform(X)`です。

そして[LogisticRegression()](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html)を使います。インスタンス化して`fit(X, y)`メソッドを呼ぶだけで、インスタンス内部の重みベクトルが学習されます。ハイパラはインスタンス化時に設定します。`X`は行列っぽいもの、`y`はリストっぽいもので、それぞれ長さが一致してれば大丈夫です。

なおロジスティック回帰を知らない人は一度勉強しておきましょう。第8章の布石にもなってます。

学習したモデルは[Model persistence](https://scikit-learn.org/stable/modules/model_persistence.html)を参考に保存します。[joblib.dump()](https://joblib.readthedocs.io/en/latest/generated/joblib.dump.html)を使うときはオプション引数`compress`を指定しないとファイルが大量に生成されるので気を付けましょう。

このとき、特徴量名とインデックスのマッピングも保存しておかないと推論時に困ります。`DictVectorizer`のインスタンスごとdumpしましょう。


```q52.py
import argparse
import json


import numpy as np
from sklearn.feature_extraction import DictVectorizer
from sklearn.linear_model import LogisticRegression
import joblib


def argparse_imf():
    parser = argparse.ArgumentParser()
    parser.add_argument('-i', '--input')
    parser.add_argument('-m', '--model')
    parser.add_argument('-f', '--feats')
    args = parser.parse_args()
    return args


def load_xy(filename):
    X = []
    y = []
    with open(filename) as f:
        for line in f:
            dic = json.loads(line)
            y.append(dic.pop('**LABEL**'))
            X.append(dic)
    return X, y


def main():
    args = argparse_imf()
    X_train, y_train = load_xy(args.input)
    
    vectorizer = DictVectorizer()
    X_train = vectorizer.fit_transform(X_train)
    y_train = np.array(y_train)
    clf = LogisticRegression(random_state=0, max_iter=1000, verbose=1)
    clf.fit(X_train, y_train)
    
    joblib.dump(clf, args.model, compress=3)
    joblib.dump(vectorizer, args.feats, compress=3)

    
if __name__ == '__main__':
    main()

```

    Overwriting q52.py



```python
!python q52.py -i train.feature.txt -m train.logistic.model -f train.feature.joblib
```



# 53. 予測
> 52で学習したロジスティック回帰モデルを用い，与えられた記事見出しからカテゴリとその予測確率を計算するプログラムを実装せよ．

この問題のいう「与えられた記事見出し」は上で作ったテストデータのことを言っているわけでは無く、任意の記事見出しから予測を行うよう言っているように感じます。

保存したモデルを読み込み、`predict(X)`を呼び出せばラベルが、`predict_proba(X)`を呼び出せば予測確率が出てきます。この`X`は入力から特徴量の辞書を作り、52で保存した`Dictvectorizer()`で変換することで得ることができます。

２つのタイトルを入力して`predict_proba()`を適用するとこのような`numpy.ndarray`が返ってきます。




```python
>>> y_proba
```




    array([[0.24339871, 0.54111814, 0.10059608, 0.11488707],
           [0.19745579, 0.69644375, 0.04204659, 0.06405386]])



全ラベルに対する予測確率が出てきていますが、欲しいのは最大値だけだと思います。どうしますか？`ndarray`は[max()メソッド](https://numpy.org/doc/stable/reference/generated/numpy.ndarray.max.html)をもっているようですが…


```python
>>> y_proba.max()
```




    0.6964437549683299




```python
>>> y_proba.max(axis=0)
```




    array([0.24339871, 0.69644375, 0.10059608, 0.11488707])



あとは頑張りましょう。以下、解答例です。


```q53.py
import argparse
import json
import sys


import numpy as np
from sklearn.feature_extraction import DictVectorizer
from sklearn.linear_model import LogisticRegression
import joblib


from q51 import make_feats_dict
from q52 import argparse_imf, load_xy


def predict_label_proba(X, vectorizer, clf):
    X = vectorizer.transform(X)
    y_proba = clf.predict_proba(X)    
    y_pred = clf.classes_[y_proba.argmax(axis=1)]
    y_proba_max = y_proba.max(axis=1)
    return y_pred, y_proba_max


def main():
    args = argparse_imf()
    vectorizer = joblib.load(args.feats)
    clf = joblib.load(args.model)
    X = list(map(make_feats_dict, sys.stdin))
    y_pred, y_proba = predict_label_proba(X, vectorizer, clf)
    for label, proba in zip(y_pred, y_proba):
        print('%s\t%.4f' % (label, proba))

    
if __name__ == '__main__':
    main()

```

    Overwriting q53.py



```python
!echo 'I have a dog.' | python q53.py -m train.logistic.model -f train.feature.joblib
```

    e	0.5441


# 54. 正解率の計測
> 52で学習したロジスティック回帰モデルの正解率を，学習データおよび評価データ上で計測せよ．

手で実装してもよいですが、私は[sklearn.metrics.accuracy_score()](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.accuracy_score.html)に任せます。

`scikit-learn`を習得する上で最も大事なのはこれまでの流れです。

1. 特徴量を抽出し`dict`（を要素とするリスト）型に変換
2. `Dictvectorizer.fit_transform()`で行列に変換
3. `LogisticRegression`など機械学習モデルを選んでインスタンス化
4. `fit(X_train, y_train)`で学習
5. `predict(X_test)`で推論
6. 何らかの方法で評価

これをしっかり押さえておきましょう。



```q54.py
import argparse
import json


import numpy as np
from sklearn.feature_extraction import DictVectorizer
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score
import joblib


from q52 import argparse_imf, load_xy


def predict(args):
    X_test, y_true = load_xy(args.input)
    
    vectorizer = joblib.load(args.feats)
    X_test = vectorizer.transform(X_test)
    y_true = np.array(y_true)
    
    clf = joblib.load(args.model)
    y_pred = clf.predict(X_test)
    
    return y_true, y_pred


def main():
    args = argparse_imf()
    y_true, y_pred = predict(args)
    accuracy = accuracy_score(y_true, y_pred) * 100
    print('Accuracy: %.3f' % accuracy)

    
if __name__ == '__main__':
    main()

```

    Overwriting q54.py



```python
!python q54.py -i train.feature.txt -m train.logistic.model -f train.feature.joblib
```

    Accuracy: 99.897



```python
!python q54.py -i test.feature.txt -m train.logistic.model -f train.feature.joblib
```

    Accuracy: 87.275


# 55. 混同行列の作成
> 52で学習したロジスティック回帰モデルの混同行列（confusion matrix）を，学習データおよび評価データ上で作成せよ．

[sklearn.metrics.confusion_matrix()](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.confusion_matrix.html)に任せます。


```q55.py
from sklearn.metrics import confusion_matrix

from q52 import argparse_imf
from q54 import predict


def main():
    args = argparse_imf()
    y_true, y_pred = predict(args)
    labels = ('b', 'e', 't', 'm')
    matrix = confusion_matrix(y_true, y_pred, labels=labels)
    print(labels)
    print(matrix)

    
if __name__ == '__main__':
    main()
```

    Overwriting q55.py



```python
!python q55.py -i train.feature.txt -m train.logistic.model -f train.feature.joblib
```

    ('b', 'e', 't', 'm')
    [[4499    1    3    0]
     [   2 4252    0    0]
     [   3    1 1206    0]
     [   0    1    0  716]]



```python
!python q55.py -i test.feature.txt -m train.logistic.model -f train.feature.joblib
```

    ('b', 'e', 't', 'm')
    [[529  26  10   0]
     [ 13 503   2   0]
     [ 37  36  89   1]
     [ 19  26   0  45]]


# 56. 適合率，再現率，F1スコアの計測
> 52で学習したロジスティック回帰モデルの適合率，再現率，F1スコアを，評価データ上で計測せよ．カテゴリごとに適合率，再現率，F1スコアを求め，カテゴリごとの性能をマイクロ平均（micro-average）とマクロ平均（macro-average）で統合せよ．

[sklearn.metrics.classification_report()](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.classification_report.html)に任せます。多クラス（シングルラベル）分類で、全クラスについてマイクロ平均を求めると正解率と一致します（[参考](https://analysis-navi.com/?p=553)）。




```q56.py
from sklearn.metrics import classification_report

from q52 import argparse_imf
from q54 import predict


def main():
    args = argparse_imf()
    y_true, y_pred = predict(args)
    print(classification_report(y_true, y_pred, digits=4))

    
if __name__ == '__main__':
    main()
```

    Overwriting q56.py



```python
!python q56.py -i test.feature.txt -m train.logistic.model -f train.feature.joblib
```

                  precision    recall  f1-score   support
    
               b     0.8846    0.9363    0.9097       565
               e     0.8511    0.9710    0.9071       518
               m     0.9783    0.5000    0.6618        90
               t     0.8812    0.5460    0.6742       163
    
        accuracy                         0.8728      1336
       macro avg     0.8988    0.7383    0.7882      1336
    weighted avg     0.8775    0.8728    0.8633      1336
    


# 57. 特徴量の重みの確認
> 52で学習したロジスティック回帰モデルの中で，重みの高い特徴量トップ10と，重みの低い特徴量トップ10を確認せよ．

属性`coef_`に重みがあるのですが、多クラス分類なので重みがクラス数×特徴量ラベル数分あります。4クラスだけなら全て出力してしまいますか。


```q57.py
import joblib
import numpy as np


from q52 import argparse_imf


def get_topk_indices(array, k=10):
    unsorted_max_indices = np.argpartition(-array, k)[:k]
    max_weights = array[unsorted_max_indices]
    max_indices = np.argsort(-max_weights)
    return unsorted_max_indices[max_indices]

def show_weights(args):
    vectorizer = joblib.load(args.feats)
    feature_nemes = np.array(vectorizer.get_feature_names())
    
    clf = joblib.load(args.model)
    coefs = clf.coef_
    y_labels = clf.classes_
    for coef, y_label in zip(coefs, y_labels):
        max_k_indices = get_topk_indices(coef)
        print(y_label)
        for name, weight in zip(feature_nemes[max_k_indices],  coef[max_k_indices]):
            print(name, weight, sep='\t')
        print('...')
        min_k_indices = get_topk_indices(-coef)
        for name, weight in zip(feature_nemes[min_k_indices],  coef[min_k_indices]):
            print(name, weight, sep='\t')
        print()

def main():
    args = argparse_imf()
    show_weights(args)

    
if __name__ == '__main__':
    main()
```

    Overwriting q57.py



```python
!python q57.py -i test.feature.txt -m train.logistic.model -f train.feature.joblib
```


`coef_`全体を`sort`するのではなく上位・下位だけ欲しいと思って、このような回りくどい方法になっています。numpyに`topk()`的な関数は無く、`argpartition()`でソートされていない上位のインデックスを取得するしか無いからです。

# 58. 正則化パラメータの変更
> ロジスティック回帰モデルを学習するとき，正則化パラメータを調整することで，学習時の過学習（overfitting）の度合いを制御できる．異なる正則化パラメータでロジスティック回帰モデルを学習し，学習データ，検証データ，および評価データ上の正解率を求めよ．実験の結果は，正則化パラメータを横軸，正解率を縦軸としたグラフにまとめよ．






```python
import argparse
import json


import numpy as np
from sklearn.feature_extraction import DictVectorizer
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score
import joblib
import matplotlib.pyplot as plt
from tqdm import tqdm


from q52 import load_xy


def get_accuracy(clf, X, y_true):
    y_pred = clf.predict(X)
    return accuracy_score(y_true, y_pred)


X_train, y_train = load_xy('train.feature.txt')
X_valid, y_valid = load_xy('valid.feature.txt')
X_test, y_test = load_xy('test.feature.txt')

vectorizer = DictVectorizer()
X_train = vectorizer.fit_transform(X_train)
X_valid = vectorizer.transform(X_valid)
X_test = vectorizer.transform(X_test)

train_accuracies = []
valid_accuracies = []
test_accuracies = []

for exp in tqdm(range(10)):
    clf = LogisticRegression(random_state=0, max_iter=1000, C=2**exp)
    clf.fit(X_train, y_train)
    train_accuracies.append(get_accuracy(clf, X_train, y_train))
    valid_accuracies.append(get_accuracy(clf, X_valid, y_valid))
    test_accuracies.append(get_accuracy(clf, X_test, y_test))


```


```python
cs = [2**c for c in range(10)]
plt.plot(cs, train_accuracies, label='train')
plt.plot(cs, valid_accuracies, label='valid')
plt.plot(cs, test_accuracies, label='test')
plt.legend()
plt.show()
```



# 59. ハイパーパラメータの探索
> 学習アルゴリズムや学習パラメータを変えながら，カテゴリ分類モデルを学習せよ．検証データ上の正解率が最も高くなる学習アルゴリズム・パラメータを求めよ．また，その学習アルゴリズム・パラメータを用いたときの評価データ上の正解率を求めよ．


アルゴリズム・ハイパラ選択は検証データ上で行い、テストセットチューニングするなということです。ただ前の問題もハイパラ探索でしたし、ここでいう「学習アルゴリズム」とは最適化手法のことというわけではなく、「ロジスティック回帰」ではない手法も検討してねってことですかね？大変なので適当に`sklearn.ensemble.GradientBoostingClassifier`でも使って終わりにします...。


```python
from sklearn.ensemble import GradientBoostingClassifier


clf = GradientBoostingClassifier(random_state=0, min_samples_split=0.01,
                                 min_samples_leaf=5, max_depth=10, 
                                 max_features='sqrt', n_estimators=500, 
                                 subsample=0.8)
clf.fit(X_train, y_train)
valid_acc = get_accuracy(clf, X_valid, y_valid) * 100
print('Validation Accuracy: %.3f' % valid_acc)
test_acc = get_accuracy(clf, X_test, y_test) * 100
print('Test Accuracy: %.3f' % test_acc)
```

    Validation Accuracy: 88.997
    Test Accuracy: 88.548


流行りのGBDTだ！と思ってやってみたものの、ハイパラでかなり性能が変わるので大変でした。時間があれば真面目にgrid searchするのですが...。

ともかく、Pythonで機械学習と言えばscikit-learnということが学べた章だと思います。