nyaaaaaam’s

データサイエンスの技術のこと,仕事のこと,育児のこととかを書いていきます

大阪に引っ越しました

来年から別の会社で働くために,大阪に引っ越してきた.有給消化で12月の第1週から,20足ぐらい早く正月休みに入ったが,そのほとんどが引っ越しとか引越し後の諸々の手続きとかで消滅して,現在に至る.特にトラブルもなく,体調を崩すこともなく,無事に引っ越しを終えることができてよかったと思う.家族もみんな元気です.

大阪は学生時代に住んでいたが,その時と比べて街のつくりとか建物や地理がだいぶ変わっている.写真は,昨日家族で行ったてんしば(天王寺の芝生公園)からの一枚.遠くに通天閣が見える.昔はここはホームレスがたくさんいて,カラオケ大会なんかをやっていたらしいけど,再開発でだいぶ様変わりした,と妻が話していた.

f:id:nyaaaaaam:20181227104112p:plain
遠くに見える通天閣

また,学生時代には気づかなかったが,子供に話しかける人がすごい多かったり,自転車乗りがすごい多かったり,東京とは人の性格もだいぶ違うなぁと思った.

冒頭に書いたが,退職するので,そのうち流行っている退職エントリーも書いてみようと思う.

Geron(2018) Scikit-learnとTensorFlowによる実践機械学習 第3章

標記の本を読んでいますが,身に付けるには読むだけでなくコードを書きながら読むことが一番ということで,演習問題を自分でコード書きながらやってみました.備忘録的にその結果を残します.

scikit-learnとTensorFlowによる実践機械学習

scikit-learnとTensorFlowによる実践機械学習

なお,著者のGeronさんがgithubに文中のコードや演習問題のコードをアップしています.そこを参考にしました(というか,後半のスパム分類はほぼそこのコピペです).

github.com

計算は以前構築したGCP環境でjupyter notebookで実施してます.

nyaaaaaam.hatenablog.com

Q1:MNISTにKNNを試そう

まず必要なモジュールをインポートします.

from sklearn.datasets import fetch_mldata
import numpy as np
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import GridSearchCV, cross_val_score
from sklearn.ensemble import RandomForestClassifier
from sklearn.multiclass import OneVsOneClassifier
from sklearn.linear_model import SGDClassifier
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.decomposition import PCA

from scipy.ndimage.interpolation import shift
import matplotlib.pyplot as plt

不足してるものもあるかもしれないです.次にMNISTデータをフェッチして,訓練・検証データに分割します.フェッチは,上でロードしたfetch_mldataでいけます.今回は著者に習って,あえてtrain_test_splitを使わずやってみました.

#MNISTデータをフェッチする
mnist = fetch_mldata("MNIST original")

#訓練セットとテストセットに分ける
X, y = mnist["data"],mnist["target"]
X_train, X_test, Y_train, Y_test = X[:60000], X[60000:], y[:60000], y[60000:]
shuffle_index = np.random.permutation(60000)
X_train, Y_train = X_train[shuffle_index], Y_train[shuffle_index]

ここでKNNを試そうと思いましたが,なぜだかすごい時間がかかったので諦め,RandomForestを使いました.適当にチューニングしてます.

param_grid = {
    "max_depth":[5,10,100],
    "random_state":[9999]
}

rf_cv = GridSearchCV(RandomForestClassifier(),param_grid=param_grid,n_jobs=-1,verbose=2)
rf_cv.fit(X_train,Y_train)

cross validationでaccuracyの平均とSDを求めます.

rf_cv_accuracy = cross_val_score(rf_cv.best_estimator_,X_train,Y_train,scoring="accuracy",cv=10) #かなりよい
print("CV accuracy (mean): {}".format(np.mean(rf_cv_accuracy)))
print("CV accuracy (sd): {}".format(np.sqrt(np.var(rf_cv_accuracy))))

>CV accuracy (mean): 0.9465169612074671
>CV accuracy (sd): 0.002701355399227736

本で言及されていたのですが,データをスケーリングするとうまく行きやすいらしいので,スケーリングを含むパイプラインを作ってみました.初めて作りましたが,簡単です.

std_rf_pipeline = Pipeline([
    ("scaler", StandardScaler()),
    ("clf", GridSearchCV(RandomForestClassifier(),param_grid=param_grid,n_jobs=-1,verbose=2))
])

std_rf_pipeline.fit(X_train,Y_train)

先ほどと同様にCVでaccuracyの平均とSDを求めます.

rf_cv_accuracy = cross_val_score(std_rf_pipeline,X_train,Y_train,scoring="accuracy",cv=10) #少し向上した
print("CV accuracy (mean): {}".format(np.mean(rf_cv_accuracy)))
print("CV accuracy (sd): {}".format(np.sqrt(np.var(rf_cv_accuracy))))

>CV accuracy (mean): 0.9464666780536921
>CV accuracy (sd): 0.00265692646724053

本当にちょびっとだけ向上しました.

Q2: データの水増し(data augmentation)

訓練データを上下左右4方向に1ピクセルずつずらし,データを水増しする関数を作り,水増しする前の結果と比べてみます.

def shift_image(image, dx, dy):
    image = image.reshape((28, 28))#MNISTの配列を28*28の行列の形に戻す
    shifted_image = shift(image, [dy, dx], cval=0, mode="constant") #ずらす値dx, dyを与えてshift
    return shifted_image.reshape([-1]) #入力とおなじ形に直して返す

こんな感じでシフトできました.

f:id:nyaaaaaam:20181120115541p:plainf:id:nyaaaaaam:20181120115615p:plain

次にこの関数を使って訓練データを水増しします.

X_train_aug_left = [shift_image(image,-1,0) for image in X_train]
X_train_aug_right = [shift_image(image,1,0) for image in X_train]
X_train_aug_up = [shift_image(image,0,1) for image in X_train]
X_train_aug_down = [shift_image(image,0,-1) for image in X_train]

X_train_aug = np.r_[X_train,X_train_aug_left,X_train_aug_right,X_train_aug_down,X_train_aug_up]
Y_train_aug = np.r_[Y_train,Y_train,Y_train,Y_train,Y_train]

先ほど作ったパイプラインを水増しした訓練データで訓練し,結果を比較してみます.

rf_clf_aug = std_rf_pipeline.fit(X_train_aug,Y_train_aug)

aug_rf_cv_accuracy = cross_val_score(rf_clf_aug,X_train_aug,Y_train_aug,scoring="accuracy",cv=10) #結構向上した
print("CV accuracy (mean): {}".format(np.mean(aug_rf_cv_accuracy)))
print("CV accuracy (sd): {}".format(np.sqrt(np.var(aug_rf_cv_accuracy))))

>CV accuracy (mean): 0.9560493405604934
>CV accuracy (sd): 0.00910255937880674

少しだけ向上しました.モデルを色々いじったりするよりは,データを増やしたり,根本的な部分を調整するほうが結果の向上幅は大きいようです.

ついでに,パイプラインにPCAも加えてみます.成分数は適当に100としました.

pca_aug_std_rf_pipeline = Pipeline([
    ("scaler", StandardScaler()),
    ("pca",PCA(n_components=100)),
    ("clf", GridSearchCV(RandomForestClassifier(),param_grid=param_grid,n_jobs=-1,verbose=2))
])

rf_clf_aug_pca = pca_aug_std_rf_pipeline.fit(X_train_aug,Y_train_aug)

aug_pca_rf_cv_accuracy = cross_val_score(rf_clf_aug_pca,X_train_aug,Y_train_aug,scoring="accuracy",cv=10) #結構向上した
print("CV accuracy (mean): {}".format(np.mean(aug_pca_rf_cv_accuracy)))
print("CV accuracy (sd): {}".format(np.sqrt(np.var(aug_pca_rf_cv_accuracy))))

>CV accuracy (mean): 0.9101188311011883
>CV accuracy (sd): 0.01848845832144993

かえって大幅に低下しましたね・・・適当に主成分を選んだのがまずかったと思います.

Q3: Titanicデータセットに挑戦

これはすでにやったことがあるので飛ばしました.

Q4: スパム分類器を作りなさい(難易度高)

今まで扱ったことのないデータなので,Githubのコードを解読しながら進めました.

まずデータを準備します.urllibで取得し,SPAMとHAM(迷惑メールでない)に分けてディレクトリに保存します.

import os
import tarfile
from six.moves import urllib

DOWNLOAD_ROOT = "http://spamassassin.apache.org/old/publiccorpus/" #Apache SpamAssassinのURL
HAM_URL = DOWNLOAD_ROOT + "20030228_easy_ham.tar.bz2"#HAMのファイルの場所
SPAM_URL = DOWNLOAD_ROOT + "20030228_spam.tar.bz2" #SPAMのファイルの場所
SPAM_PATH = os.path.join("datasets", "spam") #データを保存するディレクトリ

def fetch_spam_data(spam_url=SPAM_URL, spam_path=SPAM_PATH): #フェッチする関数 このあたりは決まり手なんだろうな パターン
    if not os.path.isdir(spam_path): #ディレクトリがなければ作る
        os.makedirs(spam_path)
    for filename, url in (("ham.tar.bz2", HAM_URL), ("spam.tar.bz2", SPAM_URL)): #HAMとSPAMのメールについて繰り返し
        path = os.path.join(spam_path, filename)
        if not os.path.isfile(path):
            urllib.request.urlretrieve(url, path) #データを取得
        tar_bz2_file = tarfile.open(path) #書き込み
        tar_bz2_file.extractall(path=SPAM_PATH)
        tar_bz2_file.close()

fetch_spam_data()

次に,具体的なファイル名を取得していきます.なぜか20文字以上のファイル名を抽出していますが,なぜかはわかりませんでした.

HAM_DIR = os.path.join(SPAM_PATH, "easy_ham")#このeasy_hamというディレクトリは展開すると勝手にあるものみたい
SPAM_DIR = os.path.join(SPAM_PATH, "spam")
ham_filenames = [name for name in sorted(os.listdir(HAM_DIR)) if len(name) > 20] #なぜ20以上にするのか?
spam_filenames = [name for name in sorted(os.listdir(SPAM_DIR)) if len(name) > 20]
print("number of hum mails: {}".format(len(ham_filenames)))
print("number of spam mails: {}".format(len(spam_filenames)))

>number of hum mails: 2500
>number of spam mails: 500

ファイル名はこんなかんじです.

ham_filenames[1:10]#こんな感じのファイル名

>['00002.9c4069e25e1ef370c078db7ee85ff9ac',
> '00003.860e3c3cee1b42ead714c5c874fe25f7',
> '00004.864220c5b6930b209cc287c361c99af1',
> '00005.bf27cdeaf0b8c4647ecd61b1d09da613',
> '00006.253ea2f9a9cc36fa0b1129b04b806608',
> '00007.37a8af848caae585af4fe35779656d55',
> '00008.5891548d921601906337dcf1ed8543cb',
> '00009.371eca25b0169ce5cb4f71d3e07b9e2d',
> '00010.145d22c053c1a0c410242e46c01635b3']

次に,保存したemailを読み込んでいきます.

import email #emailを送ったりできるパッケージだが,ここではemailをパースするために使う
import email.policy

def load_email(is_spam, filename, spam_path=SPAM_PATH):
    directory = "spam" if is_spam else "easy_ham" #SPAMかどうかでみるディレクトリを変える
    with open(os.path.join(spam_path, directory, filename), "rb") as f: #ファイルを開いて
        return email.parser.BytesParser(policy=email.policy.default).parse(f) #パースして返す

ham_emails = [load_email(is_spam=False, filename=name) for name in ham_filenames]
spam_emails = [load_email(is_spam=True, filename=name) for name in spam_filenames]

print(ham_emails[200])

>Return-Path: <ilug-admin@linux.ie>
>Delivered-To: zzzz@localhost.netnoteinc.com
>Received: from localhost (localhost [127.0.0.1])
>    by phobos.labs.netnoteinc.com (Postfix) with ESMTP id BD9B044156
>    for <zzzz@localhost>; Wed, 28 Aug 2002 05:47:26 -0400 (EDT)
>Received: from phobos [127.0.0.1]
>    by localhost with IMAP (fetchmail-5.9.0)
>    for zzzz@localhost (single-drop); Wed, 28 Aug 2002 10:47:26 +0100 (IST)
>Received: from lugh.tuatha.org (root@lugh.tuatha.org [194.125.145.45]) by
...

メール本文だけではなく,メールアドレス受信時刻なんかのデータも入っています.

メール本文には,plain textもあれば,htmlで書かれたメールもあります.これをコンテンツタイプと呼びましょう.また,メールの中にメールが入れ子になっていて(例えば返信した時とか)一つのメールが複数のコンテンツタイプを有している場合もあります.それを解析する関数を定義します.

#メールの構造を解析する関数 関数の中でこの関数を再帰的に呼び出す
def get_email_structure(email):
    if isinstance(email, str):#型をチェックする
        return email
    payload = email.get_payload()#付加的情報を除いたもの
    if isinstance(payload, list):#リストだったら(つまりmultipartからなるメール),その要素一つ一つにこの関数をapply. 最終的に文字列になるまでやる
        return "multipart({})".format(", ".join([
            get_email_structure(sub_email)
            for sub_email in payload
        ]))
    else:
        return email.get_content_type()

get_email_structure(ham_emails[1]) #singlepartのメール
> 'text/plain'

get_email_structure(spam_emails[23])#これはmultipart
>'multipart(text/plain, text/html)'

この関数を使ってコンテンツタイプを数え上げ,特徴量とすることを考えます.

from collections import Counter

#コンテントタイプを集計する関数
def structures_counter(emails):
    structures = Counter()
    for email in emails:
        structure = get_email_structure(email)
        structures[structure] += 1
    return structures

ham_emails_types = structures_counter(ham_emails).most_common()
ham_types_val = [obj[1] for obj in ham_emails_types]
ham_types_type = [obj[0] for obj in ham_emails_types]

spam_emails_types = structures_counter(spam_emails).most_common()
spam_types_val = [obj[1] for obj in spam_emails_types]
spam_types_type = [obj[0] for obj in spam_emails_types]

集計結果を描画してみましょう.スパムはマルチコンテンツで,かつHTMLで書かれたメールが多いようです.

#HAM: ほとんどテキスト,署名がある
#SPAM: HTMLがおおい書名があるものはない 画像とかを含んでる
fig, (ax1, ax2) = plt.subplots(ncols=2,figsize=(15,8))
ax1.bar(ham_types_type,height=ham_types_val)
ax1.tick_params(rotation=90)
ax1.set_title("HAM")
ax2.bar(spam_types_type,height=spam_types_val)
ax2.tick_params(rotation=90)
ax2.set_title("SPAM")

f:id:nyaaaaaam:20181120131005p:plain

これらの準備を済ませた上で,学習の準備を進めていきます.まずtrain_test_splitです.

import numpy as np
from sklearn.model_selection import train_test_split

X = np.array(ham_emails + spam_emails)
y = np.array([0] * len(ham_emails) + [1] * len(spam_emails)) #+でつなげる

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

次にhtmlをプレーンテキストに直します.これは,後に文章のステミングをするときに必要になります.

import re
from html import unescape

def html_to_plain_text(html): #HTMLをテキストに戻す
    text = re.sub('<head.*?>.*?</head>', '', html, flags=re.M | re.S | re.I) #ヘッダを落とす
    text = re.sub('<a\s.*?>', ' HYPERLINK ', text, flags=re.M | re.S | re.I) #ハイパーリンクを文字列に
    text = re.sub('<.*?>', '', text, flags=re.M | re.S) #これはなにをしているのか?
    text = re.sub(r'(\s*\n)+', '\n', text, flags=re.M | re.S) #改行をとる
    return unescape(text)

html_spam_emails = [email for email in X_train[y_train==1]
                    if get_email_structure(email) == "text/html"]
sample_html_spam = html_spam_emails[7]
print(sample_html_spam.get_content().strip()[:1000])#変換前

><HTML><HEAD><TITLE></TITLE><META http-equiv="Content-Type" content="text/html; charset=windows-1252">><STYLE>A:link {TEX-DECORATION: none}A:active {TEXT-DECORATION: none}A:visited {TEXT-DECORATION: >none}A:hover {COLOR: #0033ff; TEXT-DECORATION: underline}</STYLE><META content="MSHTML 6.00.2713.1100" >name="GENERATOR"></HEAD>
><BODY text="#000000" vLink="#0033ff" link="#0033ff" bgColor="#CCCC99"><TABLE borderColor="#660000" cellSpacing="0" >cellPadding="0" border="0" width="100%"><TR><TD bgColor="#CCCC99" valign="top" colspan="2" height="27">
...

print(html_to_plain_text(sample_html_spam.get_content())[:1000])#変換後

>OTC
> Newsletter
>Discover Tomorrow's Winners 
>For Immediate Release
>Cal-Bay (Stock Symbol: CBYI)
>Watch for analyst "Strong Buy Recommendations" and several advisory newsletters picking CBYI.  CBYI has filed to be traded >on the OTCBB, share prices historically INCREASE when companies get listed on this larger trading exchange. CBYI is >trading around 25 cents and should skyrocket to $2.66 - $3.25 a share in the near future.
...

さらに,上の関数の出力をテキストに直す関数を定義します.

#さらに上の出力をテキストに直す
def email_to_text(email):
    html = None
    for part in email.walk():
        ctype = part.get_content_type()
        if not ctype in ("text/plain", "text/html"):
            continue
        try:
            content = part.get_content()
        except: # in case of encoding issues
            content = str(part.get_payload())
        if ctype == "text/plain":
            return content
        else:
            html = content
    if html:
        return html_to_plain_text(html)

次に,ntlkを使ってステミングを実行していきます.これはステミングの例です.

import nltk

stemmer = nltk.PorterStemmer()
for word in ("Computations", "Computation", "Computing", "Computed", "Compute", "Compulsive"):
    print(word, "=>", stemmer.stem(word))

>Computations => comput
>Computation => comput
>Computing => comput
>Computed => comput
>Compute => comput
>Compulsive => compuls

ステミングにより,活用形や複数形から,単語の根幹(ステム)を取り出せます.

上の作業をクラスにまとめていきます.パイプラインを上で組んだ時と同じように,ダックタイピングの特性を活かして,fitとtransformのメソッドを作りますが,前処理ではtransformしか必要ないので,fitメソッドはなにもしないように定義してやります.

from sklearn.base import BaseEstimator, TransformerMixin

class EmailToWordCounterTransformer(BaseEstimator, TransformerMixin):
    def __init__(self, strip_headers=True, lower_case=True, remove_punctuation=True,
                 replace_urls=True, replace_numbers=True, stemming=True):
        self.strip_headers = strip_headers 
        self.lower_case = lower_case
        self.remove_punctuation = remove_punctuation
        self.replace_urls = replace_urls
        self.replace_numbers = replace_numbers
        self.stemming = stemming
    def fit(self, X, y=None): #何もしない
        return self
    def transform(self, X, y=None):
        X_transformed = []
        for email in X: #emailに関して繰り返し
            text = email_to_text(email) or "" #メールをテキストにする
            if self.lower_case:
                text = text.lower() #すべて小文字にする
            if self.replace_urls and url_extractor is not None:
                urls = list(set(url_extractor.find_urls(text))) #メール中のurlを抽出する
                urls.sort(key=lambda url: len(url), reverse=True)
                for url in urls:
                    text = text.replace(url, " URL ")
            if self.replace_numbers: #数値を置き換えるか
                text = re.sub(r'\d+(?:\.\d*(?:[eE]\d+))?', 'NUMBER', text)
            if self.remove_punctuation: #句読点を取るかどうか
                text = re.sub(r'\W+', ' ', text, flags=re.M)
            word_counts = Counter(text.split())
            if self.stemming and stemmer is not None:#単語の数え上げ+ステミング
                stemmed_word_counts = Counter() 
                for word, count in word_counts.items():
                    stemmed_word = stemmer.stem(word)
                    stemmed_word_counts[stemmed_word] += count
                word_counts = stemmed_word_counts
            X_transformed.append(word_counts)
        return np.array(X_transformed)

これのtransformメソッドで,次のような単語の数え上げが返ってきます.このメールは3つのメールを内包したマルチコンテントなので,3つの結果が返ってきます.

X_few = X_train[:3]
X_few_wordcounts = EmailToWordCounterTransformer().fit_transform(X_few)
X_few_wordcounts #3つのmultimailを内包しているので,3つのCounterがかえってくる

>array([Counter({'chuck': 1, 'murcko': 1, 'wrote': 1, 'stuff': 1, 'yawn': 1, 'r': 1}),
>       Counter({'the': 11, 'of': 9, 'and': 8, 'all': 3, 'christian': 3, 'to': 3, 'by': 3, 'jefferson': 2, 'i': 2, 'have': 2, 'superstit': 2, 'one': 2, 'on': 2, 'been': 2, 'ha': 2, 'half': 2, 'rogueri': 2, 'teach': 2, 'jesu': 2, 'some': 1, 'interest': 1, 'quot': 1, 'url': 1, 'thoma': 1, 'examin': 1, 'known': 1, 'word': 1, 'do': 1, 'not': 1, 'find': 1, 'in': 1, 'our': 1, 'particular': 1, 'redeem': 1, 'featur': 1, 'they': 1, 'are': 1, 'alik': 1, 'found': 1, 'fabl': 1, 'mytholog': 1, 'million': 1, 'innoc': 1, 'men': 1, 'women': 1, 'children': 1, 'sinc': 1, 'introduct': 1, 'burnt': 1, 'tortur': 1, 'fine': 1, 'imprison': 1, 'what': 1, 'effect': 1, 'thi': 1, 'coercion': 1, 'make': 1, 'world': 1, 'fool': 1, 'other': 1, 'hypocrit': 1, 'support': 1, 'error': 1, 'over': 1, 'earth': 1, 'six': 1, 'histor': 1, 'american': 1, 'john': 1, 'e': 1, 'remsburg': 1, 'letter': 1, 'william': 1, 'short': 1, 'again': 1, 'becom': 1, 'most': 1, 'pervert': 1, 'system': 1, 'that': 1, 'ever': 1, 'shone': 1, 'man': 1, 'absurd': 1, 'untruth': 1, 'were': 1, 'perpetr': 1, 'upon': 1, 'a': 1, 'larg': 1, 'band': 1, 'dupe': 1, 'import': 1, 'led': 1, 'paul': 1, 'first': 1, 'great': 1, 'corrupt': 1}),
>      Counter({'url': 5, 's': 3, 'group': 3, 'to': 3, 'in': 2, 'forteana': 2, 'martin': 2, 'an': 2, 'and': 2, 'we': 2, 'is': 2, 'yahoo': 2, 'unsubscrib': 2, 'y': 1, 'adamson': 1, 'wrote': 1, 'for': 1, 'altern': 1, 'rather': 1, 'more': 1, 'factual': 1, 'base': 1, 'rundown': 1, 'on': 1, 'hamza': 1, 'career': 1, 'includ': 1, 'hi': 1, 'belief': 1, 'that': 1, 'all': 1, 'non': 1, 'muslim': 1, 'yemen': 1, 'should': 1, 'be': 1, 'murder': 1, 'outright': 1, 'know': 1, 'how': 1, 'unbias': 1, 'memri': 1, 'don': 1, 't': 1, 'html': 1, 'rob': 1, 'sponsor': 1, 'number': 1, 'dvd': 1, 'free': 1, 'p': 1, 'join': 1, 'now': 1, 'from': 1, 'thi': 1, 'send': 1, 'email': 1, 'your': 1, 'use': 1, 'of': 1, 'subject': 1})],
>      dtype=object)

最後に上の結果を特徴量ベクトルとして格納するクラスを作ります.

from scipy.sparse import csr_matrix

#上の結果をベクトルに変換するクラス
class WordCounterToVectorTransformer(BaseEstimator, TransformerMixin):
    def __init__(self, vocabulary_size=1000): #引数のvocabulary_sizeを増やせばさらに多くの単語を扱える
        self.vocabulary_size = vocabulary_size
    def fit(self, X, y=None):
        total_count = Counter()
        for word_count in X:
            for word, count in word_count.items():
                total_count[word] += min(count, 10)
        most_common = total_count.most_common()[:self.vocabulary_size]
        self.most_common_ = most_common
        self.vocabulary_ = {word: index + 1 for index, (word, count) in enumerate(most_common)}
        return self
    def transform(self, X, y=None):
        rows = []
        cols = []
        data = []
        for row, word_count in enumerate(X):
            for word, count in word_count.items():
                rows.append(row)
                cols.append(self.vocabulary_.get(word, 0))
                data.append(count)
        return csr_matrix((data, (rows, cols)), shape=(len(X), self.vocabulary_size + 1))#sparse matrixにする ほぼ0なので

こんな感じで動きます.

vocab_transformer = WordCounterToVectorTransformer(vocabulary_size=10)#例のために単語数を少なくする
X_few_vectors = vocab_transformer.fit_transform(X_few_wordcounts)

X_few_vectors.toarray() #10単語が3メール分ある

>array([[ 6,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0],
>       [99, 11,  9,  8,  1,  3,  3,  1,  3,  2,  3],
>       [65,  0,  1,  2,  5,  3,  1,  2,  0,  1,  0]], dtype=int64)

vocab_transformer.vocabulary_ #数え上げた単語の結果

>{'the': 1,
> 'of': 2,
> 'and': 3,
> 'url': 4,
> 'to': 5,
> 'all': 6,
> 'in': 7,
> 'christian': 8,
> 'on': 9,
> 'by': 10}

それでは待ちに待った学習を開始します.まずパイプラインを作りましょう.パイプラインでは

  • EmalToWordCounterTransformerでメールのパース,ステミング,単語の数え上げを行う
  • WordCounterToVectorTransformerで,上の結果(dict形式)をベクトルになおす

を実施します.この前処理パイプラインは,訓練データにのみfitさせることに注意します.

from sklearn.pipeline import Pipeline

preprocess_pipeline = Pipeline([
    ("email_to_wordcount", EmailToWordCounterTransformer()),
    ("wordcount_to_vector", WordCounterToVectorTransformer()),
])

X_train_transformed = preprocess_pipeline.fit_transform(X_train)
X_test_transformed = preprocess_pipeline.transform(X_test) #fit_transformではないことに注意 はまった

ロジスティック回帰を適用し,各種指標を計算します.

from sklearn.linear_model import LogisticRegression, LogisticRegressionCV
from sklearn.model_selection import cross_val_score
from sklearn.metrics import precision_score, recall_score
from sklearn.metrics import roc_curve, auc

cv_log = LogisticRegressionCV()
cv_log.fit(X_train_transformed,y_train)

ml_log = LogisticRegression(C=cv_log.C_[0])#best parameter
ml_log.fit(X_train_transformed,y_train)

score_log = cross_val_score(ml_log, X_train_transformed, y_train, cv=3, verbose=3)

y_pred_log = ml_log.predict(X_test_transformed)
fpr_log, tpr_log, _ = roc_curve(y_test,y_pred_log)

print("CV mean:       {}".format(score_log.mean()))
print("CV sd:         {}".format(np.sqrt(np.var(score_log))))
print("Test accuracy: {}".format(ml_log.score(X_test_transformed,y_test)))
print("Precision:     {}".format(precision_score(y_test,y_pred_log)))
print("Recall:        {}".format(recall_score(y_test,y_pred_log)))
print("AUC:           {}".format(auc(fpr_log,tpr_log)))

>CV mean:       0.9858333333333333
>CV sd:         0.004823265376162628
>Test accuracy: 0.9916666666666667
>Precision:     0.96875
>Recall:        0.9789473684210527
>AUC:           0.9865033871808233

結構いい感じです.XGBoostも試します.

from xgboost import XGBClassifier

ml_xgb = XGBClassifier()
ml_xgb.fit(X_train_transformed,y_train)

score_xgb = cross_val_score(ml_xgb, X_train_transformed, y_train, cv=3, verbose=3)

y_pred_xgb = ml_xgb.predict(X_test_transformed)
fpr_xgb, tpr_xgb, _ = roc_curve(y_test,y_pred_xgb)

print("CV mean:       {}".format(score_xgb.mean()))
print("CV sd:         {}".format(np.sqrt(np.var(score_xgb))))
print("Test accuracy: {}".format(ml_xgb.score(X_test_transformed,y_test)))
print("Precision:     {}".format(precision_score(y_test,y_pred_xgb)))
print("Recall:        {}".format(recall_score(y_test,y_pred_xgb)))
print("AUC:           {}".format(auc(fpr_xgb,tpr_xgb)))

>CV mean:       0.9816666666666666
>CV sd:         0.00523741878749025
>Test accuracy: 0.9916666666666667
>Precision:     1.0
>Recall:        0.9473684210526315
>AUC:           0.9736842105263157

全然チューニングしていませんが,素晴らしい結果です.


かなり長くなりましたが,MNISTデータの水増し,標準化,PCAをパイプラインで構築,加えてスパム分類器の構築を通して自然言語の前処理から機械学習モデルの適用までを学びました.特に後者の自然言語処理は,細かいことはわかってませんが,大まかな作業の流れがつかめたような気がします.英語だと単語の分割が非常に簡単ですが,日本語だと難しそうですね.

Goodfellow et al (2016)のDeeplearningが良かった

f:id:nyaaaaaam:20181010221326p:plain

少し前から,Goodfellow, Bengio & Courville (2016)のDeeplearningを読んでいます.ほとんどの方がよく知っている本かもしれません.表紙が悪夢のようなイラストの本です.

本のサイトはこちらhttps://www.deeplearningbook.org

GithubでPDF版も配布されています(僕はこちらを読んでいます). github.com

こういった技術書は紙のほうが読みやすいという人は,印刷版を買うのも良いと思います.

Deep Learning (Adaptive Computation and Machine Learning)

Deep Learning (Adaptive Computation and Machine Learning)

もともと,この本はよく通勤中やランニング中に聞いているポッドキャストのひとつ,misreading chatで紹介されていたもので,興味を持っていました.ひょんなことから会社で印刷版を買うことになり,自宅に持ち帰ったりするのが面倒だったので,結局PDF版を読んでいます.

先日,第1部のApplied math and Machine Learning Basicsを読み終わりまして,結構感動したことが多かったので,その感想を雑にまとめておきます.

数学的道具の説明が必要かつ十分

本書では,機械学習全般,特にDeeplearningの歴史を振り返りつつ,近年の動向を概観した後に,線形代数の基礎的な話が始まります.大学院の研究で普段線形代数をめちゃくちゃ良く使うので,こういう線形代数の話は耳にタコができるほど聞いてきたのですが,シンプルかつ明快な線形代数のイントロがとても心地よかったです.特に,固有値分解の一般化としての特異値分解,ムーアペンローズ一般逆行列などの説明が好みでした.最後に,学んだことを使って主成分分析を導出するエンディングも,素敵です.

その後,線形代数と同じく重要な数学的道具として,確率論や情報理論,最適化の話が続きます.だいたいこういった数学的道具の説明は「小難しくて理解できない割に,本論ではあまり使わない」ものが多いように感じていましたが,説明も簡潔かつ直感的に理解しやすい,必要十分な説明に留められていると感じました.

新しく学んだこと1: CapacityとVC次元

5.2では,機械学習モデルのUnder/Overfittingの考え方を紹介する際に,モデルの複雑性=capacityという概念を導入して説明しています.良いモデルとは,データに対して十分なcapacityを持つものだと考えられますが,そのcapacityを評価するための重要な概念として,VC次元(Vapnic-Chervonenkis dimension)を紹介しています.これは「任意のラベル付けで完全に分離できるデータポイントの数」であり,たとえば2次元空間上のデータポイントなら,3点までならどんなラベル付けをしたとしても完全に分離できるデータポイントは3です.よってこのときのVC次元は3になります.このVC次元は,無限のサイズの仮説集合に対して,まぁまぁいい学習を行うために必要な訓練データの数を導く上で必要になるらしいです(参考:計算論的学習理論入門 -PAC学習とかVC次元とか-).

新しく学んだこと2: ベイズ推定と正則化

f:id:nyaaaaaam:20181011181315j:plain 線形回帰モデルのベイズ推定を解説している箇所があります.ここではいろいろと式を展開してパラメータの事後分布を求めていくと,Ridge回帰のclosed formの解に事後平均が一致することが示されています.つまり,ベイズ推定とは本質的には正則化である,という結論らしく,この点は今まで知らなかったポイントでした.

機械学習アルゴリズムの構築

データに対してどのような機械学習アルゴリズムを適用するかを表現するためには,次の4つのbuilding blockを明らかにせよ,という主張です. - どのデータを使うか? - 適用しようとしているモデルは何か? - 損失関数は何か? - 最適化アルゴリズムは何か? たとえば,同じデータであっても,モデルを変えればまったく違ったアルゴリズムになりますし,最適化アルゴリズムを変えれば違った解が得られるかもしれません.逆に言えば,4つのうちどれかがはっきりと決まっていないと,ほかの人にアルゴリズムを伝えるときに,正しく伝えられない=再現性が取れない可能性があります.この辺は,実務家として理解しておくべきことだと思いました.


書きながら「理解が甘いなぁ」というところも再発見しましたが,何度も書いているように,必要十分な知識を簡潔に得ることのできる本(というかまだイントロですが)だと感じました.これから本題のdeeplearningに入っていきますが,楽しみです.

また感想がまとまったら書きます.

参考文献

Goodfellow, I., Bengio, Y., Courville, A., & Bengio, Y. (2016). Deep learning (Vol. 1). Cambridge: MIT press.

普段の時間の使い方と研究の記録

f:id:nyaaaaaam:20180914230525p:plain

僕は普段会社員として働くかたわら,2年前からとある大学の博士課程に入学して,博士学生として研究もやっています.そのような二足のわらじを履いてはや二年,いろいろ試行錯誤しつつたどり着いた時間の使い方などについて書きたいと思います.

一日の流れ

だいたいこんな感じで仕事と研究を両立させています.

  • 7:00 起床(たまに早く起きてジョギングなど運動)
  • 8:00 出社
  • 18:00 帰宅
  • 21:00 食事や入浴を済ませ研究開始
  • 23:30 研究終了
  • 24:00 一日のまとめなどして就寝

割とホワイトな会社かつ「遅くまで会社にいないやつ」と周囲に印象づけることに気をつけて会社員生活を歩んできたので,結構早く帰宅しています.また,専攻は多変量解析法の理論的研究(主成分分析や因子分析などの多変量解析法をいろいろ魔改造する研究)なので,基本的に自宅で研究できます.先生とはたまに連絡を取りつつ,主に論文を書いたりしています.

今はこんな感じの時間の使い方で落ち着いているのですが,こういうパターンを発見するまでにいろいろと紆余曲折がありました.その中で学んだことを以下にまとめます.

遅くまで研究しすぎない

院生になりたての頃は「研究頑張るぞ」と息巻いて,深夜2時や3時まで研究をがんばる夜も多くありました.しかし今年で30代に入った事もあってか,遅くまで研究すると次の日の仕事にかなり響きます.何日か繰り返すと慣れるかなと思ったのですが,全く慣れず,それどころか次の日の夜は早く寝たい気持ちになってしまうので,結局研究に避ける時間が少なくなる事に気づきました.そんなこともあって,最近は(よほど何かの締切に追われていない限りは)24時までには研究を終えて床につくようにしています.「細く長く」続けていくスタイルが僕にはあっているようです.

子供が生まれると研究に割ける時間が増える

子供が生まれる前は,夕食後妻とダラダラしていたので,結局研究を開始できるのは妻が眠って以降の22時頃からでした.しかし子供が生まれて以降,妻と子供が遅くとも21時頃には,規則正しく就寝するようになったので,その分長く研究時間が取れています.周囲には「子供が生まれると自分の時間がなくなるよね」と言われていましたが,実際生まれるとそんなことはないなぁと思いました.もちろん,夜泣きがひどいときにはあやすのを手伝ったりもします.

ちゃんと研究してると仕事にもいい影響が出る

夜やっている研究と昼の仕事とが比較的近いせいもあり,夜しっかり研究していると,その時の気づきを仕事に応用したり,反対に仕事の気づきを研究に応用したりと,二足のわらじならではのメリットも大きいなと思います.頭も切り替わっていい感じです.

記録を取ることは大事

研究生活を始めた当初から,画像みたいな感じで研究ノートをEvernoteにつけています.ほぼ日記みたいな感じです. f:id:nyaaaaaam:20180914231047p:plain 毎日やったことや,やり残したこと,感想などを書き付けておくと,達成感もあるし,タスク管理も簡単にできるので,重宝しています.1年前の研究ノートとか見ると,結構懐かしい気持ちになります.WunderlistとかAsanaとか,タスク管理やプロジェクト管理のアプリを使ってみたこともあったけど,結局めんどくさくてあんまり使わず,Evernoteに落ち着いています.


また思いついたら書きます(ブルーライトカットメガネを会社に忘れてきて,目がしょぼしょぼして疲れる...).

【GCP+Bitbucket】kaggleのための環境構築

先日,kaggle(http://kaggle.io)のための環境を構築しました.その時の作業をまとめておきます. 計算機環境としてのGoogle Cloud PlatformGCP)と,バージョン管理のためのBitbucketを使っています.これらは

  • GCP : 見た目がおしゃれ,300USDのクレジットがついててしばらく無料で使える
  • Bitbucket:プライベートレポジトリが無料で使える,こちらもおしゃれ

ってことで選びました.

構築にあたっては,Takumi Satoさんの動画を参考にしました.Youtubeで説明してくださっていて,とてもわかりやすかったです.

細かい手順は上の動画を参考にしてください.

GCPを準備する

上のGCPのリンクから,VMインスタンスを作成します.けちってあまり豪華な構成にはしていないです. f:id:nyaaaaaam:20180903163952j:plain

OSは仕事でも使うことの多いUbuntuを使ってみました.セキュリティは大事なので,公開鍵認証にしています.

Bitbucketを準備する

Githubを使うと,kaggleの利用規約には「チーム内の情報シェアなどを除き,kaggleのフォーラムを通して情報をシェアしなければならない」というものがあるらしく,Githubの公開レポジトリを使うと,この規約に違反する可能性があります.そのため,無料でプライベートレポジトリが使えるBitbucketを使います. 普通にIDをを取得してレポジトリを作成しておきましょう.

GCPへの接続

毎度GCPのページにアクセスして,そこからターミナルを立ち上げてログインすることもできますが

  • 微妙にレイテンシがある
  • ログインに時間がかかる

といった理由があるので,ターミナルから接続することをおすすめします.ターミナルでつぎのようにしましょう.

ssh -i (秘密鍵のパス) (ユーザ名)@(VMインスタンスの外部IPアドレス) -p (ポート番号,あれば)

こうすれば,ターミナルから接続できて便利です.

Jupyter notebookを使えるようにする

データサイエンスでは,対話的にコードが書けるJupyter Notebookを使う人が多いんじゃないかと思います.僕もそうなので,kaggle用の環境には,Jupyter notebookが必須です.

GCPインスタンスを立ち上げたあと,Anacondaを以下のようにしてインストールします.

wget = https://repo.anaconda.com/archive/Anaconda3-5.2.0-Linux-x86_64.sh
sh Anaconda3-5.2.0-Linux-x86_64.sh

これでJupyter notebookのインストールは完了です.

あとはこの記事に従って,ポートを開放したりすれば,Jupyter notebook環境の出来上がりです.

ディレクトリ構成

データサイエンスプロジェクトのディレクトリ構成をどうするか問題に関して,今回は以下の様な構成をとっています.

├── README.md
├── input: データセットを入れておく.ファイルは絶対に編集しない.
│   ├── train.csv
│   ├── test.csv
│   └── ...
├── notebooks: 実験したnotebookを入れておく
│   ├── 20180901_EDA.ipynb
│   └── ...
├── results: 予測結果とか,submit用のcsvを入れておく
│   ├── 20180901_predict_LASSO.csv
│   └── ...
└── trained_models: 訓練済みのモデルを.pklで入れておく
    ├── 20180901_LASSO.pkl
    └── ...

まずinput/には,kaggleなどからDLしたデータそのものを入れておきます.特徴量エンジニアリングなんかでいろいろとデータを弄る場合は,このファイルをいじるのではなく,別のディレクトリを作って編集して,もとのデータセットを編集しないことをおすすめします(あとで元データが何なのかわからなくなるため).また,trained_modelsには,joblib.dumpで保存したモデルを入れておきます.CVなどですごく時間がかかる場合,結果を保存しておくといろいろとあとで時間の節約になるので.

こういったディレクトリ構成ではcoockiecutterというモジュールを使うと良いっぽいですが,そこまで細かいディレクトリ構成は必要ないと感じたので(今のところは),今回はオリジナルの構成で組んでいます.

github.com

今取り組んでるコンペ

今は,不動産の特徴からその価格を予測する,House Pricesというコンペに取り組んでいます.データセットはある程度構造化されていて(欠損値などはもちろんありますが),回帰をいろいろと試してみる問題です.なにかいい結果が出たら,Kernelでも書こうかなと思っています.

https://www.kaggle.com/c/house-prices-advanced-regression-techniques

追記

GCPのストレージからVMにファイルを転送する.

gsutil -m cp gs://ストレージ上のパス VM上のパス

はじめてのエントリー

f:id:nyaaaaaam:20180902153232p:plain

はじめまして,はてなブログでブログを始めてみました.

ブログで書いていくこと

以下のことについて,不定期で書いていきます.

  • データサイエンス関連の技術の話:Python, R, Kaggle...
  • 仕事の話:研究や開発,転職のことなど
  • 育児の話
  • 趣味の話 今後追加したり削除したりするかもしれません.

ブログを書こうと思い立ったきっかけ

現在僕はとあるメーカーの研究所で働いています.研究所内でデータサイエンティスト的なことをしています.しかしいろいろと思うところもあり(この辺の事情はまた別のエントリーで書こうと思います),2019年の1月からテック系に転職して,本格的にエンジニアとして仕事を始めることにしました.

「エンジニアってなんとなくブログを書くもの」とぼんやり思っていたのですが,セルフブランディングとか,アウトプットの場を持つことで勉強が捗るなど,いろいろとエンジニアにとっていいことがあるということなので,これを機にはじめてみることにしました.

続くかもしれないし続かないかもしれませんが,よろしくおねがいします!