Perplexityとは

　自然言語処理で良く使われるメトリクスのひとつであるperplexityについて説明します。

定義

分布pに従って生成されるデータxがある。xの分布をqでモデリングする。
このとき確率モデルqのperplexityは、

と表される。

ただし、xiはpから生成したデータで、Nはxiの個数である。

直感的な説明

• 指数部はクロスエントロピーになっている
• モデルにデータを見せたときに、どれくらい"驚く"かを表す
• いいモデルほどperplexityは小さい

サンプル

サイコロの目の分布をモデリングしてみて、そのモデルのperplexityを計算してみます。
まずサイコロを用意します。サイコロは六面で各面が出る確率は一様とします。

真の分布はわからないので、30回サイコロを振ってみて観測化された結果から分布を推定してみます。このモデルのperplexityをテストデータ1,000個を使って計算してみると、6.358となります。

import numpy as np
from collections import Counter

np.random.seed(1234)


N = 1000     # test data size
M = 30       # train data size

test_data = np.random.randint(6, size=N)
train_data = np.random.randint(6, size=M)

p = Counter(test_data)    # test distribution
q = Counter(train_data)   # proposed distribution

exponent = 0
for i in range(6):
  exponent -= p[i] / N * np.log2(q[i] / M)

print (2 ** exponent)

次に訓練データを増やしてモデリングしてみます。サイコロを1,000回振ってみます（上記のコードのMを1000に変更してプログラムを実行）。
perplexityは5.998となり、モデルの良さが改善されたことがわかります。

オントロジーとタクソノミーの違い

オントロジー（Ontology）とは

• type, properties, relationshipの集合から構成される世界を記述するモデルのこと。
• Grakn Knowledge Modelというモデルではオントロジーを以下の4つのコンポーネントで記述する。
• entity: もの。例: 人、男、女
• relation: entity同士の関連。例: 親子関係
• role: relationにおける役割。例: "結婚"というrelationには、"夫"、"妻"というroleがある
• resource: entity/relationに紐づく属性。 例:名前、日付

タクソノミー（Taxonomy）とは

• entityを親子関係を使って階層的に並べて表現するモデル。
• 例: 植物 - 被子植物 - 単子葉類 - ユリ目 - ユリ科 - チューリップ属

両者の違い

• Ontologyの方がTaxonomyより複雑な表現をすることができる。
• Taxonomyは木構造、Ontologyはグラフ構造になっているとイメージしておくとよさそう。
• Taxonomyはentityの分類や類似度を測るときに使われるのに対して、Ontologyは知識表現に使われる。Ontologyを使って知識を表現しておくと、マシンを使って新しい知識を推論することができる。

参考
https://www.quora.com/Whats-the-difference-between-an-ontology-and-a-taxonomy
https://blog.grakn.ai/what-is-an-ontology-c5baac4a2f6c

はじパタ第5章 k最近傍法（kNN法）

　平井有三先生の「はじめてのパターン認識」第5章を読んだ。
kNN自体はもちろん知っていたが、理論的な考察についてはほとんど知らなかったので勉強になった。

面白かったところ

• 式を使ったボロノイ図の定義
• kNNとベイズ誤り率の関係
• kNNにおける次元の呪いの話

理解度チェックリスト

• 最近傍法とは？
• k最近傍法とは？
• ボロノイ図とは？
• kNNを使うときに問題となる次元の呪いについて説明せよ。
• 予測データと学習データの距離はどうなる？
• kNNの計算量は？
• kNNの計算量を減らすための工夫をいくつか説明せよ。

はじパタ第4章 確率モデルと識別関数

　平井有三先生の「はじめてのパターン認識」を読んでいる。今日は第4章を読んだ。


悩んだところ

• ピマ・インディアンデータのROCカーブの説明が間違ってるのでは？と思ったら、病気の人は陰性、健康な人は陽性というデータ定義だったらしい（つまり病気の検査をするときの一般的な陽性・陰性の定義と逆だった）。よく見ると注釈に書いていた。

面白かったところ

• 共分散行列の固有ベクトル方向にデータを回転して、固有値でスケーリングする技は知っていたが、「白色化」という名前を初めて知った。
• 独立であることと、相関がないことは同値ではないらしい。独立なら相関はないが、相関がないからといって独立とは限らない。
• マハラノビス距離は知っていたが多次元正規分布の指数の部分に現れているのは気づかなかった。確かに指数の部分はxが現れる確率密度を表すものなのでマハラノビス距離が遠いものは現れにくいということに対応している。
• 最尤推定を使った正規表現のパラメータ推定は1回は自分でちゃんと計算しておいた方がよい。
• 正規分布を使った識別関数が、条件を加えることで、二次局面、線形、最近傍法と変化する話が面白かった。
• 2つのクラスの共分散行列は同じと仮定しないといけない場合、一般にはそうでないことが多い。そのような場合は、クラスをまとめたものを共通の共分散行列として使うといいらしい。

試してみたこと

Pythonで白色化をしてみた。

import numpy as np
import pandas
from matplotlib import pyplot as plt

mean = [0.0, 0.0]
sigma = [[1, 2], [2, 5]]
data = np.random.multivariate_normal(mean, sigma, 300)
x, y = data.T
plt.plot(x, y, 'o')

df = pandas.DataFrame(data)
u = df.mean()
sigma = df.cov()
lam, S = np.linalg.eig(sigma)

assert np.linalg.norm(np.dot(np.dot(S.T, sigma), S) - np.diag(lam)) < 1e-9

whitening = np.dot(np.diag(lam**(-1/2)), S.T)
data_whitened = np.dot(whitening, (df - u).T).T
x, y = data_whitened.T
plt.plot(x, y, 'o')
Out[88]:

理解度チェックリスト

• データの前処理についてそれぞれの手法を説明せよ。
• 標準化
• 無相関化
• 白色化
• パラメトリックモデルとノンパラメトリックモデルの違いは？
• 多次元正規分布に従うn個のデータを使って、分布のパラメータを最尤推定せよ。

はじパタ第3章 ベイズの識別規則

　平井有三先生の「はじめてのパターン認識」を読んでいる。今日は第3章を読んだ。


 
面白かったところ

• ベイズのMAP推定のところで、事後確率 = 修正項 × 事前確率というのは当然知っていたけど、修正項のところの概念的な意味づけが面白かった。
• 図3.1 識別境界とベイズ誤り率、図3.6 ROC曲線と等損失直線群の図が考えさせられた。特に図3.6では最適動作点を探す方法を新しく知れたのでよかった。
• 尤度比のところの話は単純に式変形としてみるのではなく、何を意味しているのかを考えてみると面白かった。

チェックリスト

• ベイズの最大事後確率に基づく識別規則を式を書いて説明せよ。
• 尤度比を使って識別規則を表現せよ。
• ベイズの識別規則の誤り率を計算せよ。（2クラス識別で）
• 誤りを犯すことによって発生する危険性がクラスによって異なる事例を述べよ。
• そういう場合はどうすればよいか？
• リジェクトとはなにか？
• ROCカーブはどのようにプロットするか説明せよ。
• 以下の語句の意味を説明せよ。
• false positive rate
• true positive rate
• precision
• recall
• accuracy
• F-measure
• AUCとは何か？
• とりうる値は？どういう値だといい？
• ROCカーブから最適動作点はどうやって選べばよいか？

はじパタ第2章 識別規則と学習法の概要

　平井有三先生の「はじめてのパターン認識」を読んでいる。今日は第2章「識別規則と学習法の概要」を読んだので学習メモを残しておく。

理解度チェック

• アフィン関数とは？
• アフィン関数を線形関数にするためによく使われる方法は？
• バイアスとは？
• バリアンストは？
• バイアスとバリアンスの関係は？
• 次の手法について説明せよ。
• holdout法
• cross validation法
• leave one out法
• bootstrap法

悩んだところ
バイアス項とバリアンス項（分散項）の導出のところで悩んだ。
一般にE[X²] = E[X]E[X]が成り立つと思い込んでしまっていて、バリアンス項って消えるんじゃない？何これ？となっていた..。
XとYが独立な場合はE[XY] = E[X]E[Y]ですが、XとXは独立じゃないので、バリアンスの項は消えませんね。

基本中の基本を忘れてました。
例えばサイコロに{-2, -1, 0, 0, 1, 2}という数字を書いて、2回投げたときにそれぞれ出る値をX、Yとすると、E(XY) = E(X)E(Y) = 0となりますが、E(XX)は0にはなりませんね。

試してみたこと
N個のデータからN回復元抽出（sampling with replacement）すると、あるデータが少なくとも1回は抽出される確率は1-e^-1となるってのが面白かったので試してみた。1-e^-1は 劣モジュラ関数の最適化や回路の時定数などでも出現する値で興味深い。

import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt


def bootstrap_unique_freq(n):
 x = np.random.randint(0, n, n)
 return len(set(x)) / n


if __name__ == '__main__':
 t = 100  # iteration
 ns = []  # sample size
 xs = []  # unique sample freq

 for n in [10, 100, 1000, 10000]:
  ns.extend(n for _ in range(t))
  xs.extend(bootstrap_unique_freq(n) for _ in range(t))

 sns.boxplot(x=ns, y=xs)
 plt.show()

箱ひげ図を書いてみた。

x軸は元のデータ数、y軸はブートストラップサンプルしたときのユニークなデータ数。

はじパタ第1章 はじめに

　平井有三先生の「はじめてのパターン認識」を読み始めた。パッと目次を見た感じだとだいたい知ってることが多そうだったが、基礎は大事なので丁寧に読み進めたい。
　今日は第1章「はじめに」を読んだ。100次元の超立方体は毬栗のような形をしているという話がおもしろかった。
　以下に理解度チェックリストを作っておく。

• よいパターン認識とは何かいかのキーワードを使って述べよ。
• 特徴抽出
• 識別規則の学習
• 般化能力
• 以下の尺度の定義、具体例を述べよ。
• 名義尺度
• 順序尺度
• 間隔尺度
• 比例尺度
• ダミー変数とはなにか？
• 次元の呪いとは？
• d次元超立方体について以下を計算せよ。
• 頂点の数はいくつ？
• 辺の数は？
• 超立方体を構成するm次元超平面(1 < m < d)の個数はいくつ？

生成モデルと識別モデルの違い

生成モデル（generative model）と識別モデル（discriminative model）の違いが理解できたのでメモっておく。

よくある定義
生成モデルは

識別モデルは

のようにyを予測する。
これだけ聞いても何のことやらよく分からん。

直感的な説明
識別モデルはクラスの境界を学習する。与えられたデータxが生成される確率は特に考慮せず、単に境界面のみを学習するモデル。SVM、決定木などは識別モデル。

生成モデルは、それぞれのクラスにおけるデータの分布を学習するモデル。生成モデルはグラフィカルモデルを使って表される。

これは分かりやすい。これ聞いた後だと式の意味も分かる。

直感的な説明と式を結びつけて考える
識別モデルは、

のように予測を行うモデル。例えば決定木の場合はデータxが与えられると、終端nodeのクラス分布を見ることで、そのxがクラスyに属する確率が分かる。そして尤もらしいクラスを予測値として用いる。このようにp(y|x)を最大化しているが、データxが生成される過程や確率に関しては何もモデリングしていない。

これに対して生成モデルは、

のように予測を行う。例えばベイジアンネットワークやマルコフ確率場などを用いたモデルの場合はp(y)、p(x|y)を定義し、ベイズの定理を使って事後確率が最大となるようなyを予測値として用いる。つまり、データxの生成過程やクラスごとのxの分布というものをモデリングすることで、結果を予測するようなモデルになっている。

sparseなindicator vectorからユーザの類似度を測る

　ユーザがadをclickした/しないでindicator vectorを作ってユーザ間の類似度を測りたかった。僕の直感だと{0, 1}のベクトルの距離なのでハミング距離だろうと思ったけど、コサイン距離のほうがいいらしい。

　確かにハミング距離だと全く嗜好の違うユーザだろうが、嗜好がほとんど同じユーザだろうが、x個の次元が異なると距離がx/|dim|になって嬉しくない。これに対してコサイン距離を使うと、嗜好が似ているユーザの場合は似ている方向に角度が寄るので距離が緩和される。

　直感的にもそうなるのは理解できるが、実際に計算して確かめてみた。

from scipy.spatial.distance import hamming, cosine

x1 = [1,1,1,0,0,0,0,0,0,0]
y1 = [0,1,1,1,0,0,0,0,0,0]

y2 = [0,1,0,0,0,0,0,0,0,0]
x2 = [1,0,0,0,0,0,0,0,0,0]

print(hamming(x1, y1))  # 0.2
print(hamming(x2, y2))  # 0.2

print(cosine(x1, y1))  # 0.333333333333
print(cosine(x2, y2))  # 1.0

信頼区間付きの線形回帰予測

　pythonで線形回帰の信頼区間付きpredictionをしたかったけど、sklearnには該当する機能が無かった。Gradient Boosting Regressionだと出来るのだけど[2]、単純なLeast Squareでやりたかったのでどうやるのか調べていたら、stasmodelsを使えばいい[1]らしい。

サンプルソース

import numpy as np
import statsmodels.api as sm
from statsmodels.sandbox.regression.predstd import wls_prediction_std
import matplotlib.pyplot as plt

np.random.seed(1234)
nsample = 30
x1 = np.linspace(0, 20, nsample)
x = np.column_stack((x1, x1**0.5))
x = sm.add_constant(x)
beta = [1, 0.5, 0.5]
y = np.dot(x, beta) + np.random.normal(size=nsample)
res = sm.OLS(y, x).fit()
prstd, iv_l, iv_u = wls_prediction_std(res)

plt.plot(x1, y, 'o', label="Data")
plt.plot(x1, res.fittedvalues, 'r--')
plt.plot(x1, iv_u, 'r--', label="Ordinary Least Squares")
plt.plot(x1, iv_l, 'r--')
plt.legend(loc="best");
plt.show()

wls_prediction_stdは、学習データに対する予測値の有意水準5%信頼区間を返してくれる。有意水準は引数alphaで自由に設定できる。また、学習データ以外のデータに対する予測を行いたい場合は、引数exogにデータを渡せばよい[3]。

結果

参考
[1] https://www.reddit.com/r/MachineLearning/comments/3raivl/code_to_calculate_confidence_interval_for_linear/
[2] http://scikit-learn.org/stable/auto_examples/ensemble/plot_gradient_boosting_quantile.html
[3] https://github.com/statsmodels/statsmodels/blob/master/statsmodels/sandbox/regression/predstd.py

多項ロジットの理論と実装

問題設定
訓練データ(x_n, y_n), x_n ∈ R ^D, y_n ∈ {1,2,...C}, n=1,2,...,Nを使って、マルチクラス識別問題を解くモデルを考えます。

多項ロジットのモデル
多項ロジットでは以下のモデルを使います。

分子は訓練データxの特徴量の線形結合をロジスティック関数に入れたものです。
分子は非負の値を取ります。

上の式はxがクラスyに属する確率を表します。
よって、xがクラス1に属する確率、xがクラス2に属する確率、...をすべて足すと1.0になっていて欲しいです。この条件を満たすように分母の値で割って正規化を行います。

パラメータの学習
多項ロジットのパラメータは、

です。ここでCはクラスの個数、Dはデータの次元数です。
つまりクラスの数だけ重みベクトルが存在します。

これらのパラメータを学習するために、以下の対数尤度関数を最大化します。

対数尤度関数の偏微分を計算します。
まず第一項の偏微分は以下のようになります。

[]は中の条件が真のときは1、偽のときは0になるような関数です。

第二項の偏微分は以下のようになります。

よって全体の偏微分は、

となります。

偏微分が求まったので、θを適当な値に初期化して、勾配法で対数尤度関数を最適化することができます。

実装
Pythonで実装しました。

import numpy as np

class MultinomialLogit:

 def __init__(self):
  self.C = 0
  self.N = 0
  self.D = 0
  self.theta = None

 def prob(self, X, y):
  return np.dot(self.theta[y], X) / sum(np.dot(self.theta[c], X) for c in range(self.C))

 def logLiklihood(self, X, y):
  return sum(np.log(self.prob(X[n], y[n])) for n in range(self.N))
 
 def fit(self, X, y, iter=200):
  self.C = len(set(y))
  self.N = len(X)
  self.D = len(X[0])+1   # bias term is added
  self.theta = np.random.rand(self.C, self.D);

  bias = np.ones((len(X), 1))
  X = np.hstack((X, bias))

  def delta(a, b):
   return 1.0 if a == b else 0.0
  
  for i in range(iter):
   for c in range(self.C):
    gradient = sum((delta(c, y[n]) - self.prob(X[n], c)) * X[n] for n in range(self.N)) 
    self.theta[c] += gradient
   print ("i=%d %f" % (i, self.logLiklihood(X, y)))
    
 def predict(self, X):
  bias = np.ones((len(X), 1))
  X = np.hstack((X, bias))
  return [np.argmax([self.prob(x, c) for c in range(self.C)]) for x in X]

 def score(self, X, y):
  Ypred = self.predict(X)
  comp = Ypred == y
  return sum(comp) / len(comp)

Irisで実験
Irisのデータを使って実験してみました。
テスト識別率は80%後半程度です。sklearnのLogisticRegressionだと100%が出るので全然ダメな感じですが、

• iterationごとに対数尤度関数が上昇すること
• 対数尤度関数が上がるとテスト識別率も上がること

を確認できたのでそれっぽい実装は出来てそうです。学習率を設定したり、 ヘッセ行列を計算してニュートン法系の最適化を行ったり、確率的勾配法を使ったりなど改善の余地はたくさんありそうです。

if __name__ == '__main__':
 from sklearn import datasets
 from sklearn import cross_validation
 
 iris = datasets.load_iris()
 X = iris.data
 y = iris.target

 Xtrain, Xtest, Ytrain, Ytest = cross_validation.train_test_split(X, y, test_size=0.3)
 logit = MultinomialLogit()
 logit.fit(Xtrain, Ytrain, iter=1000)
 print (logit.score(Xtest, Ytest))

最急降下法とニュートン法の比較

autogradを使って、最急降下法（steepest descent method）とニュートン法（Newton's method）の比較をしてみた。

最急降下法のおさらい

1. 適当な初期点xを選ぶ
2. xでの勾配ベクトルを求めて勾配ベクトル方向（最も急な方向）に点を移動させる
3. 収束するまで2.を繰り返す

ニュートン法のおさらい
関数fの点x_nのまわりでのテイラー展開は以下のようになる。

ただし、3次以上の項は無視している。

最適化問題で求めたいのはfの停留点である。

ということで、上のテイラー展開を使って微分が零になる点を近似的に計算してみる。

なので、fをxで微分したものと、fをΔxで微分したものは同じ。

よって、

となる。この値が0になるようなΔxを求めると、

である。

つまりこの方向に進めば、fの微分が0になりそうな場所に向かって、xを逐次的に更新することができる。

ニュートン法のアルゴリズムを以下に示す。

1. 適当な初期点xを選ぶ
2. 上記の手順でx = x + Δxと更新する
3. 収束するまで2.を繰り返す

実験
autogradを使って、最急降下法とニュートン法を実装する。
必要なライブラリのimport。

import autograd
import matplotlib.pyplot as plt
import numpy as np

最急降下法の実装。
グラフに描画したいので、xの最適値だけじゃなくて、xの更新履歴も返すようにしている。

def steepest_descent(f, x0, alpha=0.01, eps=1e-9, maxitr=500):
  x = x0
  history = [np.array(x)]
  for _ in range(maxitr):
    g = autograd.grad(f)
    if (np.linalg.norm(g(x)) < eps):
      break
    x -= alpha * g(x)
    history.append(np.array(x))
  return x, history

ニュートン法の実装。

def newton(f, x0, alpha=0.01, eps=1e-9, maxitr=500):
  x = x0
  history = [np.array(x)]
  for _ in range(maxitr):
    g = autograd.grad(f)
    H = autograd.hessian(f)
    if (np.linalg.norm(g(x)) < eps):
      break
    x += alpha * np.linalg.solve(H(x), -g(x))
    history.append(np.array(x))
  return x, history

軸がずれた楕円形の関数を使って、両手法を比較してみる。

def f(x, y):
  s = 2*x + y
  t =  -x + y
  return 4 * s**2 + 3 * t**2

# draw contour of f(x, y)
delta = 0.01
x = np.arange(-2.2, 0.2, delta)
y = np.arange(-2.2, 0.2, delta)
X, Y = np.meshgrid(x, y)
Z = f(X, Y)
plt.contour(X, Y, Z, 80, colors='black')

# execute steepest descent
x0 = np.array([-1.5, -2.0])
wrapper = lambda xs: f(*xs)
optimal, history = steepest_descent(wrapper, x0, alpha=0.001)
for x in history:
  plt.plot(x[0], x[1], 'ro', markeredgewidth=0.0)

# execute newton method
x0 = np.array([-1.5, -2.0])
wrapper = lambda xs: f(*xs)
optimal, history = newton(wrapper, x0)
for x in history:
  plt.plot(x[0], x[1], 'bo', markeredgewidth=0.0)
  
plt.show()

結果は以下のとおり。
最急降下法（赤点）の場合は、等高線と直交する方向にxが更新されているのが分かる。
それに対して、ニュートン法（青点）の場合は、近似的に計算した関数の停留点に向かって真直ぐxが更新されているのが分かる。

TensorFlow入門（2）Deep MNIST for Experts

　TensorFlowのtutorialをやってみた。MNISTをDeep Learningで学習させる回。識別率99%を超えるモデルを作ることが出来る。重要そうなところをメモしておく。

元ネタ
Deep MNIST for Experts

TensorFlowについて
重い計算はC++のバックエンドで行う。
バックエンドとのコネクションを"セッション"と呼ぶ。

ソースコード

import tensorflow as tf
import input_data

def weight_variable(shape):
  initial = tf.truncated_normal(shape, stddev=0.1)
  return tf.Variable(initial)

def bias_variable(shape):
  initial = tf.constant(0.1, shape=shape)
  return tf.Variable(initial)

def conv2d(x, W):
  return tf.nn.conv2d(x, W, strides=[1, 1, 1, 1], 
                      padding='SAME')

def max_pool_2x2(x):
  return tf.nn.max_pool(x, ksize=[1, 2, 2, 1],
                        strides=[1, 2, 2, 1], padding='SAME')

mnist = input_data.read_data_sets('MNIST_data', one_hot=True)

sess = tf.InteractiveSession()
x = tf.placeholder("float", shape=[None, 784])
y_ = tf.placeholder("float", shape=[None, 10])

x_image = tf.reshape(x, [-1,28,28,1])
W_conv1 = weight_variable([5, 5, 1, 32])
b_conv1 = bias_variable([32])
h_conv1 = tf.nn.relu(conv2d(x_image, W_conv1) + b_conv1)
h_pool1 = max_pool_2x2(h_conv1)

W_conv2 = weight_variable([5, 5, 32, 64])
b_conv2 = bias_variable([64])
h_conv2 = tf.nn.relu(conv2d(h_pool1, W_conv2) + b_conv2)
h_pool2 = max_pool_2x2(h_conv2)

W_fc1 = weight_variable([7 * 7 * 64, 1024])
b_fc1 = bias_variable([1024])
h_pool2_flat = tf.reshape(h_pool2, [-1, 7*7*64])
h_fc1 = tf.nn.relu(tf.matmul(h_pool2_flat, W_fc1) + b_fc1)

keep_prob = tf.placeholder("float")
h_fc1_drop = tf.nn.dropout(h_fc1, keep_prob)

W_fc2 = weight_variable([1024, 10])
b_fc2 = bias_variable([10])
y_conv=tf.nn.softmax(tf.matmul(h_fc1_drop, W_fc2) + b_fc2)

cross_entropy = -tf.reduce_sum(y_*tf.log(y_conv))
train_step = tf.train.AdamOptimizer(1e-4).minimize(cross_entropy)
correct_prediction = tf.equal(tf.argmax(y_conv,1), tf.argmax(y_,1))
accuracy = tf.reduce_mean(tf.cast(correct_prediction, "float"))

sess.run(tf.initialize_all_variables())
for i in range(20000):
  batch = mnist.train.next_batch(50)
  if i%100 == 0:
    train_accuracy = accuracy.eval(feed_dict={
        x:batch[0], y_: batch[1], keep_prob: 1.0})
    print "step %d, training accuracy %g"%(i, train_accuracy)
  train_step.run(feed_dict={x: batch[0], y_: batch[1], keep_prob: 0.5})

print "test accuracy %g"%accuracy.eval(feed_dict={
    x: mnist.test.images, y_: mnist.test.labels, keep_prob: 1.0})

以下、上記ソースコードで記述したモデルについて説明。

Weight Initialization
wはtruncated normal distributionの乱数。
bは定数。

Convolution and Pooling

Convolutionのstrideは1*1。padding='SAME'とすることで、入力と出力のサイズが同じになる。

Poolingは2*2のブロックのMax Pooling。

First Convolutional Layer

28 * 28 * 1チャネルの入力を

5*5のパッチでconvolutionして28 * 28 * 32 チャネルにした後、

Max Poolingして、14 * 14 * 32チャネルにする。

Second Convolutional Layer

14 * 14 * 32チャネルの入力を

5*5のパッチでconvolutionして14 * 14 * 64チャネルにした後、

Max Poolingして、7 * 7 * 64チャネルにする。

Densely Connected Layer

7 * 7 * 64チャネルの入力をflattenしてベクトルにした後、

1024次元のベクトルにする。

Dropout
オーバーフィッティング対策のために出力層の前にdropoutを入れる。
keepするneuronの割合をplaceholderにしておくことで、学習時と識別時でdropoutのあり/なしを切り替えることが出来る。

Readout Layer

softmaxする。

TensorFlow入門（1）MNIST For ML Beginners

　TensorFlowのtutorialをやってみた。
機械学習入門者向けのtutorialでMNISTを例題にNNを学習させるやり方が記載されていた。
NNはnum of hidden layer=0、activation function=softmaxというシンプルなモデルを使っていた。自分用にメモを残しておく。

元ネタ
MNIST For ML Beginners

softmaxについて

• 対象がそれぞれのクラスに属する確率を求めたいときに使う
• NNの最終レイヤーはsoftmaxを使うことが多い

TensorFlowについて
numpyは重い数値計算をpythonの外で行う
=> 結果をpythonに戻すときにオーバーヘッドがある
=>TensorFlowはpythonの外で実行したい計算処理の相互作用をグラフで記述できる

placeholder: 計算開始時に入力として渡す値 （データの特徴ベクトル、教師信号用）
Variable: 計算中時に利用する変更可能な変数（学習させるパラメータ用）

ソースコード

import input_data
import tensorflow as tf

mnist = input_data.read_data_sets("MNIST_data/", one_hot=True)

x = tf.placeholder(tf.float32, [None, 784])
W = tf.Variable(tf.zeros([784, 10]))
b = tf.Variable(tf.zeros([10]))
y = tf.nn.softmax(tf.matmul(x, W) + b)

y_ = tf.placeholder(tf.float32, [None, 10])
cross_entropy = -tf.reduce_sum(y_*tf.log(y))

train_step = tf.train.GradientDescentOptimizer(0.01).minimize(cross_entropy)
init = tf.initialize_all_variables()

sess = tf.Session()
sess.run(init)
for i in range(1000):
  batch_xs, batch_ys = mnist.train.next_batch(100)
  sess.run(train_step, feed_dict={x: batch_xs, y_: batch_ys})

correct_prediction = tf.equal(tf.argmax(y,1), tf.argmax(y_,1))
accuracy = tf.reduce_mean(tf.cast(correct_prediction, "float"))
print sess.run(accuracy, feed_dict={x: mnist.test.images, y_: mnist.test.labels})

中国料理店過程

中国料理店過程とは？

• 離散時刻で確率的に変化する過程
• 中国料理店に客がやってきたときに、どの円卓に座るかを考える
• 店には無数の円卓があり、円卓には無数の人が座れるものとする
• 最初の客は空の円卓に座る
• 2番目以降の客はランダムに円卓を選ぶが、座っている人の数が多い円卓を好む
• 数学的には[1,2,3,...,n]のpartitionがランダムに作成されることになり、数学者はこのpartitionの確率分布に興味を持つらしい
• ノンパラメトリックなベイズメソッドに応用される

モデリング

より一般化されたモデルもあるみたいだが、一番簡単なモデルを考える。

最初の客は空の円卓に座る。

n番目（n > 1）の客は、|b_i|/nの確率で円卓b_iに座る。もしくは、1/nの確率で空の円卓に座る。ただし、|b_i|は円卓b_iに座っている客の人数を表す。

サンプルプログラム

これだけだと何が嬉しいかわからない。ベイズにどうやって応用するか勉強したい。

from random import random
from bisect import bisect_left
import numpy as np

n = 100       # number of customers
block = []    # partition

for i in range(n):
    if not block:
        block.append([i])
    else:
        x = random()
        p = np.cumsum([1.0*len(b)/(i+1) for b in block])
        pos = bisect_left(p, x)
        if pos >= len(block):
            block.append([i])
        else:
            block[pos].append(i)
    print block

データサイエンス系ボキャブラリーまとめ

　データサイエンス系（その他もろもろ）のボキャブラリーが乏しいので、メモ書き程度に書き溜めていこうと思う。

アドホック分析
特定の目的のために実施される単発の分析。ad-hocはその場限りの、特定の目的のための、という意味。

BlueKai
あらゆるオンライン上のデータを集め、DSP/SSPなどに販売することをビジネスモデルとしたスタートアップ。DMPベンダー。2014年オラクルに買収された。

デモグラフィック
性別、年齢、職業、年収などのユーザー属性情報。
デモグラと略される。

ジオグラフィック
国、都道府県、都市などのユーザーの地理的情報。
ジオグラと略される。

サイコグラフィック
性格、嗜好、価値感などのユーザーの心理的情報。

ランディングページ
ネット上の広告やリンクをクリックしたときに表示されるページ。
略称でLP、ランペとも呼ばれる。

コンバージョン
WEBサイト上で獲得できる最終的な成果。
商品購入、会員登録、資料請求など。

CPA
Cost Per Acquisition。顧客獲得単価。
1コンバージョンあたりの広告コスト。

フィジビリティスタディ
実行可能性調査。F/S。
事業、プロジェクトが実現可能か事前に調査すること。

クリエイティブ
広告素材全般のこと。
バナー広告の画像、テキスト広告のテキスト、キャッチコピーなど。

インプレッション
広告の露出（掲載）回数のこと。

CTR
Click Through Rate。クリック率。
クリック数 / インプレッションで算出する。

MAU
Monthly Active User。
ある月に実際にサービスを利用したユーザーの数。

DAU
Daily Active User。
ある日に実際にサービスを利用したユーザーの数。

ROI
Return on Investment。
利益 / 投資額で算出。

入稿
広告主からLP、バナーなどのクリエイティブを受け取ること。

Cookie Sync
複数のドメインでクッキー情報を共有し、サイトをまたいだユーザーの特定を可能とする仕組みのこと。

Pixel Piggyback
1*1ピクセルの透過gifを使ってCookie Syncを行う技術。
　1. サイトAに、サイトBから透過gifを読むタグを設置
　2. ユーザがサイトAを訪れると、サイトBにサイトAでのユーザ情報（ID）が送られる
　3. サイトBはクッキーを発行し、発行したIDとサイトAのIDを結びつける

pixel piggybackの例：
<img alt="" src="http://xxxx/yyy.gif?uid=zzzzzzzz" height="1" width="1" />
xxxxドメイン側では、uidパラメータの取得、クッキーの発行（xxxx側のuid発行）、両方のuidの紐付けを行ったあと、空のgifを返す。

TD
Title and Description。
広告文のタイトルと説明文のこと。

CVR
Conversion Rate。
コンバージョン数 / 広告のクリック数で算出。

直帰率
1ページのみ閲覧してサイトを去った人の割合。
そのページがセッションの始まりだった場合、何割の人がすぐにセッションを終了（他サイトに移動、ブラウザを閉じるなど）しているかを表す。

離脱率
そのページがセッションの最後になった割合。
そのページのPV数のうち、そのページがセッションの最後になった割合。

UU
Unique User数。
特定の期間内にウェブサイトを訪問したユニークなユーザ数。
"ユニークな"の測定はIPアドレス、IPアドレスとUser Agentの組み合わせで行われることが多い。

LTV
Life Time Value。顧客生涯価値。
顧客ひとりが生涯で商品・サービスに支払った金額の合計値。

CRM
Customer Relationship Management。顧客関係管理。
顧客の属性、企業との接点を管理・可視化・分析し、顧客の長期的価値を最大化するマーケティング手法。

リテンション
人材領域の文脈では人材の維持・確保、マーケティング領域の文脈では既存顧客維持を意味する。

リスティング広告
検索エンジンの検索結果に連動して表示される広告のこと。

ディスプレイ広告
画像形式の広告のこと。バナー広告。

インフィード広告
SNSやキュレーションサイトなどで、フィードの間に表示される広告。ネイティブ広告の一種。

SEM
Search Engine Marketing。
具体的には、検索エンジンの表示順をあげるSEOやキーワード連動広告への広告出稿など。

コンテンツマーケティング
ブログ、ポッドキャスト、動画などのコンテンツで見込み客をサイト内に引きつけ、資料請求や商品購入を行ってもらうマーケティング手法。

エンゲージメント
親和性。
サービスの内容が、ユーザの求めているものと一致している場合、「ユーザとのエンゲージメントが高い」という。

オセロのAIを作る（2）

　Monte Carlo Tree Searchを使ったオセロAIを実装しました。

試しに対戦してみましたが、私（オセロ初心者）よりは少し強い感じです。
また、以前実装した中で最強プレイヤーだったDfsStoneNumPlayer（評価関数を石の数としたゲーム木探索プレイヤー）よりも強いです。

　せっかくなので、市販のオセロAIと対戦させてみよう。
ということで、iphoneアプリのHAYABUSAと対戦させてみました。

まずは最弱レベルであるLOCAL-EASYを相手にしてみました。
MctsPlayerというのが私が作ったAIの名前です。

MctsPlayer ●44 - ○20 HAYABUSA(LOCAL-EASY)
MctsPlayer ○50 - ●14 HAYABUSA(LOCAL-EASY)
MctsPlayer ●53 - ○11 HAYABUSA(LOCAL-EASY)
MctsPlayer ○41 - ●23 HAYABUSA(LOCAL-EASY)
MctsPlayer ●41 - ○23 HAYABUSA(LOCAL-EASY)

最弱モードと言えども、圧勝できると嬉しいものです。
相手を連続パスに追い込むという場面も多々見られ、なかなかいい感じです。

次にLOCAL-NORMAL（LOCAL-EASYより１つ上のレベル）と対戦させてみました。
ちなみに私が自分でプレーした戦績は0勝2敗でした。こいつに勝てると嬉しいですね。

MctsPlayer ○20 - ●44 HAYABUSA(LOCAL-NORMAL)
MctsPlayer ●15 - ○49 HAYABUSA(LOCAL-NORMAL)
MctsPlayer ○24 - ●40 HAYABUSA(LOCAL-NORMAL)
MctsPlayer ●15 - ○49 HAYABUSA(LOCAL-NORMAL)
MctsPlayer ○16 - ●48 HAYABUSA(LOCAL-NORMAL)

はい、完敗〜。
完膚なきまでに打ちのめされてしまいました。
打ちたくない場所にいいように打たされて、最終的にひっくり返されるというパターンで負けていました。

相手の動きを観察すると、

• 敵が打てる場所を減らす
• 敵に空マスの隣を取らせる

ということをやっているようです。
おそらく、disc count、mobility、frontierの加重和を評価関数にしてα-β探索でもやってるのかなと思います。もっと複雑なことやってるかもしれませんが。

オセロのAIを作る（1）

　先日MCTSを実装してから、オセロのAI作りたいなーという気持ちになっていた。

で、盆休みに２日くらいで作ったのがこれ。
オセロの定石知らない私でも簡単に勝てるような弱小AIしか実装できてませんが、これからちょっとずつ強化していこうと思います。

Kaggle参加記：Search Results Relevance

　KaggleのSearch Results Relevanceに参加しました。

About competition

このコンテストでは、検索クエリとそれに対する検索結果がデータとして与えられます。教師ありデータではクエリに対する検索結果のマッチ度を1-4で評価しています。テストデータではこのマッチ度を予測します。

What I did

特徴エンジニアリングに利用するデータがテキストなので、NLPまわりの知識が必要になります。知識ゼロベースからのスタートだったので、とりあえず論文やらチュートリアルやらライブラリの使い方やらを読み漁るところから始めました。

　TF-IDFベクトルのコサイン類似度、Okappi BM25あたりが検索ワード-検索結果のランキングアルゴリズムでよく使われることが分かったので、この２つをC++でカスタム実装してみました。特徴量抽出後はsklearnのランダムフォレストで識別。いい結果出ず。テキストの前処理が甘いと思い、大文字⇒小文字変換、Stop-words除去、stemmingなどを行うもスコアが思ったように伸びず。

　こういうときは、Forumだ！ということで、Forumに上がっているスレッドを全部読みました。starter scriptなるものがあり、それを実行すると自分が書いたものよりいい結果が出て凹む。
しかもそのスクリプトがTF-IDFを抽出した後、次元圧縮し、SVMで識別するだけというシンプルなモデルだったため余計に凹む。しかし何故こんなモデルでうまくいくのだと悩んでいると同じようなことを質問しているユーザがいたため、何か手がかりになるはずと考え理由を考えました。

　そして検索対象がマッチしやすいかそうでないオブジェクトに分けるとうまくいくのではという結論に達しました。実際にデータを見てみると同じ検索クエリが複数回使われていて、平均マッチ度が高いクエリと低いクエリが存在することが判明しました。もっと早めに気づけよという感じです。「特徴エンジニアリングする前にデータをよく知る」ということが出来ていなかったです。

　この知見をベースに、クエリ毎に識別器を作ったり、始めにクエリをk-meansでクラスタリングした後にクラスタ毎に識別器を作ったりして、スコアもまあまあ上がりました。cross validationでいい結果の出たモデルをいくつか選んでスタッキングしてサブミットで試合終了となりました。

What I didn't do but top kagglers did

コンテスト終了後にトップランカーたちがやっていたことを調べました。その中でこれは！というものを書いておきます。1.-4.はテクニック的な話、5.はそもそもの心構え的な話です。

1. コサイン類似度以外にも複数の距離（編集距離、圧縮距離など）を特徴量として使う
2. 検索結果=4のものをクエリに足し合わせて拡張クエリを作る（マッチ度が高い検索結果 = クエリをより詳細に記述したものであるはず）
3. 識別問題ではなく、回帰問題として扱う（評価基準がkappa estimatorなので、ラベルは離散値だけど実測値からの誤差を最小化するとうまくいく）
4. 回帰モデルの予測値をソートし、"1","2","3","4"のスコアを訓練データと同じ割合で付与する（単純に回帰モデルの予測値を整数化してもいい結果はでない）
5. validation時にスコアの悪いものを1つ1つ確認し、なぜスコアが悪いか考えモデルを改良していく

My model after the competition

ということで、上の反省をふまえてモデリングをやり直し、一から実装しなおしました。
まあこんな感じのものです。

ソースコードはこちら。
private score = 0.66830でした。これでも180位くらいなのでまだまだですね。

scikit-learn(5) Gradient Boosting

　以前から気になっていたGradient Boostingについて勉強した。

Kaggleのトップランカーたちを見ていると、SVM、Random Forest、Neural Network、Gradient Boostingの4つをstackingして使っていることが多い。SVMとRFは使えるので、そろそろGradient Boostingも使えるようになっておきたところ。

Gradient Boostingの基本

いいビデオがあった。これを見ると基本的なアイディアはつかめる。

Copyright Information

Prof. Alexander Ihler, Ensembles (3): Gradient Boosting

License: Standard YouTube License

Boostingとは？

・アンサンブル学習の一種

・弱学習器を直列に繋ぐ

・各段の学習器は前段でエラーとなったサンプルに重点を置く

Gradient Boostingとは？

・Boostingモデルの一つ

・以下の手続きを繰り返す

　1. 弱学習器を学習

　2. 残差（観測値と予測値の差）を測定

　3. 残差にフィットするように次の弱学習器を学習

なぜGradient Boostingと呼ばれるか？

残差を使って学習器の予測値を逐次的に更新する

⇒Gradient DescentでMSEを最小化することと等価

※MSEの勾配 ~ 残差

Gradient Boostingの実用的な話

scikit-learnのdeveloperの人のビデオに他のモデルとの比較やパラメータチューニングの話があった。

Copyright Information

Peter Prettenhofer, Gradient Boosted Regression Trees in scikit-learn

License: Standard YouTube License

なぜGradient Boosting？
・Gradient Boostingは有名なデータサイエンス系のコンペティションで頻繁に用いられている

Gradient Boostingの利点
・不均一なデータを扱える（特徴量のスケールが違ってもOK）
・いろいろな損失関数をサポート
・特徴量の非線形な相互作用を自動検出

scikit-learnで使用する際のパラメータチューニング
・n_estimators（boostingのステージ数）=3000で他のハイパーパラメータをgride searchでチューニング
・n_estimatorsを大きくし、learning_rate（学習率；弱学習器を結合するときの係数）をチューニング

サンプルコード

scikit-learnでボストンの家の価格を推定するサンプルコードを書いてみた。
パラメータチューニングはPeter氏が推奨している方法でおこなった。

from sklearn.ensemble import GradientBoostingRegressor
from sklearn.grid_search import GridSearchCV
from sklearn.datasets import load_boston

data = load_boston()

x = data['data']
y = data['target']

# tune hyper parameters first
model = GradientBoostingRegressor(n_estimators=3000)
parameters = {'learning_rate' : [0.1, 0.05, 0.02, 0.01], 
              'max_depth': [4, 6],
              'min_samples_leaf': [3, 5, 9, 17],
              'max_features': [1.0, 0.3, 0.1]}

gscv = GridSearchCV(model, parameters, verbose=10, n_jobs=-1, cv=4)
gscv.fit(x, y)
print "best score=", gscv.best_score_

# tune learning rate with higher n_estimators
model = gscv.best_estimator_
model.set_params(n_estimators=100000)
parameters = {'learning_rate' : [0.1, 0.05, 0.02, 0.01]}
gscv = GridSearchCV(model, parameters, verbose=10, n_jobs=-1, cv=4)
gscv.fit(x, y)
print "best score=", gscv.best_score_