自然言語処理×ナイーブベイズ分類器で羽生さんと羽生くんを分類してみた

# はじめに

機械学習についての理解を促進するため、
データから分類モデルを自動で構築する古典的な方法である、
ナイーブベイズ分類器を実装してみました。
最近は[CloudVisionAPI](https://cloud.google.com/vision/)など専ら画像解析が流行っていますが、
初学者には敷居が高そうだったため、まずは自然言語処理です。

今回は、TwitterAPIを利用してbotアカウントから名言を収集し、
[羽生さん](https://twitter.com/habuyoshibot)と[羽生くん](https://twitter.com/Yuzuru_Hanyubot)を分類する**羽生分類器**を実装しました。

尚、APIクライアントはRuby、分類器はPythonで実装し、
形態素解析にはMeCabを利用しています。

また本来ならばどちらも羽生さんとお呼びしなければならない所、
便宜上フィギュアスケーターの羽生さんを羽生くんと表記させて頂く事をお許し下さい。

# Tweetデータの収集
APIClientは、[Twitter Ruby Gem](https://github.com/sferik/twitter)を使って下記の通り実装しています。
TwitterAPIには[Rate Limits](https://dev.twitter.com/rest/public/rate-limiting)という時間あたりの実行回数制限があり、
大量のツイートを取得しようとすると15分ほどインターバルを挟む必要があります。
適宜コーヒーブレイクをとってください。

## API認証設定
```rb

require 'twitter'

client = Twitter::REST::Client.new do |config|
  config.consumer_key = 'XXX'
  config.consumer_secret = 'YYY'
  config.access_token = 'hoge'
  config.access_token_secret = 'fuga'
end
```

## ツイートの取得
```rb
client.user_timeline('TwitterUserID',{ count: 150}).each_with_index do |tl,i|
	tw = client.status(tl.id)
	tweet = tw.text
        
	# 重複を排除
	if !tweets.include?(tweet)
		puts tweet
	end
end

```

# 形態素解析
さて、先ほど取得したデータを加工します。
後述しますが、今回**ツイートに羽生さんや羽生くんが言いそうな単語が含まれるかどうか**を数量化して分類するため、まずは収集したツイートを品詞分解して、2人のツイートに頻出する単語をそれぞれ上位50件、計100件を分類に用いる変数としてピックアップしています。(実際は重複もあったため91件です。)

## MeCabでツイートを品詞分解
まず今回は、**名詞、動詞、形容詞**のみをカウントの対象としました。
(動詞、形容詞の活用は基本形に直しています)

### 除外ワードリスト

下記の通り人物固有の語彙に無関係と思われる形式名詞や、
動詞化した名詞・形容詞の分かち書きによって発生したと思われる、
カウント対象外とする単語のリストをいくつか設定しています。

```py
ng_noun = ["こと", "の", "もの", "それ", "とき", "、", "､", "。", "｡", "(", ")", "."]
ng_verb = ["する", "いる", "なる", "ある"]
ng_adjective = ["よう"]
```

### 頻出単語のカウント

カウンター付きリスト(タプル)生成には *collections* パッケージが便利です。
また、PythonとMeCabのバインディングには[natto](https://pypi.python.org/pypi/natto-py)を利用しました。

```py

import collections

from sets import Set
from natto import MeCab

def mostFrequentWords(file, num):
  words = collections.Counter()

  f = open(file)
  line = f.readline()
  while line:
    # 名詞:surface="スケート", feature="名詞,一般,*,*,*,*,スケート,スケート,スケート"
    # 動詞:surface="滑れ", feature="動詞,自立,*,*,一段,未然形,滑れる,スベレ,スベレ"
    for node in mecab.parse(line, as_nodes=True):
      features = node.feature.split(",")
    
      if features[0] == "名詞" and node.surface not in ng_noun:
        words[node.surface] += 1
      elif features[0] == "動詞" and features[6] not in ng_verb:
        words[features[6]] += 1
      elif features[0] == "形容詞" and features[6] not in ng_adjective:
        words[features[6]] += 1

    line = f.readline()
    return words.most_common(num


words["hanyu"] = mostFrequentWords("hanyu_train.txt", 50)
words["habu"] = mostFrequentWords("habu_train.txt", 50)

tpl = words["hanyu"] + words["habu"]
vocabulary = set([])
for word in tpl:
  vocabulary.add(word[0])
```

# ナイーブベイズ分類器

ここから若干、数学的背景の解説です。

## ベイズの定理
まず、ナイーブベイズ分類器は確率に基づいた分類器です。
今回求めたいのは、ある文書(ここでは一つ一つのツイート) *d* が与えられた時、
それがどのクラス(羽生さんor羽生くん) *c* に属する確率か高いかどうかです。
これはツイートが与えられた時の条件付き確率として *P(c|d)* と表す事が出来ます。
ただしこの事後確率を直接求めることは難しいため、**ベイズの定理**を用いて計算します。

```math
P(c|d) = \frac{P(c)P(d|c)}{P(d)}

```

ここで、右辺を各クラス、つまり羽生さん／羽生くんそれぞれについて計算し、
ツイートがどちらに属する確率が高いかどうかを求めます。
ただし分母 *P(d)* は分類器を構築してしまえば、クラスとは無関係に一定のため、
**分子のみ計算すればよかろうなのだ**という事になります。

### *P(c)*
今回は羽生さんと羽生くんそれぞれのツイートを100件ずつ取得し、70件を分類器の構築のための訓練データ、
30件を分類器の精度検証のためのテストデータとしています。
*P(c)* はそのクラスの生起確率なので、全訓練データ中の、それぞれの割合になります。
今回訓練データの件数は同じなので、
*P(羽生さん)*、*P(羽生くん)*はどちらも *70 / 140 = 0.5* です。

### *P(d|c)*
さて、この条件付き確率 *P(d|c)* の意味を真っ当に考えると羽生さん、
羽生くんそれぞれが発言しうる単語の種類のその組み合わせで以って、
一つ一つのツイートが生起する確率を求めなければいけないわけですが、到底無理な話です。
そこで、文書分類に適した簡略化したモデルを利用して表します。
ここでは **羽生さんや羽生くんが言いそうな単語の集合 *V* について、
それらが分類に掛けられたツイートに含まれるか含まれないかどうか** を考えます。

#### 多変数ベルヌーイモデル
言った／言わないのような2値をとる確率変数の分布と言えば、**ベルヌーイ分布**です。

```math
{P_{w,c}}^{\delta_{w,d}}(1-{P_{w,c}})^{1-\delta_{w,d}}
```

この指数部分をデルタ関数といい、*w* = *d* の時に0、そうでなければ1が出力されます。
上手く考えられてますね。

ここでは集合 *V* に属するそれぞれの単語 *w* についてベルヌーイ分布を考える、
**多変数ベルヌーイモデル**で *P(d|c)* を表します。

```math
\prod_{w \in V}{P_{w,c}}^{\delta_{w,d}}(1-{P_{w,c}})^{1-\delta_{w,d}}
```

ちなみに上記から多変数ベルヌーイモデルの特徴として下記の2点が読み取れます。

- 文書中の単語の生起回数は考慮されていない
- 文書中に単語が「生起しなかった」という事象が重要視される

### 最尤法
まとめると下記の通り表せるので、右辺を羽生さんと羽生くんそれぞれについて計算し、
値の大きい方 = **どちらを仮定するとデータが生成される確率が高いかどうか** を判定します。
この「仮説cのもとで観察dが生じる確率の積」の事を**尤度**といい、この尤度が最大となる最も尤もらしいcを見つけるアプローチを**最尤法**と言います。

```math
P(D) = {P(c)P(d|c)} = p_c\prod_{w \in V}({P_{w,c}}^{\delta_{w,d}}(1-{P_{w,c}})^{1-\delta_{w,d}})
```

数式の記述が趣旨ではないので途中式を大分端折りますが変形すると、

```math
\log P(D) = \sum N_c \log p_c + \sum_c \sum_{w \in V} N_{w,c} \log p_{w,c} + \sum_c \sum_{w \in V} (N_c - N_{w,c}) \log(1 - p_{w,c})
```
 
こうなります。
デルタ関数がどこに行ったんだという感じですが、前述のとおり *w* = *d* の時に0、そうでなければ1となる性質があるため、単語 *w* と クラス *c* の共起回数を掛ける事で表しています。
なんだか読むのも嫌になりそうですが、ポイントは分布が**Pw,c、Pcという2つのパラメータで決定される**という点です。

分類するデータを与えた際に上記のlogP(D)を最大化するcを見つければ良いわけですが、
ここでは**世の中のツイートは全て羽生さんか羽生くんによって書かれている事を前提**としているため、
それぞれのクラスに分類される確率を足し上げると1になるという、下記の数式で表される制約を満たす必要があります。

```math
\sum_c p_c = 1
```

(ここも本題ではないので端折りますが)これは等式制約付き凸計画問題といい、
ラグランジュの未定乗数法という方法で定義したラグランジュ関数に対して、
それぞれのパラメータについて偏微分を取る事で最大値が下記の通り求められます。

```math

p_{w,c} = \frac {N_{w,c}} {N_c} , p_c = \frac {N_c} {\sum_c N_c}

```

## 分類器の実装

パラメータの求め方がわかった所で、いよいよ実装に入ります。

### データの生成

さて、形態素解析に掛けたツイートから、下記のような訓練データを生成しました。

```py

cls = ["habu", "hanyu"]

# 都合でお見せできないためイメージです。前述の通りツイートを形態素解析に掛けて生成します。
vocabulary = ["スケート", "プルシェンコ", "勝負", "神の一手"] 

# 同様
documents["habu"] = [["タイトルホルダー”， "70", "人", "半分", "羽生"],[...]]
documents["hanyu"] = [[”素晴らしい”， "4回転", "成功した", "優勝"],[...]]
```

### 同時確率 *p(w,c)* の計算

上記のデータと求めた数式から、各クラス毎に単語が生起する同時確率 *p(w,c)* を計算します。

```py
def train(cls, vocabulary, documents):

  # 各訓練文書の生起回数
  n_cls = {}
  total = 0.0
  for c in cls:
    n_cls[c] = len(documents[c])
    total += n_cls[c]

  # 各訓練文書の生起確率
  p_cls = {}
  for c in cls:
    p_cls[c] = n_cls[c] / total

  # 各クラス毎の単語の生起回数
  for c in cls:
    for d in documents[c]:
      for word in vocabulary:
        if word in d:
          n_word[c][word] += 1

　# 各クラス毎の単語の生起確率
  for c in cls:
    p_word[c] = {}
    for word in vocabulary:
      p_word[c][word] = \
        (n_word[c][word] + 1) / (n_cls[c] + 2)

```

### スムージング
ちょっと余談です。
各クラス毎の単語の生起確率を求めている部分で分子に1、分母に2を足していますが、
これは尤度が確率の積であり、語彙 *V* の単語がたまたまツイートに出現しなかった時
積算の結果が確率が0になってしまう事を防ぐためです。
(極小さな値になってしまうため実装上対数をとって和の形になっていますが、
対数の定義域にも0は存在しないのでmath domain errorでプログラムがコケます)

そのため通常は単語の出現しやすさに、0を取りにくいディリクレ分布という確率分布を仮定します。
これは最尤法によって出がちな極端な値をやわらかくする働きがあるため、**スムージング**と呼ばれています。

また、厳密なデータの出現しやすさでなく、事前分布を加味してデータが与えられた後の確率を最大化しようとするこのアプローチを**MAP推定**と言います。

# 実行結果
さて、やっと分類器が構築できたため、いよいよ実行してみます。

## 分類関数

構築した分類器を用いて、与えたツイートが羽生さんと羽生くん、
どちらによって書かれた文書であるか分類する関数です。

```py

def classify(data):
  # 各クラス毎にlogP(D)を求める
  pp = {}
  for c in cls:
    pp[c] = math.log(p_cls[c])
    for word in vocabulary:
      if word in data:
        pp[c] += math.log(p_word[c][word])

      else:
        pp[c] += math.log((1 - p_word[c][word]))

 # 求めたlogP(D)の内、どれが最も大きいか判定
 for c in cls:
   maxpp = maxpp if 'maxpp' in locals() else pp[c]
   maxcls = maxcls if 'maxcls' in locals() else c

   if maxpp < pp[c]:
     maxpp = pp[c]
     maxcls =c
                
   return (maxcls, maxpp)

```

## モデルの精度検証
取得したツイートの内、精度検証用にとっておいた30件×2人分のツイートを分類器に掛けてみます。

```py
def test(data, label):
  i = 0.0
  for tweet in data:
    if nb.classify(tweet)[0] == label:
      i += 1
  return (i / len(data))

# bags_of_wordsは各ツイートを品詞分解した２次元配列を返します
test(bags_of_words("hanyu_test.txt"), "hanyu")
test(bags_of_words("habu_test.txt"), "habu")
```

|クラス|①テストデータ件数|②正答数|正答率(②/①)|
|:-----------|------------:|------------:|------------:|
|羽生さん|30|28|93.33%|
|羽生くん|30|28|93.33%|

かなり高精度で判別できていますが、
これは羽生さんと羽生くんにそれぞれに固有の語彙が多く含まれていたからと思われます。
本来は同じ語彙で、頻度が違うような分布のデータを分類するところに意義があるため、
その点ではテストデータがあまりよくなかったかもしれませんね。

## 今後
次は[羽生さんと羽生くんの画像](https://www.google.co.jp/search?q=google+%E7%BE%BD%E7%94%9F&source=lnms&tbm=isch)の解析をやってみたいです。

## 参考文献

[言語処理のための機械学習入門](http://www.amazon.co.jp/gp/product/4339027510/ref=as_li_tf_tl?ie=UTF8&camp=247&creative=1211&creativeASIN=4339027510&linkCode=as2&tag=chordprog-22)
[羽生結弦と羽生善治の違い](http://sonicch.com/24676.html)