Obey Your MATHEMATICS.

機械学習関連の純粋数学や実験など

EdwardによるDeep Beta Distribution(深層ベータ分布)モデル

こんにちは。お久しぶりの投稿です。

来週末に開催される db analytics show case Sapporo

www.db-tech-showcase.com


と言うイベントで講演する事になってまして、ベイズ統計やMCMCの基本的なところからEdwardのデモまでやっていく予定なのですが

ただ基本的なモデルを紹介するのもおもしろくないので、僕なりの新しい深層学習+確率モデリングなモデルを考えましたので紹介したいと思います。

EdwardやMCMCの基本的なところについては弊社のブログに寄稿した以下の記事が詳しいので合わせてそちらも御覧ください↓

data.gunosy.io

§1. お気持ち

例えば、ユーザーが付ける商品のレーティングを予測したいと考えます。

これは通常ターゲット { y } が閉区間{ [ 0, 1 ] } に値を取る回帰問題として定式化されます。*1

与えられた特徴量に対して、予測値 { y }を出力するわけですが、それは決定論的に決まるものなのか疑問があります。


ユーザーのレーティングであれば、気分・時間・各ユーザーが商品に対する持つ偏見(例えば◯◯製の製品は嫌い)等、ランダムネスがかなりあるように思われます。


つまり、ユーザーの特徴量と商品の特徴量が与えられた時に、閉区間{ [ 0, 1 ] }上の確率分布が定まりそこからレーティングがサンプリングされるようにモデリングしたほうが適切なのではないのかと考えるわけです。

§2. 定式化

深層学習の優秀さは言うまでもないので、それと上の確率モデリングの問題を上手く融合させます。

区間{ [ 0, 1 ] }上の連続値確率分布として、代表的なベータ分布があると思いますが

2つの正のパラメータを変化させる事によって様々な閉区間{ [ 0, 1 ] }上の確率分布を得る事が出来る事を思い出しましょう。


これを上手く利用して、特徴量 {x } に対して ベータ分布の2つのパラメータが出力されるようなニュールネットを作ります。

これにより、特徴量が与えられた際に様々な形状の閉区間{ [ 0, 1 ] }上の確率分布が出力されるモデルが得られます。

*2としてはこんな感じです。


f:id:mathetake:20170624134023p:plain

§3. 簡単な実験

今回のデータセットは 特徴量{x=(x_i)_i \in \mathbb{R}^{X_{DIM}} } に対してベータ分布のパラメタ {\alpha_x, \beta_x > 0}

{\alpha_x = \dfrac{1}{X_{DIM}}\sum \dfrac{1}{ 1 + \exp(- x_i) }}
{\beta_x = \dfrac{1}{X_{DIM}}\sum |x_i| }

により定まりそこからターゲット {y}がサンプリングされていると言う非線形なお遊びデータセットを使い、4層のニューラルネットを使ってモデリングしていきます。



まず必要なライブラリをimportして、定数とデータセットを準備します。

import numpy as np
import tensorflow as tf
import edward as ed
from edward.models import Normal, Empirical, 

# consts
N = 10000  # training sample size
N_test = 1000 # test sample size
X_DIM = 500  # input dim
MINI_BATCH_SIZE = 100  # batch_size
N_ITER = 10000 # number of iteration of MCMC
np.random.seed(seed=1)

# true_alpha =  softmax(x).mean()
# true_beta = |x|.mean()
def generate_samples(n):
    x = np.random.normal(size=[n, X_DIM])
    true_alpha = (1.0 / (1.0 + np.exp(-x))).mean()
    true_beta = np.absolute(x).mean()
    y = np.random.beta(a=true_alpha, b=true_beta, size=[n])
    return x, y

X_data, Y_data = generate_samples(N)
X_test, Y_test = generate_samples(N_test)

サンプルサイズと次元から、通常のMCMCでは歯が立たないのでSGMCMCを使っていくことにします。
そのためにunbiasedなミニバッチを生成する関数を定義しておきます。

# a function for generating unbiassed mini batches
def next_batch(mini_batch_size=MINI_BATCH_SIZE):
    indexes = np.random.randint(N, size=mini_batch_size)
    return X_data[indexes], Y_data[indexes]

ここからEdwardでモデリングを開始します。
まずニューラルネットのパラメタを準備して事前分布を定義し、またそれらを近似するための経験分布を定義します。

# define prior distributions over parameters
W_0 = Normal(loc=tf.zeros([X_DIM, 50]), scale=tf.ones([X_DIM, 50]))
W_1 = Normal(loc=tf.zeros([50, 10]), scale=tf.ones([50, 10]))
W_2 = Normal(loc=tf.zeros([10, 2]), scale=tf.ones([10, 2]))
b_0 = Normal(loc=tf.zeros(50), scale=tf.ones(50))
b_1 = Normal(loc=tf.zeros(10), scale=tf.ones(10))
b_2 = Normal(loc=tf.zeros(2), scale=tf.ones(2))

# define empirical distributions for prameter inference
q_W0 = Empirical(params=tf.Variable(tf.random_normal([N_ITER, X_DIM, 50])))
q_W1 = Empirical(params=tf.Variable(tf.random_normal([N_ITER, 50, 10])))
q_W2 = Empirical(params=tf.Variable(tf.random_normal([N_ITER, 10, 2])))
q_b0 = Empirical(params=tf.Variable(tf.random_normal([N_ITER, 50])))
q_b1 = Empirical(params=tf.Variable(tf.random_normal([N_ITER, 10])))
q_b2 = Empirical(params=tf.Variable(tf.random_normal([N_ITER, 2])))


次にモデルを構築します

def build_deep_Beta_distribution(x):
    h_1 = tf.tanh(tf.matmul(x, W_0) + b_0)
    h_2 = tf.tanh(tf.matmul(h_1, W_1) + b_1)
    h_3 = tf.nn.softplus(tf.matmul(h_2, W_2) + b_2)
    alpha = h_3[:, 0]
    beta = h_3[:, 1]
    return Beta(concentration1=alpha, concentration0=beta)

# prepare placeholders for data points
x_ph = tf.placeholder(tf.float32, [MINI_BATCH_SIZE, X_DIM])
y_ph = tf.placeholder(tf.float32, [MINI_BATCH_SIZE])

# build our model
y = build_deep_Beta_distribution(x_ph)


今回は、推定アルゴリズムとしてSGHMCを使います。SMHMCはHMCを単純に確率勾配にするだけではなく(この場合収束が保証されない)、補助パラメータを導入してサンプリングしていく手法です。
[1402.4102] Stochastic Gradient Hamiltonian Monte Carlo

# initialization of parameter inference
latent_vars = {W_0: q_W0, W_1: q_W1, W_2: q_W2,
               b_0: q_b0, b_1: q_b1, b_2: q_b2}
SGHMC = ed.SGHMC(latent_vars=latent_vars, data={y: y_ph})
SGHMC.initialize(scale={y: float(N) / MINI_BATCH_SIZE}, step_size=0.0001, n_iter=N_ITER)

init = tf.global_variables_initializer()
init.run()
for _ in range(N_ITER):
    X_batch, Y_batch = next_batch(MINI_BATCH_SIZE)
    info_dict = SGHMC.update(feed_dict={x_ph: X_batch, y_ph: Y_batch})
    SGHMC.print_progress(info_dict)

§4. 精度評価

Edwardでの精度評価は簡単で、テストデータ用のplaceholderを用いて edward.copy メソッドを使って推定されたパラメータの平均値を使ったグラフを新たに構築した後 edward.evalutate メソッドにデータを食わせるだけです。

# criticize our model
x_test_ph = tf.placeholder(tf.float32, [N_test, X_DIM])
y_test = build_deep_Beta_distribution(x_test_ph)
dict_swap = {W_0: q_W0.mean(), W_1: q_W1.mean(), W_2: q_W2.mean(),
             b_0: q_b0.mean(), b_1: q_b1.mean(), b_2: q_b2.mean()}
y_post = ed.copy(y_test, dict_swap=dict_swap)
log_loss = ed.evaluate(metrics='mean_squared_error', data={x_test_ph: X_test, y_post: Y_test})
print('[deep Beta distribution] mean_squared_error on test data: ', log_loss)


簡単に他のモデルと比較するために、線形回帰でもモデリングしておきます:

from sklearn.linear_model import SGDRegressor
from sklearn.metrics import mean_squared_error
linear = SGDRegressor(loss='squared_loss')
linear.fit(X_data, Y_data)
linear_loss = mean_squared_error(Y_test, linear.predict(X_test))
print('[linear regression] mean_squared_error on test data: ', linear_loss)


結果はこんな感じでした。

[deep Beta distribution] mean_squared_error on test data:  0.108737
[linear regression] mean_squared_error on test data:  0.132299

§5. おわりに

ものすごく雑に書いてしまいましたが、このようにEdwardを使えば自由自在に深層学習と確率モデリングを組み合わせる事ができます。

日本人のコミュニティはほとんどない状態ですので、この記事を読んで少しでもおもしろいなあと思ったら Edwardのチュートリアルをやってみてください↓

Edward – Getting Started


それでは。

*1:正直なところ普通はマルチクラス分類にするんでしょうが、分かりやすい例を出したかったのでこうしてます

*2:http://www.st.keio.ac.jp/learning/1403.html, https://ja.wikipedia.org/wiki/%E3%83%99%E3%83%BC%E3%82%BF%E5%88%86%E5%B8%83

コサイン類似度&L^2ノルムの変動を用いた特徴語抽出

こんにちは

Twitter就活コンサルタントと巷で噂のマスタケです

学生の頃みたいなクソみたいな記事は書いてはイケないと言うプレッシャーがあり、しばらくご無沙汰してましたが、書きます。(今回の記事がクソ記事ではないと言う意味ではありません。クソ記事です)

§1. はじめに

今日のブログネタは次の2つの論文

[1]抽出型文書要約における分散表現の学習 —文書と要約の距離最小化—

[2]Summarization Based on Embedding Distributions


を読んで思いついた事の実験の報告と言った感じです。

これらの論文では、(word2vecなどの)単語の分散表現{\mathcal{M}: \{ \mbox{単語}\} \rightarrow \mathbb{R}^d } を用いて次のように文章ベクトルを定義し*1

{ \mathcal{F}:\{\mbox{文章}\} \rightarrow \mathbb{R}^d   , \ \  \mathcal{D} \mapsto \dfrac{1}{|\mathcal{D}| }\sum_{w \in \mathcal{D}} \mbox{tf-idf}(w) \mathcal{M}(w)   }

コサイン類似度を用いて2つの文章間の類似度を測り、それをベースに自動要約のアルゴリズムを提案しています*2:

f:id:mathetake:20170504175842p:plain


この文章ベクトルの良いところは、tf-idfベースの類似度算出の欠点である、互いに出現しない単語の類似度が考慮される点にあります。



これらの論文を見つける前に、最近仕事でこんな感じの似たようなベクトルを作って要約とは別の目的で使っていたのですが、なにかこの路線でおもしろい事が出来ないかなあと考えていて

特徴語抽出にこのベクトルを使えるのでは??

とアイデアを思いついたのでその定式化と実験を報告します。

§2. 提案手法

今回のアイデアはいたってシンプルで

特徴語 = その単語を除いた文章のベクトルと、元の文章のベクトルの類似度の乖離が大きい

と考える事です。より正確に、次の2つの関数

{ \mathcal{G}_{L^2}: \mathcal{D} \mapsto \mbox{argmax}_{w \in \mathcal{D}} \| \mathcal{F}(\mathcal{D}) - \mathcal{F}(\mathcal{D} \backslash \{w\} ) \|_2 }

{ \mathcal{G}_{cosine}: \mathcal{D} \mapsto \mbox{argmin}_{w \in \mathcal{D}} \mbox{cosinesim} \biggl(\mathcal{F}(\mathcal{D}), \mathcal{F}(\mathcal{D} \backslash \{w\} ) \biggr) }


が最も特徴的な単語を抽出するであろうと考えます。1つ目の式ではL^2距離の乖離が激しい単語こそ特徴語であると考え、2つ目の式ではコサイン類似度の乖離が激しい単語こそ特徴語であると考えています。

§3. データセットと分散表現モデル

今回実験のために使ったデータセットはKaggleのAmazon Fine Food Reviewsです:

Amazon Fine Food Reviews | Kaggle


また、分散表現モデルにはGoogleが公開している学習済みのword2vecのモデルを使いました:

Google's trained Word2Vec model in Python · Chris McCormick


§4. 実験結果

実験では、それぞれの手法に対してトップ5の特徴語候補を抽出してみました。

トキトウに選んだいくつかの文章に対する結果は以下のとおりです:

元の文章: If you are looking for the secret ingredient in Robitussin I believe I have found it. I got this in addition to the Root Beer Extract I ordered (which was good) and made some cherry soda. The flavor is very medicinal.

TOP1 TOP2 TOP3 TOP4 TOP5
コサイン類似度 medicinal secret cherry soda ingredient
L^2ノルム medicinal secret cherry soda ingredient


元の文章: Good oatmeal. I like the apple cinnamon the best. Though I wouldn't follow the directions on the package since it always comes out too soupy for my taste. That could just be me since I like my oatmeal really thick to add some milk on top of.

TOP1 TOP2 TOP3 TOP4 TOP5
コサイン類似度 soupy directions cinnamon apple follow
L^2ノルム soupy cinnamon directions oatmeal oatmeal


元の文章: My English Bulldog had skin allergies the summer we got him at age 3. The vet recommended we wean him off the food his previous owner gave him (Iams Lamb and Rice) and onto a new kind. This was the second one we tried, and it has been working ever since. It's for dogs that need a limited diet who can be sensitive to additives and proteins commonly found in commercial dog food (like chicken or beef).

TOP1 TOP2 TOP3 TOP4 TOP5
コサイン類似度 wean proteins commonly additives allergies
L^2ノルム wean proteins additives commonly allergies


元の文章: To me, these are nothing like the regular Altoids and are not breath mints. They are pleasant-tasting little candies in a cute convenient tin, and that's as far as it goes. The mintiness is just not strong, and the wintergreens are definitely weaker than the peppermint minis. I'm not a dragon-breath person, but still, one of these mints is too small to have any effect on my breath. Four or five will freshen my breath for a short while - maybe 15 minutes. At this point, I think the Icebreakers Frost mints are the best as sugar free breath mints.

TOP1 TOP2 TOP3 TOP4 TOP5
コサイン類似度 freshen dragon minis weaker peppermint
L^2ノルム freshen dragon minis mints mints


元の文章: i cannot live without this citron falksalt, it is wonderful on watermelon, fish, almost any food requiring salt and in my world it all does. just try it on red juicy watermelon, wonderful experience.

TOP1 TOP2 TOP3 TOP4 TOP5
コサイン類似度 citron requiring juicy watermelon watermelon
L^2ノルム citron watermelon watermelon requiring juicy


元の文章: Kettle Brand chips used to be so good...oily, crunchy, flavorful. I suspect the company has been bought out and the recipe has been changed for the worse. Now they're no better than any other big name brand chip. Try the Good Health Kettle Style Olive Oil chips instead. They are as good as the Kettle Brand once was. R.I.P., Kettle Brand chips. :(

TOP1 TOP2 TOP3 TOP4 TOP5
コサイン類似度 oily suspect flavorful recipe worse
L^2ノルム oily flavorful suspect recipe worse



こんな感じになりました。
コサイン類似度でもL^2ノルムでもほとんど変わっていませんが、なんとなく特徴語は引っ張ってこれているかなあと言う雰囲気です。

§5. 所感

特徴語抽出は出来ていると言う雰囲気はありますが、果たして…と言う感じです。
そもそもtf-idfベースでの抽出と結果が変わらないのでは、と言う疑問がありますが、小学生の自由研究なのでめんどくさいので比較はしません。

お粗末さまでした。

*1:少し定義が違いますが、ここではこの定義を使います

*2:この図は[1]より

【Recsys2016】Adaptive, Personalized Diversity for Visual Discovery (AmazonStream)に関するメモ

この記事は次の論文:

Adaptive, Personalized Diversity for Visual Discovery

Choon Hui Teo Amazon, Palo Alto, CA, USA
Houssam Nassif Amazon, Seattle, WA, USA
Daniel Hill Amazon, Palo Alto, CA, USA
Sriram Srinivasan UC Santa Cruz, Santa Cruz, CA, USA
Mitchell Goodman Amazon, Seattle, WA, USA
Vijai Mohan Amazon, Seattle, WA, USA
S.V.N. Vishwanathan Amazon & UC Santa Cruz, Palo Alto, CA, USA


についての個人的なメモです。*1

§1. 概要と背景

この論文は、Recsys2016 のショートペーパー部門の受賞論文。*2

AmazonStreamと呼ばれるサービス(日本にはないっぽい?)を題材*3にして、推薦システムにつきまとう推薦アイテムの種類の偏りを劣モジュラ関数最適化問題に焼き直して上手くやったと言うもの。多様性を持たせた推薦システムの構築。

好きなお店にウィンドウ・ショッピングしていく感覚でアイテムをユーザーに推薦したい ---> 豊かなユーザー体験へ。

ステップは次の3つに分けられる
◯CTR予測モデル
◯多様性問題の定式化
◯ユーザーのモデリング(?)

§2. CTRのモデリングについて-その1-

基本的に次のMicrosoftの論文と同じ。

Web-Scale Bayesian Click-Through Rate Prediction for Sponsored Search Advertising in Microsoft’s Bing Search Engine - Microsoft Research

↑この論文には次のような背景があるらしく、なんか楽しい

Recognising the importance of CTR estimation for online advertising, management at Bing/adCenter decided to run a competition to entice people across the company to develop the most accurate and scalable CTR predictor. The algorithm described in this publication tied for first place in the first competition and won the subsequent competition based on prediction accuracy. As a consequence, it was chosen to replace Bing’s previous CTR prediction algorithm, a transition that was completed in the summer of 2009.

結構古い(2009年)のですが、まだまだ現役なのもまた良い。



ユーザー&アイテムの特徴量をconcatenateした特徴量ベクトルを {  { \bf x} = (x_i), \ \ x_i \in \{0, 1 \}   }, そしてclick or notを表すターゲットとなる値を{ y \in \{ -1, 1 \} } とする。(特徴量ベクトルは何故か全部カテゴリカル。)

この時クリック確率を表す確率分布を

{ p(y|x):= \Phi \left(  \dfrac{y }{\beta} {\bf  w^T x }  \right)}

{  \Phi(t) := \displaystyle \int^t_{-\infty} \dfrac{1}{\sqrt{2\pi}} \exp\left( -\dfrac{x^2}{2} \right) dx }

で定める。ここで{ {\bf w } }が最適化すべきパラメータで、ベイズ的に扱うために正規分布を事前分布を置き予測分布がclosed formで書けるようにする、が、事後分布は解析的には解けない。それを推定するためのオンライン学習アルゴリズム
www.microsoft.com

ここにあるのでそれを利用している。

そのトリックにより、サンプルが生成される度にオンラインで事前分布のパラメータを更新していく形で最適化をしていく。

§3. CTRのモデリングについて-その2-

前セクションの単純なCTRモデルの問題点↓

Our catalog consists of hundreds of thousands of items with thousands of items added or purged daily. The ex- pected click probability may unduly favor popular items over underexposed or recently added items. This is because un- derexposed items will likely be associated with undertrained model weights that potentially underestimate the true click probability.

新規アイテムの訓練データも欲しい、と言う事でThompson samplingしてランダムに新しいアイテムも晒す。*4

§4. 多様性問題の定式化

CTRモデリングによって、各ユーザーに対してアイテムのスコアリストが得られる。
k個をユーザーに表示するとして、スコア順に並べ替え大きい順にk個取ってしまうと種類が偏る-->多様性をもたせたい。



アイテムの特徴量ベクトルの集合を{ \mathcal{A} := \{ \bf a_1, a_2, a_3,... \}  (a_i \in \mathbb{R}^d ) }とする。

ここで、{ \mathcal{A} }上の劣モジュラ関数 {\rho : 2^\mathcal{A} \rightarrow \mathbb{R}} を次のように定義する:

{  \rho(A,{\bf u}) := \displaystyle \left<{\bf w}_u, \log\left( 1 + \sum_{ {\bf a } \in A}{\bf a }  \right)  \right> + \sum_{ {\bf a } \in A} s({\bf a } ) }

ここで第一項は {\mathbb{R}^d  } 上の標準内積{s} はCTRスコアを出力する関数。また、{ {\bf w}_u }{\mathbb{R}^d  } に埋め込まれたユーザーベクトル(これは次のセクションで。)

そして最適なk個のアイテムリストは次のように与えられる:

{  A^*_k({\bf w}_u)  := \displaystyle \underset{A \in 2^\mathcal{A}, |A|=k}{argmax}   \ \ \rho(A,{\bf w}_u)}

対数を噛ませる事で、一方向だけの成分を伸ばすことで最大化されていくのを防いでいる。これにより多様性が得られる。

§5. ユーザーのモデリング

この論文では、アイテムの特徴量はすべてカテゴリ値になっている事に注意する。
ユーザー {u} のクリック回数を {m_u} 、カテゴリー毎のクリック回数を表すベクトルを { {\bf c_u} }とする。

各クリックはユーザー毎に多項分布に従い、多項分布のパラメータはユーザー毎に、共通のディリクレ分布から生成されているとしてモデリングする。

{ {\bf w}_u \sim Dirichlet ({\bf \alpha} ) }

{ {\bf c}_u \sim Multinomial ({\bf w}_u) }

ここで{ {\bf \alpha} } はディリクレ分布のパラメータ。


事後分布もディリクレ分布であり、期待値{ \hat{{\bf w}}_u}はexplicitに書くことが出来る:

{  \hat{{\bf w}}_u \ = \ \left({\bf c}_u +{\bf \alpha} \right) \| {\bf c}_u +{\bf \alpha} \|^{-1}_1  }

§6. 所感

◯カテゴリ値にこだわる意味がよく分からない。

◯劣モジュラに出てくるユーザーの特徴量ですが、結局LDAっぽい事やっているだけで、もっとPersonalize出来るのでは。

◯ユーザー毎に劣モジュラ最大化を解くの、相当しんどいのでは。

◯劣モジュラにこだわらなくても、アイテムのベクトル化が上手く行っていれば、内積を使って上手くフィルタリングできそうな気がする。(劣モジュラ無知勢なので逃げたい)

*1:4月から推薦システム絡みの機械学習をすることになりそうなので、それ関係の記事が増えるかと思います。

*2:人 工 知 能 32巻1号(2017年1月)の神嶌先生の記事参照

*3:ウィンドウ・ショッピングをECでやる感じのサービス?

*4:この辺は次元圧縮してアイテムをモデリングすればやらなくても良さそう。

ただの微分幾何学徒だった僕がデータサイエンスを何故/どのように勉強したのか

こんにちは。久々の投稿です。

僕のTwitterをフォローしてくれている方はご存知かと思いますが、4月から機械学習エンジニア/データサイエンティスト(見習い)として働く事が決まりました。


良い区切りですので今回はタイトルの通り、ただの純粋数学の学生だった僕がデータサイエンスの勉強を何故/どのようにしてきたのか、についての思い出せる範囲で書こうと思います。


Disclaimer: この記事は基本的に、"What I did" に関する記事であって決して "What you should do" についての記事ではありません。そんな勉強方法おかしいとか、こうすべきだ、みたいなマサカリは一切受け付けませんのでご了承を。ただ、最後に数学科の後輩の皆様に向けてアドバイスと言うか心構えというかそんな感じの事を書いてます。そこはマサカリをTwitterでぶん投げて下さい、と言う感じです。

§1. データサイエンス・機械学習に出会ったきっかけ

ちょうど一年ほど前、とある外資系企業から内定を頂き就職活動を終えました。元々純粋数学で博士号を取るつもりで勉強・研究していたのもあって修論の内容もこの時点でほぼ固まっており、かなり時間に余裕がありました。

そんな時ふと「一度くらい応用数学を勉強してみたい。」と言う衝動に駆られ、ミーハーな僕は数理ファイナンスの勉強をはじめました。するとこれがまたおもしろい。世の中にはこんな高度な数学が応用される世界があるのか、と。当時は確率や統計も数学科の学部の講義程度の知識しかなかったので、結構苦労しました、が、ドハマリしました。

最初に読んだ本は(記憶が正しければ)

ファイナンスのための確率解析 II (連続時間モデル)

ファイナンスのための確率解析 II (連続時間モデル)

ファイナンスのための確率微分方程式―ブラック=ショールズ公式入門

ファイナンスのための確率微分方程式―ブラック=ショールズ公式入門

確率微分方程式

確率微分方程式

共立講座21世紀の数学 (27) 確率微分方程式

共立講座21世紀の数学 (27) 確率微分方程式

Python for Finance: Analyze Big Financial Data

Python for Finance: Analyze Big Financial Data

こんな感じでした。

そんなこんなで数ヶ月数理ファイナンス、特にデリバティヴの価格決定問題を勉強していくうちに一つの疑問が湧いてきました。

「原資産の確率過程のキャリブレーションや推定はどうやって行うのだろう?ここに精度がなければ価格決定の理論もクソもないのでは?」*1

ここで出会ったのが機械学習でした。

で、元々Pythonはある程度馴染みがあったこともあってどんどんハマっていき、今となっては数理ファイナンスには全く興味がないと言った感じです。笑



§2. どのように機械学習/データサイエンスを勉強していったのか

機械学習”と言うワードを知った後、一番最初に勉強したのは神嶌 敏弘先生の

1. 機械学習の Python との出会い — 機械学習の Python との出会い

でした。Pythonに馴染みがあったのでとっつきやすかったです。その後

2. Neural Networks and Deep Learning

初めて深層学習に出会い、そして初めて最適化数学(の初歩の初歩ですが純粋数学ではまず使わない内容)を知りました。*2

この頃のモチベーションはまだ数理ファイナンスであり、超絶ミーハーデータサイエンス芸人な僕は「株価のモデリングをしたい」と思っていました。そして出会ったのが

3. Deep Learning for Multivariate Financial Time Series

です。この論文(?)を読んでいく中で知ったのが(今は廃れてしまった?)制限付きボルツマンマシンによる貪欲学習で、次の2つの論文とノートを眺めながら3.をフルスクラッチで実装する事にしました*3:

4. An Introduction to Restricted Boltzmann Machines

5. LEARNING DEEP GENERATIVE MODELS


こうやって機械学習にハマっていくうちに、段々と数学徒ならではのお気持ち:

「Deep Learningも機械学習も数学テキトウすぎ……論文読んでてイライラする……」

に支配されていきました。それが次の論文リストにつながっています:

mathetake.hatenablog.com

お陰さまで人よりちょっとは深層学習の理論に詳しくなりました。

また、深層学習だけではなく今流行の統計モデリング勢のアプローチも気になりだし、情報幾何学ベイズ統計の勉強をはじめました。その時読んだのは

ベイズ統計の理論と方法

ベイズ統計の理論と方法

Algebraic Geometry and Statistical Learning Theory (Cambridge Monographs on Applied and Computational Mathematics)

Algebraic Geometry and Statistical Learning Theory (Cambridge Monographs on Applied and Computational Mathematics)

Methods of Information Geometry (Tanslations of Mathematical Monographs)

Methods of Information Geometry (Tanslations of Mathematical Monographs)

情報幾何学の基礎 (数理情報科学シリーズ)

情報幾何学の基礎 (数理情報科学シリーズ)

こんな感じでした。




かなり短いですが、ざっとデータサイエンスや機械学習について勉強した流れはこんな感じです。

合間合間でテキトーにSVMに関する論文読んだり、漸近理論の勉強したりしていましたが、基本こんなもんです。

§3. 数学科の皆様へ

正直、真面目に純粋数学をやってきた皆さんにとって機械学習や統計モデリング使う側に*4なるのはかなり楽勝です。根気と時間さえあれば。最初は知らない用語だらけ(僕は今でも知らない単語だらけ)ですが、数学的意味や背景を汲み取るスピードには自信があると思いますのですぐに最前線まで辿り着けるかと思われます。

ただ、こっちの世界に一歩踏み出す前に、

「数学的に厳密じゃない」と言う事に固執しない

と決心してください。こっちの世界で生きていくには数学的厳密さを求めいていたら無理です。不可能。

僕はもう慣れましたが、最初はかなり苦労しました。先ほどの論文リストを作ってしまうほどに固執していました。

数学的厳密さを追い求めるばかり、誰も引用しないような定理を証明する、そんな事になっては本末転倒です。*5

データ分析・機械学習モデルの構築の目的はあくまで応用であって、数学的正しさを証明する実験ではありません。もしこれが嫌であれば一生純粋数学をやっていてください。きっとビジネスの世界では生きていけないでしょう。(と言ってもまだ僕はひよっこですが。)



ただ、数学的厳密さを求める姿勢はもたなくても、数学徒特有の、数式を見て頭の中で捏ねくり回し別のアイデアを思いつく、あの感覚は間違いなく生きます。そして、その感覚は数学科出身でない人たちに対する武器になります。

その数学的センスを武器に、目の前の人の生活や人生に影響を与えるそんな仕事が目の前に、今はそんな時代です。

是非こっちの世界で一緒に頑張りましょう。




物理出身の機械学習マンは大勢いる一方で、数学系の人は少なく、正直悔しくてこんなくだらない文章を書いてしまいましたが、少しでも何か感じ取ってくれたら幸いです。

質問等ありましたら、可能な限りで答えますので Twitter↓まで御連絡ください。

マスタケ (@MATHETAKE) | Twitter


あ、最後に一つだけ


Pythonは絶対に書けるようになってください。


まっっっっっっっったくまとまりのない文章になってしまいました。すみません。

ひよっこデータ芸人の戯言、失礼しました。

*1:この辺はちょっとマサカリ飛んできそうですが、若気の至りです。許してください。

*2:恥ずかしながら当時は勾配法すら知りませんでした。

*3:とても人様に見せられるレベルのものではなかったと思います…。

*4:研究する側になったらそれはまた別の話。測度論的確率論を数学的にきっちりやっている人はざらにいます。

*5:実際DNNのフィッシャー情報行列の特異点を一般の場合に分類しようとして一ヶ月ぐらい溶かしました。あれは時間のムダだった。

ラーメン屋の名前からラーメン屋の食べログスコアを推定出来ると思って頑張ったがやっぱり無理だった話

「うわっこのラーメン屋不味そう」

「変な名前だしやめとこ」


みたいな感覚、誰しもあると思います。(僕だけ…?)



と言うわけで、ハヤシくんはやしくんさん (@tree0_0tree) | Twitter全面協力の元

日本中のラーメン屋のデータを集め

ラーメン屋の名前から食べログスコアを推定に挑戦してみました、と言う報告です。

§1. データセット

今回はスーパーエンジニアハヤシくんがデータを猿のようにかき集めて、

Columnsが
1列目: ラーメン屋の名前
2列目: 食べログスコア
3列目以降 : ラーメン屋の名前をmecab分かち書きした各単語

で構成されるcsvファイルを都道府県毎に用意してくれました。

感謝感激。

§2. w2vのモデル作成

まず、すべてのラーメン屋の名前をramen_name.txtファイルに書き出しました。
そこからコーパスを作りw2vのモデルを作成しました:

from gensim.models import word2vec
data = word2vec.Text8Corpus('ramen_models/ramen_name.txt')
model = word2vec.Word2Vec(data, size=50, min_count=1, window=3)
model.save("ramen_models/w2v_models/ramen_word2vec_size_50_min_1.model")

size=50,100,200
min_count = 1,2,3

で計9個のモデルを作りました。

§3. 前処理

◯ターゲットについて

本来であれば

ラーメン屋の名前の単語系列 --> 食べログスコア(スカラー)

と言う回帰モデルを構築したかったのですが、非自明なモデルが全く作れなかったため、

ターゲットとベクトルは食べログスコア(2.8ぐらい〜4.3ぐらいの区間)を4分割に

離散化して分類タスクにすることにしました。*1

◯特徴量について

w2vのモデルを使って

ラーメン屋の名前 --> 単語数 × w2vのモデルの出力次元

と言う時系列データに変換、そして、最大の単語数を算出し、0ベクトルを追加して系列の長さを整えました。

§4. モデル

前処理によって得られたのは

◯shape = (全国のラーメン屋の数) x (全国のラーメン屋の名前の最大単語数(系列の長さ)) x (w2vモデルの出力次元) のテンソル

◯大きさが全国のラーメン屋の数の、離散化(4分割)された食べログスコアクラスに値を持つベクトル

です。

これを学習させるためにKerasでRNNを構築しました。

より具体的に以下の図のようなモデル*2

f:id:mathetake:20170205205740p:plain

を実装しました。図ではAverageのPoolingレイヤーが入っていますが、結果的にはMaxPoolingしたほうが精度が良かったです。

コードは以下のようです:

model = Sequential()
model.add(SimpleRNN(hidden_neurons, batch_input_shape=(None,length_of_sequences, in_neurons),activation="tanh", return_sequences=True))  
model.add(MaxPooling1D(pool_length=length_of_sequences)) #出力は (None,1,hiden_neuraons)
model.add(Reshape((hidden_neurons,))) #出力を (None,hiden_neuraons)に直す
model.add(Dense(out_neurons,W_regularizer=l2(0.01)))
model.add(Activation("softmax"))
model.compile(loss="categorical_crossentropy", optimizer="RMSprop")
model.fit(x_train, y_train, batch_size=batch_size, nb_epoch=nb_epoch, validation_split=0.05)

§5. 結果

どんなに頑張ってチューニングしたりモデルチェンジしても、学習に使っていない、全データの中の20%(約6000件)に対するテストで

31±0.35%

の精度は超えられませんでした。(§4のモデルが一番精度良かったやつです。)


4値分類タスクなので、ランダムに選ぶ(25%)よりは6%ほど精度がよいので非自明っちゃ非自明ですが……。



おもしろいもの作ってウェブサービスにしましょ〜〜〜〜〜〜〜〜〜〜〜〜〜〜〜!


って意気込んでいた二人の落胆たるや、お分かりいただけますでしょうか。


いいですか皆さん


ラーメン屋さんは名前で選んではいけません。

*1:区間のサンプルの数が均一になるようにpandasのqcutで離散化しました。

*2:LSTM Networks for Sentiment Analysis — DeepLearning 0.1 documentation

"PyMC3" x "DNN" x "ADVI"で学習させて重みを可視化したら綺麗だったと言う話

こんにちは。


先日今話題沸騰中の ライブラリ Edward でVariational InferenceでBayesian DNNを学習させてみたと言う記事を書きました:

mathetake.hatenablog.com


今回の記事は、Edwardの対抗馬(?)の一つであるPyMC3 を使って

DNNの分類器をADVIでミニバッチ学習させるまでの道のりとその結果を簡単にご報告します。


例のごとく、Higgs粒子のデータセット

mathetake.hatenablog.com

を使います。

PyMCの開発者のブログ記事↓をかなり参考にしました。
Bayesian Deep Learning


コード全文はこの記事の下にGistのリンク張っておきますのでそちらを。

§1 準備

まず諸々のモジュールをimportしてから
前回の記事で作成したデータセットを読み込み、加工します。

import pymc3 as pm
from pymc3 import Model, Normal
import pandas as pd
import theano.tensor as T
import theano
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('white')


#ファイルの読み込み
train_df = pd.read_csv("outlier_removed_high_feature_train.csv", index_col=0, header=0)
test_df = pd.read_csv("high_feature_test.csv", header=0, index_col=0, names=['is_Higgs','m_jj', 'm_jjj', 'm_lv', 'm_jlv', 'm_bb', 'm_wbb', 'm_wwbb'])
train_df = train_df.reindex(np.random.permutation(train_df.index)).reset_index(drop=True)

#トレーニングデータの目的変数をまず-1を0に変換してから抽出
train_df.loc[(train_df.is_Higgs == -1 ),"is_Higgs"] = 0
Y_train = train_df.values[0::,0]

#トレーニングデータの説明変数の正規化
train_df_means = train_df.mean().values
train_df_max = train_df.max().values
train_df_min = train_df.min().values
X_train = train_df.apply(lambda x: (x-x.mean())/(x.max()-x.min()), axis=0).fillna(0).values[0::,1::]

#テストデータの目的変数もまず-1を0に変換してから抽出
test_df.loc[(test_df.is_Higgs == -1 ),"is_Higgs"] = 0
Y_test = test_df.values[0::,0]

#テストデータの説明変数を上の正規化定数を使って正規化
X_test = test_df.values[0::,1::]
for i in range(X_test.shape[1]):
  X_test[0::,i] = (X_test[0::,i] - train_df_means[i+1])/(train_df_max[i+1]-train_df_min[i+1] )

あとは、ネットワークの構成で使う定数を定義しておきます。

total_size = len(Y_train) #トレーニングデータの総数
n_feature = X_train.shape[1] #特徴量の数
n_first_layer_hidden = 20 #第一隠れ層のニューロンの数
n_second_layer_hidden = 15 #第二隠れ層のニューロンの数
mini_batch_size = 10000 #ミニバッチサイズ

§2 ミニバッチの作成

次にミニバッチを返すgeneratorを定義し、ミニバッチを作成します。

# Generator that returns mini-batches in each iteration
def create_minibatch(data, mini_batch_size=50,seed=0):
    rng = np.random.RandomState(seed)   
    while True:
        ixs = rng.randint(len(data), size=mini_batch_size)
        yield data[ixs]
minibatches = zip(
    create_minibatch(X_train, mini_batch_size), 
    create_minibatch(Y_train, mini_batch_size),
)

§3 モデルの構成

準備が出来たのでモデルの構成をします。

ですがモデル構築の前に、
そのInputとtargetになるSharedテンソルを作成しておきます。

ann_input = theano.shared(X_train)
ann_output = theano.shared(Y_train)

ただ学習させるだけであれば、sharedテンソルを用意する必要はありません*1
ですが、今回はテストデータに対する予測値も出したいので、
このX_test,Y_testを放り込むためにsharedテンソルを使います。*2


今回はEdwardの時と同じ4層(7-20-15-1)のネットワークを学習させます。*3
まず事前分布を設計してきます。

with pm.Model() as neural_network:
    #事前分布 on weights (標準正規分布で正則化)
    weights_in_1 = Normal('w_in_1', mu=0, sd=1, shape=(n_feature, n_first_layer_hidden) )
    weights_1_2 = Normal('w_1_2', mu=0, sd=1, shape=(n_first_layer_hidden, n_second_layer_hidden))
    weights_2_out = Normal('w_2_out', mu=0, sd=1, shape=(n_second_layer_hidden,))

    #事前分布 on bias (標準正規分布で正則化)
    bias_on_1 = Normal("b_on_1", mu=0, sd=1,shape=(n_first_layer_hidden,))
    bias_on_2 = Normal("b_on_2", mu=0, sd=1,shape=(n_second_layer_hidden,))

標準正規分布を事前分布として正規化の役目を果たしてもらいます。

次に実際にネットワークを構築します。

with neural_network:
    act_1 = T.tanh(T.dot(ann_input, weights_in_1) + bias_on_1)
    act_2 = T.tanh(T.dot(act_1, weights_1_2) + bias_on_2)
    act_out = T.nnet.sigmoid(T.dot(act_2, weights_2_out))

    # Binary classification -> Bernoulli likelihood
    out = pm.Bernoulli('out', act_out, observed=ann_output)

アウトプットはベルヌーイ分布です。
インプットが上で定義したsharedテンソルann_input
Bernoulliobserved 変数にsharedテンソルann_outputが渡されています。

以上でモデル構築完了です。

§3 ミニバッチADVIのセッティング

学習のセッティングに入ります。

with neural_network:
    minibatch_RVs = [out]
    minibatch_tensors = [ann_input, ann_output]

    #minibatch_ADVIで推定
    v_params = pm.variational.advi_minibatch(n=50000, 
        minibatch_tensors=minibatch_tensors, 
        minibatch_RVs=minibatch_RVs,
        minibatches=minibatches, 
        total_size=total_size, 
        learning_rate=1e-2, epsilon=1.0)

まず、minibatch_RVsで、ミニバッチ学習させる確率変数を格納します。
さらにminibatch_tensors上で定めた確率変数で使われているテンソルの中で、ミニバッチに分割するテンソルを格納します。

この2つと、そのミニバッチ化を実現する先ほど作成したジェネレーターminibatches
pm.variational.advi_minibatch に渡すことでADVIのミニバッチ学習ができるわけです。

§4 学習開始

上で構築したモデルを学習させます。
コードはただ1文で終わりです。
このコードにより、v_paramsで定められた学習が実行され、
学習終了後に事後分布から5000個のサンプルが生成され、traceに格納されます。

with neural_network:
    trace = pm.variational.sample_vp(v_params, draws=5000)

こんな感じのプログレスバーが出てくるはずです↓
f:id:mathetake:20170122155031p:plain

§5 学習結果の確認

学習の過程を確認するには以下のコードで行けます。

plt.plot(v_params.elbo_vals)
plt.ylabel('ELBO')
plt.xlabel('iteration')

f:id:mathetake:20170122155123p:plain

なんだか上手い具合に学習してそうな感じがしますね。


次にパラメータのサマリを確認します。

pm.traceplot(trace);

すると次のように重みの分布が可視化されます。
f:id:mathetake:20170122155630p:plain
f:id:mathetake:20170122155637p:plain
各WeightやBiasの中の各成分が色分けされています。
左側の列は密度推定で可視化された分布図です。右側はサンプリングの回数を横軸に取ったパラメータの分布です。
デルタ関数みたいになっているので、綺麗に最適化(or局所最適化)されていると考えられるでしょう。

綺麗!

ちなみにミニバッチサイズ=2,000で、ADVIのパラメータ更新の回数を1/10にした際のサマリー(一部)は
こんな感じでした↓笑

f:id:mathetake:20170122160738p:plain
全然学習できていません。


さて最後になりましたが、テストデータに対する予測を算出しましょう。

まずsharedテンソルにテストデータを放り込みます。

ann_input.set_value(X_test)
ann_output.set_value(Y_test)

次に事後分布からサンプルを生成します。

#(パラメータの)サンプルを事後分布から生成
ppc = pm.sample_ppc(trace, model=neural_network, samples=50)

最後にppcを使ってAccuracyを算出します。

#アウトプットのサンプルを取り出し、サンプルの平均を取る。
#そして0.5より大きければTrueを取るテストデータサイズのベクトルをpredに格納。
pred = ppc['out'].mean(axis=0) > 0.5

#Accuracyの算出
print('Accuracy = {}%'.format((Y_test == pred).mean() * 100))

得られた結果は
f:id:mathetake:20170122160422p:plain

と言う感じです。
頑張ってハイパーパラメータ調節をして通常の最適化手法で学習した際Accuracyはだいたい70%ほどでしたので、
一発でこの結果が出たと言うのはまずまずと言った感じです。

もっとしっかりと調節すれば良い結果が得られるかもしれません。


楽しい!!!楽しいぞベイズ学習!!



今回の記事は以上です。

PyMCは日本語の詳しい情報がほとんど存在しませんが、
慣れれば簡単にDNNを作れますので、ぜひ一度DNNをベイズ学習させてみてはどうでしょうか?


コード全文↓
PyMC3でDNNをADVIで学習 · GitHub

*1:この場合モデルにデータをそのまま放り込めば良いです。

*2:詳しくは開発者のブログを見て下さい。

*3:全然ディープじゃないやんけ!って方。御尤も。

アヒル本のモデルをいくつかPyMCで実装しました. 詰まったところ.

今話題のアヒル本

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)


の後半の方にあるモデルをいくつかPyMC3で実装しました:

github.com


特に一番重要であろうChapter8は全部実装してあります。

間違いや、こうしたらどうですか?みたいなコメントあったらTwitterまで御連絡ください。


この本は”統計モデリングとは”から始まり懇切丁寧にベイズ統計の実践方法を解説してあり、とてもためになります。

データ分析プロセス↓に並んで、Rユーザー以外にもおすすめ出来ます。*1

データ分析プロセス (シリーズ Useful R 2)

データ分析プロセス (シリーズ Useful R 2)



と、言う報告で終わらせたかったのですがPyMC3で実装していて詰まった所をメモがてら残しておこうと思います。


§1. クラス値ベクトルの代入

ahiru_book_pymc/model8-2_pymc.py at master · mathetake/ahiru_book_pymc · GitHub
model 8-2 について考えます。
ここで考えるデータの行は会社員一人一人のサンプルに対応していて、
各サンプルは変数として年齢-給与-会社 と言う変数を持っています。

そしてここでは、
会社(Class)ごとに、給与(Y)年齢(X)を説明変数として線形回帰するものです。

コードを見てみましょう:

df = pd.read_csv("data-salary-2.txt")
X_data = df.values[:,0]
Y_data = df.values[:,1]
Class_data  = df.values[:,2]-1
n_Class = len(df["KID"].unique())


basic_model = Model()

with basic_model:
    a = Normal('a', mu=0, sd=10, shape=n_Class)
    b = Normal('b', mu=0, sd=10, shape=n_Class)
    epsilon = HalfNormal('sigma', sd=1)

    #likelihood 
    mu = a[Class_data] + b[Class_data]*X_data
    Y_obs = Normal('Y_obs', mu=mu, sd=epsilon, observed=Y_data)
    trace = sample(2000)
    summary(trace)

さて、会社毎(Class)に切片と傾きを求める必要があります。つまりn_Class個の切片と傾きがいるわけです。

ですのでまず a(切片) とb(傾き) に対する事前分布に shape=n_Classを渡しています。
ここまでは良いでしょう。

そのあと尤度を計算する際に

a[Class_data]
b[Class_data]

と言う記述があります。Class_data大きさがn_data のnumpyベクトルで、
会社のクラスの値(0~n_Class-1)を取っています。*2

このベクトルを a や b に代入すると、shape = n_data のvariableになり、各成分がクラスに対応した確率変数になっています。

うーーーん。難しい。いや、理解はしてるんですが、なんとなく気持ち悪い。

開発者のブログポストのコメント欄でも、counterintuitiveであると本人が述べています:

The Best Of Both Worlds: Hierarchical Linear Regression in PyMC3





model11-8はもっと複雑で、上述の操作の行列に代入バージョンをやってます。。。。


§2. 変化点検出の収束が異常に遅い


model12-7はコーシー分布を使って変化点検出をするんですが、収束が異常に遅い。

1/4のデータ数でなんとかサンプリングしましたが、収束してるのかどうか怪しい…。

パラメータのサマリ↓
f:id:mathetake:20170122000323p:plain

実際のデータと予測のプロット↓
f:id:mathetake:20170122000331p:plain

どうにかしたいけど………。

しっかりと逆関数法を使って再パラメータ化もしてますが……。


追記:

著者の@berobero11さんからリプをもらいました:

たぶん原因はこれです. 気が向いたらまた実験します.

*1:僕自身Rユーザーではありません。

*2:0~n_Class-1 に直すためにデータを取り出した段階で -1 しています。1~n_Classだとたしかエラーがでました。