Abnormally Distributed

統計解析担当のライフサイエンス研究者 -> データサイエンティスト@コンサル

Googleのデータサイエンスブログ(Misadventures in experiments for growth)

The Unofficial Google Data Science Blogに統計学まわりのなかなか興味深い記事が書いているので、軽く要約してみる。 www.unofficialgoogledatascience.com

大規模な実験はプロダクトの成長に必要不可欠である。ただし、確立したプロダクトではうまくいく標準的な実験手法が成長中のプロダクトに対しては有効でないことがある。 注意深く実験を進めないと、プロダクトの成長を助けるどころか誤った方向に導いてしまう。

Fledgling product

fledglingは「駆け出しの」と言う意味。
established productは確立されたユーザー層がおり、拡大することは好ましいが存続のために必須というわけではない。 一方、fledgling productはまだマーケットを探している段階。

実験の目的

  • To decide on incremental product improvements
    特定の指標に関する仮説検定を行い、意思決定する。established productの成長には有効。 fledgling productでは検出力が足りず小さな変化を同定できないか、小さな改善を熟考するような余裕がない。

  • To do something sizable
    sizableは「かなり大きな」
    fledgling productで行われるような大規模な変更の場合、通常はABテストは適用できない。 しかしながら、長期的なholdback experiment (一部のユーザーのみ変更を保留しておく)をする意義はある。 最終的な影響を事後測定できるため、今後何を試すべきかのヒントが得られる。

  • To roll out a change
    機能変更を展開する前に、問題がないかどうか確認するため。ランダムサンプリングしたユーザーに対して機能変更を適用する。 ランダム化することでselection biasが回避でき、少ない処置群で退化の有無を検出できる。実験インフラの有用な利用法だが、仮説検定のための実験ではない。

ユーザーの想定の間違い

ユーザー全体からランダムサンプリングし、施策の効果を検証した結果、施策Aが一番効果が高かったとする。
ただしプロダクトはまだ成長の途上であり、現段階では特定のニッチな層の利用者が多いとする。 施策Aが有効となったのは、この特定のユーザー層に対して有効だったためであり、他のユーザー層にはむしろ逆効果だったとする。 プロダクトの長期的な成長を考えると、他のユーザー層の利用を拡大する必要があるため施策Aは避けるべきだが、短期的な視点では実施すべきとなってしまう。
現在のユーザー層がターゲットとするユーザー層と異なることを認識する必要がある。

改善案

ユーザーの各セグメントの大きさに基づき、重み付けして効果検証すれば良い。 ユーザーが少数のセグメントに分割できる場合はこれで良いが、多数の区分がある場合は傾向スコアによるユーザーごとの重み付けも考えられる。

その他気をつけるべきこと

  • プロダクトに不満を抱き、周囲に拡散することでプロダクトに悪影響を与えるユーザーも存在する。
  • Early adaptorは一般人口に比べ、テクノロジーに精通している。これらの層に対して有効な施策が、一般ユーザーにとっては逆効果となることもある。
  • 一人で大量に購入するユーザーが存在する。特にプロダクトがまだ小さい場合のABテストはこうしたユーザーの影響を受けやすい。コンバージョンが独立だという仮定にも反しており、二項分布による推定では不確実性を過小評価することがある。

結論

大規模な実験を専門とする人は評価指標を気にかけがちだが、fledgling productを扱う際にはターゲットとするユーザー層にも気を配るべき。
ターゲットユーザーがこのプロダクトを欲しているのか、常に自問しよう。

Tensorflow2.0でのDNNモデル実装方法

1. Keras

一番手軽な方法。

実装の流れ

  1. tf.keras.model.Sequentialなどを使ってモデル構築
  2. model.compileで学習方法を設定 (optimizer, loss, metrics)
  3. model.fitで学習
model = Sequential()
model.add(Dense(...))

optimizer = optimizers.SGD()
model.compile(optimizer=optimizer, 
              loss='binary_crossentropy', 
              metrics=['accuracy'])
model.fit(x_train, y_train, epoch=100, batch_size=10)

# test_pred = model.predict(x_test)
loss, acc = model.evaluate(x_test, y_test)

2. Modelクラス利用

複雑なモデルの実装や細かい設定をしたい場合はこちらを利用。

実装の流れ

  1. tensorflow.keras.models.Modelを継承してモデル定義
  2. 損失関数、最適化手法を定義
  3. 学習ステップを定義(勾配計算とパラメータの更新)
  4. 学習ループを実装
from tensorflow.keras.models import Model
from tensorflow.keras import losses, optimizers, metrics
from sklearn.utils import shuffle

class MLP(Model):
    def __init__(self, hidden_dim, output_dim):
        super().__init__()
        self.l1 = Dense()
        self.l2 = Dense()

    def call(self, x):
        h = self.l1(x)
        y = self.l2(h)
        return y

model = MLP()

criterion = losses.BinaryCrossEntropy()
optimizer = optimizers.SGD()

def compute_loss(t, y):
    return criterion(t, y)

def train_step(x, t):
    with tf.GradientTape() as tape:
        preds = model(x)
        loss = compute_loss(t, preds)
    grads = tape.gradient(loss, model.trainable_variables)
    optimizer.apply_gradients(zip(grads, model.trainable_variables))
    return loss

epochs = 30
batch_size = 100
n_batches = x_train.shape[0] // batch_size

for epoch in range(epochs):
    x_, t_ = shuffle(x_train, t_train)

    for batch in range(n_batches):
        start = batch * batch_size
        end = start + batch_size
        train_step(x_[start:end], t_[start:end])

モデル定義において、callで順伝播を実装する。これにより、model(x)のように書くことができる。 Pythonでは下記のようにcall()を実装すると、model.forward()の代わりに、model(x)と記述できるようになる。

class Model(object):
    
    def __call__(self, x):
        return self.forward(x)

3. Estimator API

高レベルなAPI
学習、予測、評価、モデルのエクスポートの機能を組み込んでいる。
作成済みのEstimatorを用いることも、新しく作ることも可能。

強み

tf.kerasではまだ開発中の下記の機能が利用可能。 * Parameter server based training * Full TFX integration.

また、並列処理を安全に行うための機能がある。 あまりドキュメントが充実しておらず、TF2では下火なのかもしれない。

使い方

入力データをtf.data.Datasetに変換し、ストリーミング形式で入力パイプラインに与えられるようにするためのinput_functionが必要。

def make_input_fn(data_df, label_df, num_epochs=10, shuffle=True, batch_size=32):
  def input_function():
    ds = tf.data.Dataset.from_tensor_slices((dict(data_df), label_df))
    if shuffle:
      ds = ds.shuffle(1000)
    ds = ds.batch(batch_size).repeat(num_epochs)
    return ds
  return input_function

train_input_fn = make_input_fn(dftrain, y_train)
eval_input_fn = make_input_fn(dfeval, y_eval, num_epochs=1, shuffle=False)


#define feature columns
CATEGORICAL_COLUMNS = ['sex', 'n_siblings_spouses', 'parch', 'class', 'deck',
                       'embark_town', 'alone']
NUMERIC_COLUMNS = ['age', 'fare']

feature_columns = []
for feature_name in CATEGORICAL_COLUMNS:
  vocabulary = dftrain[feature_name].unique()
  feature_columns.append(tf.feature_column.categorical_column_with_vocabulary_list(feature_name, vocabulary))

for feature_name in NUMERIC_COLUMNS:
  feature_columns.append(tf.feature_column.numeric_column(feature_name, dtype=tf.float32))


# training and evaluation
linear_est = tf.estimator.LinearClassifier(feature_columns=feature_columns)
linear_est.train(train_input_fn)
result = linear_est.evaluate(eval_input_fn)

参考リンク

[1708.02637] TensorFlow Estimators: Managing Simplicity vs. Flexibility in High-Level Machine Learning Frameworks

Estimators  |  TensorFlow Core

Premade Estimators  |  TensorFlow Core

変分推論

前回紹介したEMアルゴリズム最尤推定の枠組みになる。 最尤推定では尤度関数の特異性に起因する問題がある。

そこで、今回はベイズ的な枠組みでモデルの推定を行う。ほぼPRMLの要約になっている。

EMアルゴリズム再訪

EMアルゴリズムの目的は対数尤度関数の最適化である。
対数尤度関数は下記のように分解できる。

 \ln p(\boldsymbol{X}|\theta) = \mathcal{L}(q|\theta) + {\rm KL}(q||p)

ただし、潜在変数の分布としてq(\boldsymbol{Z})を導入した。 ここで、KLダイバージェンスは常に非負なので、 \mathcal{L}(q|\theta)は対数尤度関数の下界となる。

そのため、対数尤度の最大化は、 \mathcal{L}(q|\theta)と同じことになる。

Eステップでは、パラメータ\thetaを固定した上で、下界をq(\boldsymbol{Z})について最大化する。 Mステップではq(\boldsymbol{Z})を固定した上で、下界をパラメータ\thetaについて最大化する。

変分推論

パラメータ、潜在変数を全てまとめてZと表記すると、EMアルゴリズムの場合と同様に、周辺分布は下記のように分解できる。

 \ln p(\boldsymbol{X}) = \mathcal{L}(q) + {\rm KL}(q||p)

ここで、下界 \mathcal{L}(q)を分布q(\boldsymbol{Z}) について最大化することは、KLダイバージェンスの最小化と同値になる。 q(\boldsymbol{Z}) = p(\boldsymbol{Z|X})の場合にKLダイバージェンスは0になるが、真の事後分布を求めることはできないとする。

そこで、分布 q(\boldsymbol{Z})のクラスを制限した上でKLダイバージェンスの最小化を考える。

近似分布のクラスの制限としては、下記の平均場近似がよく用いられる。 分布qがいくつかの因子に分解されると仮定している。

 \displaystyle q(\boldsymbol{Z}) = \prod_{i=1}^{M} q_i(\boldsymbol{Z_i})

ここで { q _ { i \neq j } }を固定した上で、下界 \mathcal{L}(q) q_j(\boldsymbol{Z_j})について最大化する。 最適解は下記の通り。

 \ln q_j^*(\boldsymbol{Z_j}) = \mathbb{E} _{i \neq j} [\ln p(\boldsymbol{X, Z})] + {\rm const}

観測データと潜在変数の同時分布の対数について、 i \neq jである他の因子全てについての期待値をとったものとなる。

全ての因子を適当に初期化した上で、1つ1つの因子を上式に従って改良していくことで、下界を最大化するq(\boldsymbol{Z})を取得することができる。

LDAによるトピックモデル ②gensimによる実装

前回の記事ではLDAの概要や関連手法、確率モデルについて書いた。
今回はPythonのgensimというライブラリを用いて、LDAを実践してみる。
その前に、前回触れていなかったトピックモデルの評価方法について説明する。

評価指標

LDAは教師なしモデルであり、精度等の指標で評価することはできない。
一般的には、PerplexityやCoherenceといった指標が用いられる。

Perplexity

予測性能に関する指標。

定義は下記の通りで、単語ごとの尤度の幾何平均の逆数である。Mは(学習データの)文書数、Ndは文書dの単語数、wdは文書dにおける各単語の出現回数を表す。
教師なし学習では、テストデータの尤度が大きくなるほど、つまりpeplexityが小さくなるほど好ましい。

 \displaystyle perplexity(D _ {test}) = \exp \left(- \frac{\sum _ {d=1} ^ {M} \log p(w _ d)} {\sum _ {d=1} ^ {M} N _ d} \right)

直感的な理解についてはこちらを参照。
トピックモデルの評価指標 Perplexity とは何なのか?

要はPeplexityは選択肢の数を表している。
ある文書の1単語を決める際、何も仮定がない状態ではvocabulary中の全単語が候補になるが、モデルの元では候補となる単語がもっと絞り込めるようになる。この絞り込んだあとの候補数を表すのがPerplexityである。

Coherence

トピックの品質の良さ、人間にとっての解釈しやすさを表す指標。様々な計算方法が提唱されている。 例として、下記に2つのトピックとそのトピックにおいて出現頻度が高い単語を示す。

Topic1: guitar album track bass drum vocal
Topic2: know call name several refer hunter

Topic1は音楽に関するトピックだと明らかだが、Topic2は何を表しているのかはっきりわからない。この場合、Topic1の方がcoherenceが高くなるべきである。

coherenceが高いのは各トピックに含まれる単語の類似度が高い場合である。そこで、coherenceの算出方法としては、各トピックの上位N個の単語を抽出し、その中の各単語ペアについて何らかの類似度を計算して全ペアでの平均値を求めることが多い。
類似度の算出方法については、いくつかの方法が提唱されている。

coherenceの評価基準としては、トピックの良し悪しについて人がつけた正解ラベルとどの程度相関があるかになる。

各手法について、詳しくはこちらを参照。
トピックモデルの評価指標 Coherence 研究まとめ #トピ本

"Exploring the Space of Topic Coherence Measures" M. Röder, A Both, A. Hinneburg, 2015
https://dl.acm.org/doi/10.1145/2684822.2685324
gensimのCoherenceModelはこちらの論文に基づいている。

gensimによる実装

Kaggleのデータセットを用いる。
約20万件のニュース記事の見出しと短い説明文を含む。また、各記事のカテゴリーも記載されており、抽出されたトピックとの関連を確認することもできそう。

コードはKaggleのNotebookとして動かせるようになっている。
LDA_gensim | Kaggle

from time import time
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from IPython.core.display import display
pd.set_option('display.max_colwidth', 200)
%matplotlib inline

# load data
df = pd.read_json('../input/news-category-dataset/News_Category_Dataset_v2.json', lines=True)
df.sample(3)
category headline authors link short_description date
134770 HOME & LIVING Exploiting Familial Ideals for a Car Alissa Fleck, Contributor\nAlissa Fleck is a reporter based in New York City. https://www.huffingtonpost.com/entry/exploiting-familial-ideal-for-car-commercials_us_5b9dc09fe4b03a1dcc8c71ce Is anyone else more than a little crushed when images of gooey, familial perfection splash dazzlingly across the screen to a backdrop of dramatic music swelling to an emotional climax, only to hav... 2014/1/2
5634 POLITICS Why You Shouldn't Freak Out Even If The Stock Market Does Zach Carter https://www.huffingtonpost.com/entry/stock-market-dont-freak-out_us_5a846bc5e4b0774f31d15aaf Everything is fine ... except your paycheck. 2018/2/14
944 ENTERTAINMENT R. Kelly Accusers Say They Had To Ask Permission To Use The Bathroom Sara Boboltz https://www.huffingtonpost.com/entry/r-kelly-accusers-bathroom-permission_us_5af1e468e4b0c4f19327a1f3 And when they did, the women said he demanded they call him "Daddy." 2018/5/8

次に、テキストの前処理をする。
headlineとshort descriptionをつなげたデータを用いる。

from gensim.parsing.preprocessing import preprocess_string

combined_text = df['headline'] + df['short_description']
processed_text = combined_text.map(preprocess_string)

gensimのpreprocess_stringを使うと、テキストデータの一般的な前処理を一通り適用できる。 デフォルトでは下記7つの処理を実行する。

  1. strip_tags(): HTMLタグ等の除去
  2. strip_punctuation(): 句読点の除去
  3. strip_multiple_whitespaces(): 複数の空白の除去
  4. strip_numeric(): 数字の除去
  5. remove_stopwords(): ストップワード(頻出単語)の除去
  6. strip_short(): 短い単語の除去(デフォルトでは3文字未満)
  7. stem_text(): 小文字にして、語幹のみ抽出(porter-stemmerの適用)

実行前
Hugh Grant Marries For The First Time At Age 57The actor and his longtime girlfriend Anna Eberstein tied the knot in a civil ceremony.

実行後
[hugh, grant, marri, time, ag, actor, longtim, girlfriend, anna, eberstein, ti, knot, civil, ceremoni]

Porter stemmerは古典的なステマーだが、age→ag、tied→tiなど元の単語が分かりにくくなるものも多く、正直微妙かもしれない。

次に、各単語に番号を割り振ったディクショナリーを作成する。この際、極端に出現頻度が少ない単語や多い単語は除外する(filter_extremes)。 さらに、ディクショナリーを用いて前処理済みのデータをbag-of-words表現に変換する。

from gensim.corpora.dictionary import Dictionary

# Create a corpus from a list of texts
dictionary = Dictionary(processed_text)
dictionary.filter_extremes(no_below=10, no_above=0.7, keep_n=100000)

corpus = [dictionary.doc2bow(text) for text in processed_text]
corpus[10]
[(7, 1),
 (8, 1),
 (101, 1),
 (102, 1),
 (103, 1),
 (104, 1),
 (105, 1),
 (106, 1),
 (107, 1),
 (108, 1),
 (109, 1),
 (110, 1),
 (111, 1),
 (112, 1)]

このように、文章は(単語id, 出現回数)のリストで表現されている。

from gensim.models import LdaModel, LdaMulticore
lda_model = LdaMulticore(corpus, num_topics=20, id2word=dictionary, random_state=42, workers=3)

print(lda_model.log_perplexity(corpus))
print()

# Get the most significant topics 
for idx, topic in lda_model.print_topics(-1): # -1 for all topics
    print('Topic: {} Word: {}'.format(idx, topic))

log-perprexityは-9.15となった。 また、各トピックにおける単語の出現確率を高い順に出力している。

-9.149718461222365

Topic: 0 Word: 0.013*"tip" + 0.013*"food" + 0.013*"cancer" + 0.011*"life" + 0.010*"sleep" + 0.007*"help" + 0.007*"studi" + 0.006*"research" + 0.006*"need" + 0.006*"risk"
Topic: 1 Word: 0.030*"dai" + 0.024*"parent" + 0.018*"children" + 0.018*"need" + 0.012*"time" + 0.011*"know" + 0.009*"help" + 0.009*"kid" + 0.009*"child" + 0.009*"learn"
Topic: 2 Word: 0.021*"marriag" + 0.011*"anim" + 0.009*"gai" + 0.008*"appl" + 0.007*"kid" + 0.006*"famili" + 0.006*"time" + 0.006*"sex" + 0.006*"peopl" + 0.005*"coupl"
Topic: 3 Word: 0.027*"new" + 0.018*"york" + 0.016*"citi" + 0.015*"yoga" + 0.013*"exercis" + 0.010*"summer" + 0.008*"practic" + 0.006*"spiritu" + 0.005*"chees" + 0.005*"want"
Topic: 4 Word: 0.026*"fashion" + 0.025*"photo" + 0.017*"year" + 0.014*"video" + 0.011*"old" + 0.010*"style" + 0.009*"men" + 0.008*"school" + 0.008*"game" + 0.007*"week"
Topic: 5 Word: 0.027*"year" + 0.011*"babi" + 0.011*"best" + 0.010*"dai" + 0.010*"weight" + 0.010*"mother" + 0.010*"old" + 0.010*"mom" + 0.009*"hotel" + 0.009*"spring"
Topic: 6 Word: 0.013*"peopl" + 0.008*"kill" + 0.008*"polic" + 0.007*"cocktail" + 0.006*"pari" + 0.006*"china" + 0.006*"year" + 0.006*"video" + 0.005*"said" + 0.005*"human"
Topic: 7 Word: 0.023*"obama" + 0.013*"recip" + 0.010*"presid" + 0.009*"flavor" + 0.009*"dish" + 0.008*"trump" + 0.007*"grill" + 0.007*"late" + 0.007*"jimmi" + 0.006*"question"
Topic: 8 Word: 0.018*"red" + 0.015*"morn" + 0.011*"decor" + 0.009*"carpet" + 0.009*"ic" + 0.009*"welcom" + 0.008*"new" + 0.007*"fruit" + 0.007*"cream" + 0.006*"video"
Topic: 9 Word: 0.025*"love" + 0.015*"divorc" + 0.015*"life" + 0.012*"live" + 0.011*"photo" + 0.010*"time" + 0.010*"dai" + 0.010*"thing" + 0.008*"wai" + 0.008*"like"
Topic: 10 Word: 0.039*"photo" + 0.024*"wed" + 0.016*"check" + 0.013*"twitter" + 0.012*"huffpost" + 0.011*"look" + 0.011*"week" + 0.011*"facebook" + 0.011*"style" + 0.010*"video"
Topic: 11 Word: 0.014*"new" + 0.014*"super" + 0.012*"bowl" + 0.010*"workout" + 0.009*"resolut" + 0.008*"state" + 0.007*"america" + 0.006*"media" + 0.006*"polit" + 0.005*"trump"
Topic: 12 Word: 0.022*"studi" + 0.014*"new" + 0.010*"tumblr" + 0.010*"diseas" + 0.008*"peopl" + 0.007*"kate" + 0.007*"heart" + 0.007*"come" + 0.006*"time" + 0.006*"american"
Topic: 13 Word: 0.017*"new" + 0.011*"dai" + 0.010*"like" + 0.010*"video" + 0.008*"dad" + 0.008*"look" + 0.007*"photo" + 0.007*"eat" + 0.007*"gift" + 0.007*"year"
Topic: 14 Word: 0.030*"health" + 0.022*"care" + 0.011*"tax" + 0.008*"patient" + 0.008*"parti" + 0.007*"work" + 0.006*"want" + 0.006*"plan" + 0.006*"help" + 0.006*"republican"
Topic: 15 Word: 0.013*"click" + 0.013*"video" + 0.009*"clinton" + 0.009*"com" + 0.008*"star" + 0.008*"new" + 0.008*"photo" + 0.008*"hillari" + 0.007*"luxuri" + 0.007*"dog"
Topic: 16 Word: 0.014*"world" + 0.012*"healthi" + 0.012*"recip" + 0.011*"food" + 0.010*"wai" + 0.008*"chang" + 0.007*"time" + 0.006*"eat" + 0.006*"wine" + 0.006*"dii"
Topic: 17 Word: 0.031*"photo" + 0.011*"dress" + 0.010*"bodi" + 0.010*"look" + 0.010*"art" + 0.008*"like" + 0.008*"design" + 0.007*"diet" + 0.007*"green" + 0.006*"wai"
Topic: 18 Word: 0.014*"medit" + 0.011*"vintag" + 0.008*"bank" + 0.008*"san" + 0.007*"new" + 0.006*"ebai" + 0.006*"home" + 0.005*"googl" + 0.005*"sandi" + 0.005*"gun"
Topic: 19 Word: 0.016*"pinterest" + 0.015*"new" + 0.009*"water" + 0.007*"loss" + 0.007*"fit" + 0.007*"american" + 0.006*"park" + 0.006*"wall" + 0.006*"drink" + 0.005*"evolut"

WordCloudでトピックを可視化する。

from wordcloud import WordCloud
from PIL import Image
 
fig, axs = plt.subplots(ncols=4, nrows=int(lda_model.num_topics/4), figsize=(15,10))
axs = axs.flatten()
 
def color_func(word, font_size, position, orientation, random_state, font_path):
    return 'darkturquoise'

 
for i, t in enumerate(range(lda_model.num_topics)):
 
    x = dict(lda_model.show_topic(t, 30))
    im = WordCloud(
        background_color='white',
        color_func=color_func,
        random_state=0
    ).generate_from_frequencies(x)
    axs[i].imshow(im)
    axs[i].axis('off')
    axs[i].set_title('Topic '+ str(t))
        
plt.tight_layout()
plt.savefig('/kaggle/working/wordcloud.png')

f:id:kibya:20200517145648p:plain

ある程度は解釈できそうなトピックが見られる。 例えばTopic0はcancer, food, sleep, research, riskなどから、生活習慣と疾病リスクに関する研究結果の記事などが想像される。
またTopic1はparent, children, kid, learnなどから子供や教育に関するトピックだと思われる。

最後に、coherenceを算出する。top_topicsメソッドを用いると、各トピックのcoherenceと出現確率の高い単語を、coherenceの高い順に表示できる。
デフォルトでは、coherenceはu-massと呼ばれる手法で算出する。

topics = lda_model.top_topics(corpus, topn=10)

print('coherence: top words')
for topic in topics:
    print('{:.3f}: {}'.format(topic[1], ' '.join([i[1] for i in topic[0]])))
coherence: top words
-2.436: photo wed check twitter huffpost look week facebook style video
-2.708: love divorc life live photo time dai thing wai like
-2.939: dai parent children need time know help kid child learn
-3.486: new dai like video dad look photo eat gift year
-3.514: tip food cancer life sleep help studi research need risk
-3.548: health care tax patient parti work want plan help republican
-3.575: fashion photo year video old style men school game week
-3.992: photo dress bodi look art like design diet green wai
-4.208: world healthi recip food wai chang time eat wine dii
-4.223: year babi best dai weight mother old mom hotel spring
-4.578: marriag anim gai appl kid famili time sex peopl coupl
-4.654: studi new tumblr diseas peopl kate heart come time american
-4.664: new super bowl workout resolut state america media polit trump
-5.238: peopl kill polic cocktail pari china year video said human
-5.576: click video clinton com star new photo hillari luxuri dog
-6.067: red morn decor carpet ic welcom new fruit cream video
-6.260: new york citi yoga exercis summer practic spiritu chees want
-6.344: obama recip presid flavor dish trump grill late jimmi question
-6.711: pinterest new water loss fit american park wall drink evolut
-8.871: medit vintag bank san new ebai home googl sandi gun

coherenceが上位のトピックは確かに関連する単語が多いが、下位のトピックは一見して無関係の単語が多くなっており、たしかに人間の感覚と一致しているようだ。

LDAによるトピックモデル ①概要と確率モデル

トピックモデルはテキストデータの潜在的意味の解析に用いられる統計モデルの総称である。 ここでは一番有名なモデルである、Latent Dirichlet Allocation (LDA)について説明する。

LDAとは

文書データをLDAで処理することで、文書のトピックの分類や、トピックごとの単語の出現頻度を知ることができる。

LDAでは、単語の羅列である文章の背後に潜在的なトピック(スポーツ、政治、教育など)があると考え、トピックに基づいて文書中の単語が生成されると考える。 例えばスポーツに関する文章なら「試合」「得点」「選手」などの単語が多く、政治に関する文章なら「議会」「法案」「選挙」などの単語が多いと考えられる。 1つ注意点として、LDAは教師なし学習なので、トピックを教師データとして与えて学習する必要はない。

また、LDAでは文章における単語の順序を考慮せず、単語の出現回数のみを考える(Bag of words)。

トピックモデル とLDAの歴史

トピックモデルの起源は1990年に発表されたLatent Semantic Indexing (LSI, LSA, 潜在意味解析)にさかのぼる。
LSIでは共起行列(文書-単語行列)の特異値分解により、トピックを表す潜在変数を取得する手法。例えばK個のトピックを得たい場合は、特異値が大きいものからK個を選択すれば良い(行列の低ランク近似)。
ただし、LSIは確率モデルではなく、文書データの生成プロセスを直接的にモデリングする手法ではない。

1999年にHoffmanらが提唱したProbabilistic Latent Semantic Analysis (PLSA) は確率モデルである。
PLSAは混合モデルであり、文書ごとのトピックの確率 p(z \mid d)を混合比率、トピックごとの単語の生成確率 p(w \mid zを混合する確率分布としている。 ただし、新規の文書dについて、 p(z \sim d)を生成することができないため、完全な生成モデルとは言えない。つまり学習に用いた文書以外にはうまく適用できる保証がない。

f:id:kibya:20200516153355p:plain
PLSA

一方、LDAは完全な生成モデルであり、2003年、David Bleiらの論文で発表された。 その後、LDAについては教師ありモデル、単語の順番を考慮したモデルなどへの拡張や、最適なトピック数の推定、変分ベイズによる効率的な推論方法など様々な研究の対象となっている。

確率モデル

まず、全体のイメージを下図で確認して欲しい。

f:id:kibya:20200506141706p:plain
LDAのイメージ

次に、具体的な確率モデルを説明する。 まず、トピックの総数をKとする(この数値は自分で設定する必要がある)。 文書はK種類のトピックがトピック比率 (topic proportion) \theta_d)と呼ばれる重みをつけて混合されたものと考える。
例えば、国会で新しい教育制度の法案が可決されたことに関するニュース記事だとすると、トピック「政治」の割合0.7、トピック「教育」の割合0.2、その他のトピックの割合0.1のようになるだろう。

 \theta _ d)はK次元のベクトルであり、各成分の和は1になる。
このような値を生成する確率分布としては、ディリクレ分布(ベータ分布の多次元版)が用いられる。

 p(\theta _ d) = {\rm Dirichlet} (\theta _ d \mid \alpha)

文書中の各単語はK種類のトピックのいずれかから生成されると考える。 まず、文書dにおけるi番目の単語 w _ {d,i}がどのトピックから生成されたかを示すトピック割り当て (topic assignment) z _ {d, i}を導入する。  z _ {d, i}はK次元のベクトルであり、いずれか1つの要素が1、他の要素は0である。

このような値を生成する確率分布としては、カテゴリ分布(ベルヌーイ分布の多次元版)が用いられる。

 p(z _ {d, i} \mid \theta _ d) = {\rm Categorical} (z _ {d, i} \mid \theta _ d)

次に、単語 w _ {d,i} z_ {d, i}で指定されたトピックから生成される。 ただし、単語 w _ {d,i}はいわゆるOne-of-K表現したベクトルである。すなわち、語彙の総数をVとすると  w _ {d,i}は長さVのベクトルであり、該当する単語のインデックスのみ1、他は0となっている。

LDAでは、単語の出現確率がトピックにより異なると仮定しており、トピックkにおける単語の出現比率はV次元ベクトル \beta _ kで表される。  \beta _ kの各成分は0から1の範囲の値であり、成分の和は1になる。
そこで、 w _ {d,i}もカテゴリ分布により生成されると考えることができる。この際、 z_ {d, i}のk番目の成分が1のとき、単語  w _ {d,i}はトピックkから生成されていることに注意する。

 p(w _ {d,i} \mid z _ {d, i}, \boldsymbol{\beta}) = \prod _ {k=1} ^{K} {\rm Categorical} (z _ {d, i} \mid \beta _ k) ^ {z _ {d, i, k}}

ただし、 \boldsymbol{\beta} = \{ \beta_1, \beta_2, \dots, \beta_k \} としている。

また、 \beta _ kの事前分布としては、ディリクレ分布を用いる。
 p(\beta _ k) = {\rm Dirichlet} (\beta _ k \mid \eta)

グラフィカルモデルで示すと、下記の通り。

f:id:kibya:20200506092410p:plain
LDAのグラフィカルモデル

実践

scikit-learnのLatentDirichletAllocationやgensimのLdaModelを使うと簡単に実装することができる。
長くなってしまうので、この記事はここまで。 Tensorflow Probabilityを使ったベイズ推論について、別の記事で説明する。

参考文献

  1. “Latent Dirichlet Allocation” D. Blei, A. Ng, M. Jordan, 2003 (http://www.jmlr.org/papers/v3/blei03a)
    LDAを提唱した元論文。

  2. Chapter5でLDAの変分ベイズによる推論方法を解説している。

  3. "Indexing by latent semantic analysis" S. Deerwester, S. Dumais, T. Landauer, G. Furnas, and R. Harshman, 1990 http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.108.8490 LSI (LSA)の元論文。

  4. "Probabilistic latent semantic indexing" T. Hoffman, 1999(https://dl.acm.org/doi/10.1145/312624.312649) PLSIの元論文。

  5. PLSAの仕組みとその有用性についての記事 www2.deloitte.com

Optunaの概要

仕事でOptunaを使う機会があったので、論文を軽く読んでまとめてみる。

arxiv.org

Optunaとは

OptunaはPFNが開発したハイパーパラメータのチューニングツールであり、近年、勾配ブースティングやニューラルネットワークのチューニングによく使われている。
既存のツールに対するOptunaの特徴として、下記3点が挙げられている。

  1. Define-by-runで動的にパラメータ探索空間を設定可能
  2. 効率的なsampling, pruningアルゴリズム
  3. 設定が簡単で、軽量な実験から分散処理による大規模な計算まで対応可能

1. Define-by-run

define-by-runは深層学習フレームワークの分類に使われる用語だが、筆者らはパラメータチューニングツールにも同様の概念を提唱している。
従来のツールはdefine-and-run、すなわちパラメータの探索範囲を最初に設定した上で、最適なパラメータの探索を行うものだった。

一方、define-by-runでは動的に探索範囲を設定することができる。これは例えばニューラルネットワークのレイヤー数とユニット数を同時にチューニングする際などに役に立つ。

2. 効率的なsampling, pruningアルゴリズム

パラメータ最適化の効率は探索戦略(探索すべきパラメータを決める)と性能評価戦略(学習曲線から現在探索中のパラメータの良さを判断し、捨てるかどうか判断する)により決まる。

論文には探索(サンプリング)アルゴリズムの詳細は記載されていないが、TPE (tree-based Parzen estimator)*1, CMA-ES(Covariance Matrix Adaptation - Evolution Strategy)*2, GP-BO(Bayesian optimization) *3などを用いているとのこと。 サンプリング手法にはパラメータの相関を利用するrelational samplingと独立にサンプルするindependent samplingがある。TPEはindependent sampling、CMA-ES, GP-BOはrelational samplingである。

また、性能評価については、見込みの悪い試行を途中で終わらせること(pruning、枝刈り)が重要である。 OptunaではAsynchronous Successive Halving (ASHA)*4と呼ばれるアルゴリズムを採用している。ASHAでは試行中のパラメータセットの暫定的な順位に基づき、分散処理の個々のworkerが非同期的にpruningを行う手法とのこと。

使い方

Optunaではパラメータの最適化をobjective functionの最小化・最大化の問題と捉える。
最適化のプロセスをstudy、objective functionの1回の評価をtrialと呼ぶ。 パラメータの探索空間はtrialオブジェクトのメソッドにより動的に構築される。

以下のコードでは(モデルのチューニングではないが)、簡単な例として2次関数の最小化を行っている。

import optuna

def objective(trial):
    x = trial.suggest_uniform('x', -10, 10)
    return (x - 2) ** 2

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

study.best_params  # E.g. {'x': 2.002108042}}



以下はより実践的な、TensorFlowでニューラルネットワークのパラメータチューニングを行う場合のコード抜粋。 完全なコードは下記リンク先参照。 https://colab.research.google.com/drive/1MQp5rn6Dxhdv4r9PuJOyGf1sdsS2PHF1

ここでは、隠れ層の数、ユニット数、weight decayといったモデルのパラメータに加え、optimizerの種類とそのパラメータも最適化している。 従来のツールでこうしたチューニングを行うためには、最初にかなり複雑なパラメータ探索範囲を記述する必要があるが、Optunaのdefine-by-runの思想により、簡潔に実装することができる。

def create_model(trial):
    # We optimize the numbers of layers, their units and weight decay parameter.
    n_layers = trial.suggest_int("n_layers", 1, 3)
    weight_decay = trial.suggest_loguniform("weight_decay", 1e-10, 1e-3)
    model = tf.keras.Sequential()
    model.add(tf.keras.layers.Flatten())
    for i in range(n_layers):
        num_hidden = int(trial.suggest_loguniform("n_units_l{}".format(i), 4, 128))
        model.add(
            tf.keras.layers.Dense(
                num_hidden,
                activation="relu",
                kernel_regularizer=tf.keras.regularizers.l2(weight_decay),
            )
        )
    model.add(
        tf.keras.layers.Dense(CLASSES, kernel_regularizer=tf.keras.regularizers.l2(weight_decay))
    )
    return model


def create_optimizer(trial):
    # We optimize the choice of optimizers as well as their parameters.
    kwargs = {}
    optimizer_options = ["RMSprop", "Adam", "SGD"]
    optimizer_selected = trial.suggest_categorical("optimizer", optimizer_options)
    if optimizer_selected == "RMSprop":
        kwargs["learning_rate"] = trial.suggest_loguniform("rmsprop_learning_rate", 1e-5, 1e-1)
        kwargs["decay"] = trial.suggest_uniform("rmsprop_decay", 0.85, 0.99)
        kwargs["momentum"] = trial.suggest_loguniform("rmsprop_momentum", 1e-5, 1e-1)
    elif optimizer_selected == "Adam":
        kwargs["learning_rate"] = trial.suggest_loguniform("adam_learning_rate", 1e-5, 1e-1)
    elif optimizer_selected == "SGD":
        kwargs["learning_rate"] = trial.suggest_loguniform("sgd_opt_learning_rate", 1e-5, 1e-1)
        kwargs["momentum"] = trial.suggest_loguniform("sgd_opt_momentum", 1e-5, 1e-1)

    optimizer = getattr(tf.optimizers, optimizer_selected)(**kwargs)
    return optimizer


def objective(trial):
    # Get MNIST data.
    train_ds, valid_ds = get_mnist()

    # Build model and optimizer.
    model = create_model(trial)
    optimizer = create_optimizer(trial)

    # Training and validating cycle.
    with tf.device("/cpu:0"):
        for _ in range(EPOCHS):
            learn(model, optimizer, train_ds, "train")

        accuracy = learn(model, optimizer, valid_ds, "eval")

    # Return last validation accuracy.
    return accuracy.result()

とりあえず、特徴量を手当たり次第に作ってみてLightGBMの投入し、Optunaでハイパーパラメータ調整する、なんて作業は 誰か他の人にやってもらいたい。。

説明変数のある時系列予測(状態空間モデル)

概要

時系列データは実務で取り扱う機会が多い。その際、予測対象の時系列データ(目的変数)と似たような変動をする他の時系列データ(説明変数)を手がかりに予測したい場合がある。

説明変数を用いて目的変数を予測する最も単純な方法として、線形回帰がある。ただ、線形回帰では時系列データの特徴を捉えきれないという問題があり、高い精度は期待できない。

この記事では時系列分析の1手法である状態空間モデルによる予測を行い、線形回帰に比べて何が優れているかを説明する。

データ

ここでは、下図のように変動するX, Yという2つの変数を考える。 Yが予測対象の変数、Xは説明変数である。 Xを利用して、Yをリアルタイムに予測したい。つまり、t時点のYを予測する際には、t-1時点までのYの値とt時点までのXの値が利用可能とする。

f:id:kibya:20200104224855p:plain

XとYの散布図は下図のようになる。

f:id:kibya:20200104225925p:plain

XとYにははっきりと相関があり、XはYの予測に役立ちそうである。

線形回帰

ここで単純な予測モデル、すなわち線形回帰することを考える。
これはt-1時点までのX, Yを用いて回帰直線を作成し(切片αと傾きβを求め)、t時点のXの値を代入することでYの予測値を算出する。

 {
\begin{eqnarray} 
y_t &=& \alpha+ \beta x_t + \varepsilon_t \\
\varepsilon_t &\sim& N(0, \sigma^2_\varepsilon) \\
\end{eqnarray} 
}

下図は、t=10までのデータ(グレー)を用いて、回帰直線を作成し、t=11のYを予測している。

f:id:kibya:20200104230756p:plain

全時点の予測結果をプロットすると下図の通り。

f:id:kibya:20200104231634p:plain

そこまで悪い予測ではないが、途中に連続して予測を大きく外してしまうケースが見られる。一期先予測誤差(MAE)は2.38程度となった。 これは線形回帰では、過去のデータ点をすべて同等に扱っており、直近のデータは類似しているというデータからはっきり読み取れる傾向を無視しているからである。
この直近のデータは類似している、という性質こそ時系列データの特徴であり、予測の際に考慮すべき情報なのだ。

状態空間モデル

時系列予測の1手法として、状態空間モデルがある。
特に今回扱うのは、一般的に「回帰成分を含むローカルレベルモデル」、と呼ばれるモデルで、式で表すと下記の通り。
線形回帰と違うのは、回帰直線の傾きと切片が、1つ前の時点の傾き、切片にノイズが加わって得られる、としている点である。

 {
\begin{eqnarray} 
y_t &=& \alpha_{t} + \beta_{t} x_t + \varepsilon_t \\
\alpha_{t} &=& \alpha_{t-1} + \eta_{\alpha,t} \\
\beta_{t} &=& \beta_{t-1} + \eta_{\beta,t} \\
\varepsilon_t &\sim& N(0, \sigma^2_\varepsilon) \\
\eta_{\alpha,t} &\sim& N(0, \sigma^2_{\eta_{\alpha}}) \\
\eta_{\beta,t} &\sim& N(0, \sigma^2_{\eta_{\beta}}) \\
\end{eqnarray} 
}

回帰直線の切片と傾きが「状態」と呼ばれる潜在変数である。各時点の状態は1時点前の状態によって決まる。
つまり予測の際は直近のデータの影響が大きくなる。

状態の推定にはカルマンフィルタと呼ばれる手法が用いられる。詳しくはこちらの記事を参照。 logics-of-blue.com

カルマンフィルタの要点としては、観測値が得られるたびに状態の推定値が補正され、補正された状態を用いて次の時点における状態と観測値の予測が行われる、ということになる。

状態空間モデルによる一期先予測の結果は下図の通り。

f:id:kibya:20200107075456p:plain

一期先予測誤差(MAE)は1.81程度となり、線形回帰よりはっきり改善している。

予測結果の比較

2つのモデルの違いを見るため、t=47, 48における予測結果を比べてみる。

まず、線形回帰の場合。この場合は過去の全ての観測値を同等に扱って予測を行う。結果、t=47, 48のいずれでも実際の観測値より大幅に小さい値を予測してしまっている。

f:id:kibya:20200107080906p:plainf:id:kibya:20200107080925p:plain

次に、状態空間モデルの場合。t=47では予測誤差が大きいのは同じだが、t=48の予測誤差は比較的小さくなっている。 これは、t=47の観測値を用いた状態の補正が行われているためである。つまり、直近の観測値は予測値よりだいぶ大きかった、という情報をモデルに与えることで、次の時点の予測値を大きくするように補正が働く

f:id:kibya:20200107081003p:plainf:id:kibya:20200107081013p:plain

アニメーションも載せておく。

f:id:kibya:20200111084459g:plainf:id:kibya:20200111084530g:plain
線形回帰(左)、状態空間モデル(右)

まとめ

時系列データの分析には状態空間モデルなど時系列モデリングの適用をしましょう。