人工知能に関する断創録

人工知能に関連したエントリーが多いです。JavaやPythonのゲームプログラミングの記事を書いてたこともあります。



ナイーブベイズを用いたブログ記事の自動分類

カイ二乗値を用いた特徴選択(2010/6/25)の続きです。今まで使ってきた20 Newsgroupsというデータは英語文書でかつ元ネタがよく分からずあまり面白くなかったので、今回はこのブログ(人工知能に関する断想録)の記事を分類してみます。このブログの各記事には私の判断でカテゴリをつけています。たとえば、この記事は[機械学習][自然言語処理]です。カテゴリのリストはこのブログの左メニューにあります。この前、少し整理したので全部で18のカテゴリがあります。新しい記事を書いたとき自動でカテゴリを割り振ることはできるのでしょうか?

(注)プログラミング言語Pythonを使っています。シリーズもので以前作ったコードを再利用してるので検索で飛んできた人はナイーブベイズを用いたテキスト分類(2010/6/13)から順に読んでください。

はてなダイアリーデータのダウンロードと整形

まず、はてなダイアリーのデータを全部ダウンロードしてみます。これは、自分のブログなら管理>データ管理>ブログのエクスポート>Movable Type形式からダウンロードできます。Movable Type形式以外にはてな日記データ形式CSV形式でもエクスポートできますが、Movable Type形式が一番整形しやすそうでした。私のブログ(人工知能に関する断想録)のエクスポートデータは

です。けっこう長年に渡って書きためてきた感じがしますが、たった3.6MB程度とはちょっとショックかも。このファイルをナイーブベイズの入力ファイル形式に変換します。すなわち、下のような形式のファイルです。

カテゴリ 単語:カウント 単語:カウント ...     <- 文書1
カテゴリ 単語:カウント 単語:カウント ...     <- 文書2
...

整形のポイントは、

  • Movable Typeのファイルからエントリごとにタイトル、カテゴリ、日付、本文を抽出する
  • ソースコード部分は意味のある単語を含まないので無視する
  • タイトルと本文から単語(名詞と固有名詞のみ)を形態素解析で抽出する
  • 形態素解析にはMeCab Pythonを用いる
  • エントリごとに単語の出現頻度をカウントする
  • エントリごとにカテゴリと単語:カウントをファイルに出力する
  • 単語はカウントの降順でソートする

参考までに変換プログラムを掲載します。この変換プログラムをaidiary.txtと同じフォルダでプログラムを実行すると

という変換されたファイルができます。ただし、このプログラムを動かすには後述する形態素解析モジュールMeCab Pythonの導入が必要です

#coding:utf-8
import codecs
import re
import sys
from collections import defaultdict

# MeCab Pythonのインストールが必要
import MeCab

"""
はてなダイアリーのブログデータの前処理
HTMLタグは除去
プログラムコードブロックは除去
"""

# ブログデータ(Movable Type形式)
# 自分のブログなら管理メニューから全文ダウンロード可
FILE = "aidiary.txt"
OUTPUT = "blog.txt"

def analysis(tagger, text):
    """textをMeCabで形態素解析して単語のリストを返す"""
    node = tagger.parseToNode(text.encode("utf-8"))
    wordList = []
    while node:
        temp = node.feature.split(",")
        if temp[0] == u"名詞":
            # 一般名詞、固有名詞のみ抽出する
            if temp[1] == u"一般" or temp[1] == u"固有名詞":
                wordList.append("%s" % node.surface)
        node = node.next
    return wordList

if __name__ == "__main__":
    # 形態素解析器
    tagger = MeCab.Tagger("-Ochasen")
    
    # 属性抽出用の正規表現
    htmlTagPattern = re.compile(r'<.*?>')
    titlePattern = re.compile(r'TITLE:\s(.*)')
    categoryPattern = re.compile('CATEGORY:\s(.*)')
    datePattern = re.compile(r'DATE:\s(.*)')
    bodyPattern = re.compile(r'BODY:')
    codeStartPattern = re.compile(r'<pre class="syntax-highlight">')
    codeEndPattern = re.compile(r'</pre>')
    
    fin = codecs.open(FILE, "r", "utf-8")
    fout = codecs.open(OUTPUT, "w", "utf-8")
    
    title = category = date = body = ""
    bodyflg = codeflg = False
    for line in fin:
        line = line.rstrip()           # 末尾改行除去
        line = line.replace("\t", "")  # タブ文字除去
        if line == "": continue
        
        # タイトルの抽出
        m = titlePattern.match(line)
        if m:
            title = m.group(1)
        
        # カテゴリの抽出(2つ以上ある場合は1つ目のみ)
        m = categoryPattern.match(line)
        if m and category == "":
            category = m.group(1)
        
        # 日付の抽出
        m = datePattern.match(line)
        if m:
            date = m.group(1)          # 02/21/2002 05:30:47 PM
            date = date.split(" ")[0]  # 02/21/2002
            month, day, year = [int(x) for x in date.split("/")[0:3]]
            date = "%04d-%02d-%02d" % (year, month, day)  # 2002-02-21
        
        # プログラムコード行が見つかったらcodeflgを立てる
        m = codeStartPattern.match(line)
        if m:
            codeflg = True
        
        # プログラムコード終了行が見つかったらcodeflgを下ろす
        m = codeEndPattern.match(line)
        if m:
            codeflg = False
        
        # bodyflg = Trueでコード行以外なら本文行なのでbodyに追加
        if bodyflg and not codeflg and line != "-----":
            # HTMLタグを除去して追加
            body = body + ' ' + htmlTagPattern.sub("", line)
        
        # 本文終了記号が見つかったら1エントリ分のデータを出力
        if bodyflg and line == "-----":
            # タイトルと本文を形態素解析して単語を抽出
            wordList = analysis(tagger, title + body)
            # 単語:カウントの形式で出力するために単語をカウント
            wc = defaultdict(int)
            for w in wordList:
                wc[w] += 1
            # カウントで降順ソート
            lst = wc.items()
            lst.sort(lambda p0,p1: cmp(p1[1],p0[1]))
            # ファイルに出力
            words = []
            for k, v in lst:
                words.append("%s:%d" % (k, v))
            print category, " ".join(words)
            fout.write("%s %s\n" % (category, " ".join(words)))
            bodyflg = False
            title = category = date = body = ""
        
        # 本文が見つかったらフラグを立てる(次の行からbodyに追加)
        m = bodyPattern.match(line)
        if m:
            bodyflg = True
    
    fout.close()
    fin.close()

MeCab Python

形態素解析にはMeCabのインストールが必要です。MeCabは高性能な形態素解析モジュールでPythonRubyPerlJavaなどさまざまな言語から使えます。Mac OS XLinuxでは簡単にコンパイルしてインストールができるんですが、WindowsではMinGWVisual Studioのインストール & コードの修正が必要でかなり面倒くさい。そこで、Pythonモジュールはid:fgshunさんがコンパイルしたバイナリを使わせてもらいました。以下、導入方法です。

  1. MeCabの本サイトでダウンロードしたWindows版のmecab-0.98.exeをインストール(辞書はUTF-8形式が無難です)
  2. 形態素解析エンジン MeCab 0.98pre3 野良ビルドからダウンロードしたlibmecab-1.dll、MeCab.py、_MeCab.pydをパッケージフォルダ(Python2.5ならC:\Python25\Lib\site-packages)にコピー。IPA辞書はmecab-0.98.exeでインストールしたので不要。

id:fgshunさんはTaggerにmecabrcファイルを指定してますが、指定しないとデフォルトでC:\Program Files\MeCab\etc\mecabrcを読みにいくようなのでmecab-0.98.exeでIPA辞書をインストールしたほうが簡単だと思います。下のサンプルプログラムが動くか確かめてみます。

#coding:utf-8
import MeCab
s = u"形態素解析にはMeCabを使っているので上のプログラムを動かすには事前にインストールが必要です"
tagger = MeCab.Tagger('-Ochasen')
node = tagger.parseToNode(s.encode('utf-8'))
wordList = []
while node:
    if node.feature.split(",")[0] == u"名詞":
        wordList.append(node.surface)
    node = node.next
for word in wordList:
    print word

実行すると名詞のみが抽出されて出力されます。

形態素
解析
MeCab
上
プログラム
事前
インストール
必要

MeCabはまた使いそうなのであとでまとめとこうかな。形態素解析にはMeCabではなく、Yahoo!形態素解析API(2009/4/15)を用いてもよいと思います。Yahoo!はMeCabより固有表現の抽出精度が高いように思います。

カテゴリ分布

では、話はナイーブベイズに戻ります。まず、カテゴリごとに各記事がどれくらいあるか統計を取ってみました。記事に複数のカテゴリがある場合は第一カテゴリのみ使っています。カテゴリ分布は下のようになりました。20 Newsgroups(2010/6/18)と違ってかなりばらつきがあります。一般的にカテゴリ間のばらつきが大きいと精度が低下することが知られています。ちょっと心配になってきた。

f:id:aidiary:20100703210346p:plain

特徴語抽出

次に各カテゴリで相互情報量(2010/6/19)が高い上位10個の単語を抽出してみます。

[人工知能]
0.0668466792518 知能
0.0605570522974 人工
0.032089567665 学会
0.0254121860825 時代
0.0253081836479 AI
0.0185970342997 技術
0.0182707784424 Norvig
0.0178956959728 機械
0.0160822792793 分野
0.013954708507 康一

[機械学習]
0.0554374306209 ベイズ
0.0487473103335 データ
0.0463269427927 式
0.0431790653499 最小
0.0431790653499 w
0.0430882201539 ベクトル
0.0362336163557 線形
0.0357786049225 パターン
0.0295890551301 誤差
0.0276829622679 t

[自然言語処理]
0.0407276053721 形態素
0.0305270697973 Google
0.0209781544664 マルコフ
0.0209781544664 コンテンツ
0.0198496993147 名詞
0.0198496993147 スパイダー
0.0176073150915 Web
0.0174213092901 連鎖
0.0174213092901 Hacks
0.0174213092901 Eliza

[ロボティクス]
0.0836949238735 ロボット
0.038853580853 AIBO
0.0300575326888 センサー
0.0278687435843 ROBODEX
0.0225409242532 ヒューマノイド
0.0173236338447 ロボカップ
0.0173236338447 インテリジェンス・ダイナミクス
0.0158608128224 日記
0.0149268606495 QRIO
0.0125190099368 産業

[哲学]
0.0696987261079 哲学
0.0466277158515 ゲーデル
0.0433229397643 バッハ
0.0433229397643 エッシャー
0.0402016652066 ホフスタッター
0.0293236094059 ドレイファス
0.0293236094059 ダグラス
0.0278300807178 コンピュータ
0.0230329975586 定理
0.0222695028808 揚

[認知科学]
0.0344178968077 心理
0.0291007471038 動物
0.0215349075037 科学
0.0158318961015 身体
0.0148852974513 自身
0.0140415540977 ヒト
0.0131966937095 相互
0.0129084130782 犬
0.0126732630125 ゲーム
0.0124394410983 報酬

[強化学習]
0.0821109858052 報酬
0.072511052999 エージェント
0.0702965495728 迷路
0.0470118931038 Reinforcement
0.0469046359785 Q
0.0381829229148 振子
0.0380810498439 状態
0.0371525790713 エピソード
0.0325025838501 ブログ
0.0309618292129 テーブル

[複雑系]
0.047560226878 乱数
0.0422215036047 現象
0.0407280367457 ランダム
0.0303254834174 擬似
0.0303254834174 スモール
0.025055107977 ワールド
0.0232314372039 カオス
0.0189085366345 最前線
0.0189085366345 エントロピー
0.0183085362125 確率

[脳科学]
0.114387617974 脳
0.0314007895288 コンピューター
0.0201049064124 患者
0.0196755440911 大脳
0.0161414204785 茂木
0.0143327777787 科学
0.0133952718291 健一郎
0.0124255845359 つながり
0.0122488069338 海馬
0.0108136899383 筑摩書房

[生物学]
0.0311204925831 バイオインフォマティクス
0.0304531477327 細胞
0.025867712687 生物
0.0160213110636 動物
0.0133887567948 生命
0.0130598253523 見方
0.0130598253523 個体
0.0104193172884 遺伝子
0.00996531603537 生き物
0.00927583404748 ヒト

おおおおお。地味だけどけっこうすごいぞ。

交差検定(Cross Validation)

最後にナイーブベイズの分類精度を評価してみます。20 Newsgroupsは文書数が20000件と多かったので訓練データとテストデータを完全に分けてから分類精度を評価しました。しかし、ブログ記事は全部で760件しかないので独立したテストデータを用意するのは少しもったいないです。そこで、交差検定(N-fold Cross Validation)という評価法を使ってみます。N-foldの代わりにK-foldを使うことが多いようですが、すでにボキャブラリの数で文字Kを使ってしまったので代わりにNを使いました。N-fold Cross Validationは全データを均等にN分割し、その中から1個をテストデータとし、残りのN-1個を訓練データとする評価法です。テストデータの選び方はN通りあるためN通りの訓練データとテストデータのそれぞれで分類精度を求め、最後に平均を取って最終精度とすることが多いです。ちなみに分割数Nをデータ数と同じにした場合はLeave-one-outと呼びます。

では、Cross Validationのコードです。Recipe 521906: K fold cross validation partition (Python)のコードを参考にしました。剰余を使ってデータを分けるのはけっこううまいやり方だなぁ。下のコードを動かすには相互情報量を用いた特徴選択(2010/6/19)で作ったナイーブベイズ(naivebayes_fs.py)と特徴選択(feature_selection.py)が必要です。ボキャブラリの数は適当に1000としました(相互情報量が高い順に1000語を使います)。

#coding:utf-8
import codecs
import sys
from naivebayes_fs import NaiveBayes  # 特徴選択付きのナイーブベイズを使う
from pylab import *

# 特徴抽出に使う単語数
K = 1000

def crossValidation(data, N=5, randomize=False):
    """N-fold Cross Validationで分類精度を評価"""
    # データをシャッフル
    if randomize:
        from random import shuffle
        shuffle(data)
    
    # N-fold Cross Validationで分類精度を評価
    accuracyList = []
    for n in range(N):  # 各分割について
        # 訓練データとテストデータにわける
        trainData = [d for i, d in enumerate(data) if i % N != n]
        testData = [d for i, d in enumerate(data) if i % N == n]
        # ナイーブベイズ分類器を学習
        nb = NaiveBayes(K)
        nb.train(trainData)
        print nb
        # テストデータの分類精度を計算
        hit = 0
        numTest = 0
        for d in testData:
            correct = d[0]
            words = d[1:]
            predict = nb.classify(words)
            if correct == predict:
                hit += 1
            numTest += 1
        accuracy = float(hit) / float(numTest)
        accuracyList.append(accuracy)
    # N回の平均精度を求める
    print accuracyList
    average = sum(accuracyList) / float(N)
    return average

if __name__ == "__main__":
    # ブログデータをロード
    data = []
    fp = codecs.open("blog.txt", "r", "utf-8")
    for line in fp:
        line = line.rstrip()
        temp = line.split()
        data.append(temp)
    fp.close()
    # N-fold Cross Validationで分類精度を評価
    average = crossValidation(data, N=10, randomize=True)
    print "accuracy:", average

実行すると

accuracy: 0.651315789474

20 Newsgroupsは80%だったので少し落ちますけど思っていたより高いなぁ。

まとめ

今回は、日本語のブログ記事のデータをナイーブベイズで分類してみました。実データらしくデータ数の偏りが大きかったですがそこそこな分類精度が出ました。今回はボキャブラリ数とか抽出する形態素の品詞選択とかかなり適当なので本来ならもっとちゃんとした検証が必要です。とは言いながら、次回はPRMLで学んだ(2010/5/3)サポートベクトルマシン(SVM)を使ってブログ記事を分類してみます。SVMはテキスト分類で高精度を叩き出すことで有名な手法です。どれくらい精度が上がるか楽しみだなぁ。

カイ二乗値を用いた特徴選択

相互情報量を用いた特徴選択(2010/6/19)のつづきです。今回は、相互情報量ではなく、カイ二乗を用いて特徴語を抽出してみます。カイ二乗検定は独立性の検定によく使いますけど、特徴語の抽出にも応用できるってのははじめて知りました。結局のところ相互情報量カイ二乗値もカテゴリと単語がどれくらい依存しているかを表す尺度なのでアプローチは似ている感じがします。IIR13.5を参考にして実装します。

カイ二乗

カイ二乗値の定義は、

f:id:aidiary:20100625000913p:plain

です。NやEが出てきますが、下のようなクロス表を用いて計算します。たとえば、単語「iPhone」とカテゴリ「IT」のカイ二乗値を求めたいとき、クロス表は下のようになります。たとえば、カテゴリがITで単語iPhoneを含む文書はデータ中にN11個あるなどと解釈します。

カテゴリがITである カテゴリがITでない
単語iPhoneを含む N11 (E11) N10 (E10) N1.
単語iPhoneを含まない N01 (E01) N00 (E00) N0.
N.1 N.0 N

N.1はN01+N11、N1.はN10+N11の略です。ここで、カテゴリITと単語iPhoneが独立であると仮定したときに期待される文書数がEです。2つの事象が独立であるとき上のようなクロス表上では各行(または各列)の出現頻度の比が等しくなります。なので、iPhoneを含むケースでも含まないケースでもカテゴリがITである割合はN.1/Nで、カテゴリがITでない割合はN.0/Nとなります。これはプログラミングのための確率統計でもくどいくらい説明されてました。こんなところで役立つとは!。

以上の結果から具体的な文書数E..は下の式で計算できます。

f:id:aidiary:20100625004119p:plain

で、さっきのカイ二乗値の式は実際の文書数Nijと単語とカテゴリが独立と仮定したときの期待される文書数Eijがどれだけ乖離しているかを表す指標になっています。この乖離が大きいほど単語とカテゴリの依存度合いが強いと解釈できるとのこと。相互情報量のときと同じくカイ二乗値が大きい単語ほどカテゴリに特徴的な単語と言ってよさそうです。

カイ二乗値が高い単語を抽出

上の定義をそのまま素直にPythonで実装してみます。20 Newsgroupsの各カテゴリからカイ二乗値が高い上位k個の単語を求めるのが目標です。クロス表を求めるところまでは相互情報量のときとまったく同じです。

#coding:utf-8
import codecs
import math
import sys
from collections import defaultdict

# feature_selection.py

def chisquare(target, data, k=0):
    """カテゴリtargetにおけるカイ二乗値が高い上位k件の単語を返す"""
    # 上位k件を指定しないときはすべて返す
    if k == 0: k = sys.maxint
    
    V = set()
    N11 = defaultdict(float)  # N11[word] -> wordを含むtargetの文書数
    N10 = defaultdict(float)  # N10[word] -> wordを含むtarget以外の文書数
    N01 = defaultdict(float)  # N01[word] -> wordを含まないtargetの文書数
    N00 = defaultdict(float)  # N00[word] -> wordを含まないtarget以外の文書数
    Np = 0.0  # targetの文書数
    Nn = 0.0  # target以外の文書す
    
    # N11とN10をカウント
    for d in data:
        cat, words = d[0], d[1:]
        if cat == target:
            Np += 1
            for wc in words:
                word, count = wc.split(":")
                V.add(word)
                N11[word] += 1  # 文書数をカウントするので+1すればOK
        elif cat != target:
            Nn += 1
            for wc in words:
                word, count = wc.split(":")
                V.add(word)
                N10[word] += 1
    
    # N01とN00は簡単に求められる
    for word in V:
        N01[word] = Np - N11[word]
        N00[word] = Nn - N10[word]
    # 総文書数
    N = Np + Nn
    
    # 各単語のカイ二乗値を計算
    chiSquare = []
    for word in V:
        n11, n10, n01, n00 = N11[word], N10[word], N01[word], N00[word]
        
        # カテゴリと単語が独立である(帰無仮説)と仮定したときの期待値を求める
        e11 = (n10 + n11) * (n01 + n11) / float(N)
        e10 = (n10 + n11) * (n00 + n10) / float(N)
        e01 = (n00 + n01) * (n01 + n11) / float(N)
        e00 = (n00 + n01) * (n00 + n10) / float(N)
        
        # カイ二乗値の各項を計算
        score = (n00-e00)**2/e00 + (n01-e01)**2/e01 + (n10-e10)**2/e10 + (n11-e11)**2/e11
        chiSquare.append( (score, word) )
    
    # カイ二乗値の降順にソートして上位k個を返す
    chiSquare.sort(reverse=True)
    return chiSquare[0:k]

if __name__ == "__main__":
    # 訓練データをロード
    trainData = []
    fp = codecs.open("news20", "r", "utf-8")
    for line in fp:
        line = line.rstrip()
        temp = line.split()
        trainData.append(temp)
    fp.close()
    
    # カイ二乗値を用いて特徴選択
    target = "comp.graphics"
    features = chisquare(target, trainData, k=10)
    print "[%s]" % target
    for score, word in features:
        print score, word

行すると、comp.graphicsカテゴリの相互情報量が高い順に10件の単語が出力されます。

[comp.graphics]
1457.58534448 graphics
612.848781005 image
607.658961581 animation
534.398607189 polygon
483.393632538 gif
389.517701949 tiff
367.564520621 pov
365.271144653 images
349.145975364 cview
340.23014597 jpeg

いかにもcomp.graphicsっぽい単語が並んでます。相互情報量と少し違った単語が上位に上がってます。どちらがいいんでしょうか?

[comp.os.ms-windows.misc]
3254.01838505 windows
863.806249156 cica
783.368406974 dos
619.487823399 ini
517.674741517 file
468.921312024 microsoft
462.793325495 win
438.326468403 ms
422.676491104 drivers
385.311002551 files

[rec.sport.baseball]
2242.51216749 baseball
1482.94560718 pitching
1028.09080232 mets
936.43915033 hitter
910.661498564 braves
882.2585267 phillies
851.746184022 season
848.192467208 pitcher
775.729371181 pitchers
750.249334866 cubs

[sci.space]
2291.74284858 space
1941.32877456 orbit
1528.74443512 nasa
1403.56063934 launch
1302.97108661 moon
1257.05550119 shuttle
1157.29780773 lunar
993.562973086 spacecraft
922.346259941 nsmca
827.767660218 flight

[talk.politics.guns]
2759.33818739 gun
2135.88183801 guns
1865.90273163 firearms
1347.63548737 weapons
1087.54646142 firearm
1025.85886288 handgun
1023.15971587 batf
851.875428063 weapon
800.696565647 waco
794.415883505 nra

[talk.religion.misc]
774.507903866 sandvik
593.936874472 jesus
590.599644857 kendig
576.386454526 newton
569.000025676 kent
498.495688011 ksand
486.967079804 alink
477.218439237 bskendig
473.722409841 christian
464.190984532 royalroads

特徴選択を考慮したナイーブベイズ

では、本題のナイーブベイズで特徴選択を考慮した分類器を作成できるようにします。news20、news20.tのデータはそのままでナイーブベイズのプログラム側でカイ二乗値が大きい上位k個のボキャブラリのみ使うようにします。NaiveBayesクラスのコンストラクタにkの値を指定します。相互情報量のナイーブベイズプログラムで

from feature_selection import mutual_information
features = mutual_information(cat, data)

の部分を

from feature_selection import chisquare
features = chisquare(cat, data)

と書き換えるだけでOKです。ではまた横軸にボキャブラリサイズの対数スケール、縦軸に分類精度としてグラフを描いてみます。

f:id:aidiary:20100625214420p:plain

相互情報量のときとあまり変わらない・・・ただこれもIIRに載ってますが、Reuters-RCV1というデータセットだと相互情報量カイ二乗値では差が出てきます。相互情報量の方が単語数が少なくても分類精度が高い。とりあえず両方試してみるのがいいのかも。

今回はなじみのない英語の文書だったので、次回からは日本語のブログ記事(このブログ)のデータを分類してみようと思います。

参考文献

相互情報量を用いた特徴選択

20 Newsgroupsで分類精度を評価(2010/6/18)のつづきです。今回は、特徴選択に挑戦してみようと思います。テキスト分類における特徴とは基本的に単語のことです。

特徴選択

前回、ナイーブベイズの出力結果で

documents: 11269, vocabularies: 53852, categories: 20
accuracy: 0.802265156562

となってました。documentsは訓練データの総文書数、categoriesは訓練データのカテゴリ数、vocabulariesは訓練データの総単語数を表します。テキスト分類において53852個の単語を考慮していることを意味します。しかし、この単語の中には分類に寄与しないばかりかノイズになって逆に性能を悪化させるような単語が含まれていることがあります。たとえば、the, in, toなどのストップワードがその一例です。その他にもすべてのカテゴリに同じくらいの頻度で出現する単語なんかもそうです。そのような単語は分類の役に立ちません。今回は、53852個のボキャブラリをさらに絞り込む(特徴を選択する)のが目標です。

相互情報量

特徴選択で代表的なのは、相互情報量(Mutual Information)という尺度を用いる手法です。IIR13.5を参考にして実装します。

相互情報量は2つの確率変数の相互依存の尺度を表す量とのこと。テキスト分類の特徴選択で用いる場合は、ある単語tの出現を表す確率変数Uとあるカテゴリcの出現を表す確率変数Cを用いて相互情報量I(U; C)を定義します。Uは1または0の値をとり、U=1のとき単語tが出現する事象、U=0のとき単語tが出現しないという事象を表します。Cも1または0の値をとり、C=1のときカテゴリがcである、C=0のときカテゴリがcでないという事象を表します。相互情報量の定義は、


I(U;C) = \displaystyle \sum_{e_t \in \{1,0\}} \displaystyle \sum_{e_c \in \{1,0\}} P(U = e_t, C = e_c) \log_2 \frac{P(U = e_t, C = e_c)}{P(U = e_t) P(C = e_c)}

です。tはterm、cはcategoryの略で具体的な単語やカテゴリが入ります。同時分布や周辺分布が出てきますが、下のようなクロス表を用いると簡単に計算できます。たとえば、単語「iPhone」とカテゴリ「IT」の相互情報量を求めたいとき、クロス表は下のようになります。

カテゴリがITである カテゴリがITでない
単語iPhoneを含む N11 N10
単語iPhoneを含まない N01 N00

ここで、

  • N11は、訓練文書中でカテゴリがIT(C=1)でかつ単語iPhoneを含む(U=1)文書数
  • N10は、訓練文書中でカテゴリがITでなく(C=0)かつ単語iPhoneを含む(U=1)文書数
  • N01は、訓練文書中でカテゴリがIT(C=1)でかつ単語iPhoneを含まない(U=0)文書数
  • N00は、訓練文書中でカテゴリがITでなく(C=0)かつ単語iPhoneを含まない(U=0)文書数

となります。このようなクロス表が求まれば、

  • P(U=1, C=1) = N11 / N
  • P(U=1) = (N10+N11) / N
  • P(C=1) = (N01+N11) / N

といった感じにすべての同時確率と周辺確率が簡単に計算できます。I(U; C)をNを使って書き直すと

f:id:aidiary:20100619120136p:plain

となります。ここで、N1.はN10+N11の略です。上のようなクロス表は、単語とカテゴリのすべての組み合わせ分だけ生成されます。たとえば、単語が50000個あってカテゴリが20個あれば、50000x20通りのクロス表ができます。相互情報量は0以上の値をとり、値が大きいほどカテゴリの特徴を表すような単語と見なすことができるとのこと。これがなぜかは少し考えてしまったけど、いくつか極端なクロス表を作って相互情報量を計算してみると納得できるかも。単語とカテゴリが完全に独立だと相互情報量は0になってしまうってことは・・・そのような単語はカテゴリの内容を表すとは言いがたいってのが直感的な理解でしょうか?

相互情報量が高い単語を抽出

上の定義をそのまま素直にPythonで実装してみます。20 Newsgroupsの各カテゴリから相互情報量が高い上位k個の単語を求めるのが目標です。

#coding:utf-8
import codecs
import math
import sys
from collections import defaultdict

# feature_selection.py

def mutual_information(target, data, k=0):
    """カテゴリtargetにおける相互情報量が高い上位k件の単語を返す"""
    # 上位k件を指定しないときはすべて返す
    if k == 0: k = sys.maxint
    
    V = set()
    N11 = defaultdict(float)  # N11[word] -> wordを含むtargetの文書数
    N10 = defaultdict(float)  # N10[word] -> wordを含むtarget以外の文書数
    N01 = defaultdict(float)  # N01[word] -> wordを含まないtargetの文書数
    N00 = defaultdict(float)  # N00[word] -> wordを含まないtarget以外の文書数
    Np = 0.0  # targetの文書数
    Nn = 0.0  # target以外の文書す
    
    # N11とN10をカウント
    for d in data:
        cat, words = d[0], d[1:]
        if cat == target:
            Np += 1
            for wc in words:
                word, count = wc.split(":")
                V.add(word)
                N11[word] += 1  # 文書数をカウントするので+1すればOK
        elif cat != target:
            Nn += 1
            for wc in words:
                word, count = wc.split(":")
                V.add(word)
                N10[word] += 1
    
    # N01とN00は簡単に求められる
    for word in V:
        N01[word] = Np - N11[word]
        N00[word] = Nn - N10[word]
    # 総文書数
    N = Np + Nn
    
    # 各単語の相互情報量を計算
    MI = []
    for word in V:
        n11, n10, n01, n00 = N11[word], N10[word], N01[word], N00[word]
        # いずれかの出現頻度が0.0となる単語はlog2(0)となってしまうのでスコア0とする
        if n11 == 0.0 or n10 == 0.0 or n01 == 0.0 or n00 == 0.0:
            MI.append( (0.0, word) )
            continue
        # 相互情報量の定義の各項を計算
        temp1 = n11/N * math.log((N*n11)/((n10+n11)*(n01+n11)), 2)
        temp2 = n01/N * math.log((N*n01)/((n00+n01)*(n01+n11)), 2)
        temp3 = n10/N * math.log((N*n10)/((n10+n11)*(n00+n10)), 2)
        temp4 = n00/N * math.log((N*n00)/((n00+n01)*(n00+n10)), 2)
        score = temp1 + temp2 + temp3 + temp4
        MI.append( (score, word) )
    
    # 相互情報量の降順にソートして上位k個を返す
    MI.sort(reverse=True)
    return MI[0:k]

if __name__ == "__main__":
    # 訓練データをロード
    trainData = []
    fp = codecs.open("news20", "r", "utf-8")
    for line in fp:
        line = line.rstrip()
        temp = line.split()
        trainData.append(temp)
    fp.close()
    
    # 相互情報量を用いて特徴選択
    target = "comp.graphics"
    features = mutual_information(target, trainData, k=10)
    print "[%s]" % target
    for score, word in features:
        print score, word

実行すると、comp.graphicsカテゴリの相互情報量が高い順に10件の単語が出力されます。

[comp.graphics]
0.0396591914417 graphics
0.018190625901 image
0.013256620866 animation
0.0122108176792 gif
0.0109647921474 polygon
0.0102937113306 images
0.00984347501837 files
0.00839170021167 format
0.00832240938732 tiff
0.00799473476391 people

いかにもcomp.graphicsっぽい単語が並んでます。他にもいくつかのカテゴリで試してみます。

[comp.os.ms-windows.misc]
0.094183661915 windows
0.0237696712358 dos
0.0184571833467 file
0.0176200422274 cica
0.0161838846441 win
0.014151739062 ms
0.0135110433857 files
0.0131853055706 drivers
0.0126673382015 driver
0.0126148229575 ini

[rec.sport.baseball]
0.0522161980385 baseball
0.031323260157 pitching
0.0241678683103 season
0.0226248072068 games
0.021565415199 mets
0.0204313593006 team
0.0198532689063 braves
0.0195487548191 hitter
0.0192974411302 game
0.0184112914927 phillies

[sci.space]
0.0648087600801 space
0.0414339611904 nasa
0.0410198092698 orbit
0.0307253416988 launch
0.0290759088133 moon
0.0266171725319 shuttle
0.0242961005478 lunar
0.0224255674403 earth
0.0208045078764 spacecraft
0.0197463695696 flight

[talk.politics.guns]
0.064023268956 gun
0.047075696447 guns
0.0372513306097 firearms
0.0321421838979 weapons
0.0213623165189 batf
0.0201845527411 handgun
0.019986956496 fire
0.018426078685 weapon
0.018311826034 fbi
0.0182303153541 waco

[talk.religion.misc]
0.0167788239223 jesus
0.0152625625906 god
0.014476187063 christian
0.013821963122 sandvik
0.0122269234191 kent
0.0121410088867 bible
0.0116866079726 christians
0.0115979054194 newton
0.010151710003 religion
0.00988892270487 christ

それっぽい、それっぽい。相互情報量でそのカテゴリの特徴語が抽出できるのはけっこう面白い。

特徴選択を考慮したナイーブベイズ

では、本題のナイーブベイズで特徴選択を考慮した分類器を作成できるようにします。news20、news20.tのデータはそのままでナイーブベイズのプログラム側で相互情報量が大きい上位k個のボキャブラリのみ使うようにします。NaiveBayesクラスのコンストラクタにkの値を指定します。

#coding:utf-8
import math
import sys
from collections import defaultdict
from feature_selection import mutual_information

# 特徴選択を行うナイーブベイズ分類器

class NaiveBayes:
    """Multinomial Naive Bayes"""
    def __init__(self, k):          # 相互情報量が大きい順にk個の単語をボキャブラリとする
        self.categories = set()     # カテゴリの集合
        self.vocabularies = set()   # ボキャブラリの集合
        self.wordcount = {}         # wordcount[cat][word] カテゴリでの単語の出現回数
        self.catcount = {}          # catcount[cat] カテゴリの出現回数
        self.denominator = {}       # denominator[cat] P(word|cat)の分母の値
        self.k = k                  # ボキャブラリ数
    
    def train(self, data):
        """ナイーブベイズ分類器の訓練"""
        # 文書集合からカテゴリを抽出して辞書を初期化
        for d in data:
            cat = d[0]
            self.categories.add(cat)
        for cat in self.categories:
            self.wordcount[cat] = defaultdict(int)
            self.catcount[cat] = 0
        
        # 特徴選択してボキャブラリを絞り込む
        L = []
        for cat in self.categories:
            features = mutual_information(cat, data)
            L.extend(features)
        L.sort(reverse=True)
        for i in range(len(L)):
            # L[i]=(score, word)なので単語はL[i][1]で取り出せる
            self.vocabularies.add(L[i][1])
            # ボキャブラリの数が指定した数に達したら終了
            if len(self.vocabularies) == self.k: break
        
        # 文書集合からカテゴリと単語をカウント
        for d in data:
            cat, doc = d[0], d[1:]
            self.catcount[cat] += 1
            for wc in doc:
                word, count = wc.split(":")
                count = int(count)
                # 単語がボキャブラリに含まれなければ無視
                if not word in self.vocabularies: continue
                self.wordcount[cat][word] += count
        # 単語の条件付き確率の分母の値をあらかじめ一括計算しておく(高速化のため)
        for cat in self.categories:
            self.denominator[cat] = sum(self.wordcount[cat].values()) + len(self.vocabularies)
    
    def classify(self, doc):
        """事後確率の対数 log(P(cat|doc)) がもっとも大きなカテゴリを返す"""
        best = None
        max = -sys.maxint
        for cat in self.catcount.keys():
            p = self.score(doc, cat)
            if p > max:
                max = p
                best = cat
        return best
    
    def wordProb(self, word, cat):
        """単語の条件付き確率 P(word|cat) を求める"""
        # ラプラススムージングを適用
        # 分母はtrain()の最後で一括計算済み
        return float(self.wordcount[cat][word] + 1) / float(self.denominator[cat])
    
    def score(self, doc, cat):
        """文書が与えられたときのカテゴリの事後確率の対数 log(P(cat|doc)) を求める"""
        total = sum(self.catcount.values())  # 総文書数
        score = math.log(float(self.catcount[cat]) / total)  # log P(cat)
        for wc in doc:
            word, count = wc.split(":")
            count = int(count)
            # 単語がボキャブラリに含まれなければ無視
            if not word in self.vocabularies: continue
            # logをとるとかけ算は足し算になる
            for i in range(count):
                score += math.log(self.wordProb(word, cat))  # log P(word|cat)
        return score
    
    def __str__(self):
        total = sum(self.catcount.values())  # 総文書数
        return "documents: %d, vocabularies: %d, categories: %d" % (total, len(self.vocabularies), len(self.categories))

横軸にボキャブラリサイズの対数スケール、縦軸に分類精度としてグラフを描いてみました。

f:id:aidiary:20100619212419p:plain

あれ?ボキャブラリサイズを減らすと精度は下がってしまいました・・・実はこういうこともあるみたいです。元ネタの論文

でも同じ結果になっています。ただ、データセットによってはボキャブラリサイズを減らすと精度が上がる場合もあるようで、IIRにはReuters-RCV1というデータの実行例が載っています。オリジナルのボキャブラリサイズが132776ですが、相互情報量が大きい順に100個に絞るとaccuracyが20%近く上がっています。こんなお得なこともあるのでとりあえずやってみた方がよいのかも。

20 Newsgroupsで分類精度を評価

ナイーブベイズを用いたテキスト分類(2010/6/13)の続きです。前回、実装したナイーブベイズの分類精度を評価してみます。テキスト分類のベンチマークとして使われるのは

といったデータセットです。今回は、ナイーブベイズの分類精度を20 Newsgroupsで評価してみたいと思います。論文は散々読んだけど自分で試すのは初めてなんだよなー。

20 Newsgroups

http://qwone.com/~jason/20Newsgroups/

Usenet*1から収集した約20000文書、20カテゴリのデータセットです。カテゴリは下の20個。まあ何となくどんなカテゴリなのかわかりますね。おおまかにcomp、rec、sci、talkに分けられるので4カテゴリとして扱うこともあるようです。

  1. comp.graphics
  2. comp.os.ms-windows.misc
  3. comp.sys.ibm.pc.hardware
  4. comp.sys.mac.hardware
  5. comp.windows.x
  6. rec.autos
  7. rec.motorcycles
  8. rec.sport.baseball
  9. rec.sport.hockey
  10. sci.crypt
  11. sci.electronics
  12. sci.med
  13. sci.space
  14. talk.politics.misc
  15. talk.politics.guns
  16. talk.politics.mideast
  17. talk.religion.misc
  18. alt.atheism
  19. misc.forsale
  20. soc.religion.christian

上のサイトに行くと英文記事のオリジナルデータもダウンロードできますが、親切なことにbag-of-wordsの形にしたデータも提供されています。今回はMatlab/Octaveで簡単に読み込める形式になっている20news-bydate-matlab.tgzというデータを使います。解凍すると下の6つのファイルができます。あとvocabulary.txtもいっしょにダウンロードしておきます。

  • train.data - 訓練データ集合、1行は「文書ID」「単語ID」「頻度」の3つ組。
  • train.label - 訓練データの各文書のラベル(カテゴリID)
  • train.map - 訓練データのカテゴリIDとカテゴリ名の対応
  • test.data - テストデータ集合、フォーマットはtrain.dataと同じ
  • test.label - テストデータの各文書のラベル(カテゴリID)
  • test.map - テストデータのカテゴリIDとカテゴリ名の対応(train.mapと同じ)
  • vocabulary.txt - 単語IDと単語の対応

テキスト分類の性能を評価するためには訓練データとテストデータにわけるのが一般的です。訓練データで分類器を学習して、テストデータのラベルを分類器で予測して精度を評価します。訓練データで学習して訓練データのカテゴリを予測するのはNGです(ものすごく精度が高くなる)。まず、簡単な統計をとってみます。train.labelとtest.labelを使って、訓練データ集合とテストデータ集合に各カテゴリの記事がいくつあるか調べてみました。結果は(画像をクリックしてオリジナルサイズを表示で拡大できます)

f:id:aidiary:20100617220155p:plain

なるほど。各カテゴリは約1000の記事からなっていて訓練とテストで6:4くらいでわけているようです。訓練データとテストデータは綺麗に同じ形をしています。ナイーブベイズは訓練データとテストデータが同じ分布から生成されたという仮定があるようですが、実際の応用ではこんなに綺麗に同じ比率で含まれることはまれだと思います(あとでブログ記事分類を取り上げます)。

ちなみに上のグラフを書くPythonプログラムです。グラフ描画のライブラリmatplotlibを使うので別途インストールしてください*2

#coding:utf-8
from pylab import *
from collections import defaultdict

# statistics.py
# 訓練データとテストデータの各カテゴリの文書数をヒストグラムで描画

# カテゴリー名をロード
category = []
fp = open("train.map")
for line in fp:
    line = line.rstrip()
    category.append(line.split()[0])
fp.close()

# 訓練データの文書数をカテゴリーごとにカウント
train_count = defaultdict(int)
fp = open("train.label")
for line in fp:
    cat_idx = int(line)
    cat_idx = cat_idx - 1  # カテゴリIDは1から始まるので1引いておく
    train_count[cat_idx] += 1
fp.close()

# テストデータの文書数をカテゴリーごとにカウント
test_count = defaultdict(int)
fp = open("test.label")
for line in fp:
    cat_idx = int(line)
    cat_idx = cat_idx - 1
    test_count[cat_idx] += 1
fp.close()

# ヒストグラムで表示
pos = arange(len(category)) + 0.5
fig = figure(1)
# サブグラフ1
subplot(2, 1, 1)
bar(pos, train_count.values(), align="center")
ylim((0, 650))
ylabel("# data")
title("training data")
grid(True)
xticks(pos)
# サブグラフ2
subplot(2, 1, 2)
bar(pos, test_count.values(), align="center")
ylim((0, 650))
ylabel("# data")
title("test data")
grid(True)
xticks(pos, category)
fig.autofmt_xdate()
show()

データの整形

配布されているMatlab形式のtrain.dataとtest.dataは少し使いにくいので以下のような形式に変換します。

[category] [word:count] [word:count] ...  <- doc1
[category] [word:count] [word:count] ...  <- doc2
                ...
[category] [word:count] [word:count] ...  <- docN

各行が1つの文書のbag-of-words表現を表しています。各列はスペースで区切られていて1列目にはカテゴリID、2列目以降は単語:出現頻度の組がずらずら並びます。このファイル形式は、多くの機械学習ライブラリで採用されています。たとえば、SVMのライブラリlibsvm岡野原さんのオンライン機械学習ライブラリollなど。ベクトル空間モデルの場合は、単語:tf-idf値のようになります。本当はカテゴリや単語はそのまま書かず、カテゴリIDや単語IDに変換して数値化しますが、Pythonだと数値だろうとテキストだろうとどっちも大差なく扱えるのでわかりやすいように単語やカテゴリ名をIDに変換せずにそのまま使いました。下のような感じのデータを作成します。

alt.atheism archive:4 name:2 atheism:10 ...
comp.graphics name:1 last:1 modified:1 addresses:2 of:14 ...
talk.religion.misc of:4 and:4 other:1 are:2 the:14 in:5 ...

たとえば、1番目の文書は、alt.atheismというカテゴリでarchiveという単語が4回、nameが2回、atheismが10回出現することを意味しています。変換プログラムです。解凍した20 Newsgroupsの7つのファイルと同じフォルダで実行します。実行すると

  • news20 - 訓練データ集合
  • news20.t - テストデータ集合

という2つのファイルができます。

#coding:utf-8
import sys
from collections import defaultdict

# trans_data.py
# news20のデータセットをナイーブベイズで扱えるデータ形式に変換する
# http://people.csail.mit.edu/jrennie/20Newsgroups/
# 出力データフォーマット
# 1列目にカテゴリ、2列目以降は単語と出現頻度の組を列挙
# [category] [word:count] [word:count] ...  <- doc1
# [category] [word:count] [word:count] ...  <- doc2

def trans_data(labelfile, datafile, outfile):
    # カテゴリをロード
    category = []
    fp = open("train.map")  # test.mapでも同じ
    for line in fp:
        line = line.rstrip()
        category.append(line.split()[0])
    fp.close()
    
    # ボキャブラリをロード
    vocabulary = []
    fp = open("vocabulary.txt")
    for line in fp:
        line = line.rstrip()
        vocabulary.append(line)
    fp.close()
    
    # ラベルをロード
    train_label = []
    fp = open(labelfile)
    for line in fp:
        line = line.rstrip()
        # 配布ファイルのカテゴリIDは1から始まるので1をひく
        idx = int(line) - 1
        cat = category[idx]
        train_label.append(cat)
    fp.close()
    
    # 総文書数
    num = len(train_label)
    
    # 変換
    train_data = []
    for i in range(num):
        train_data.append([])
    
    fp = open(datafile)
    for line in fp:
        line = line.strip()
        temp = line.split()
        docidx, wordidx, count = int(temp[0]), int(temp[1]), int(temp[2])
        # 配布データの文書IDと単語IDは1から始まるので1を引いたインデックスを使う
        word = vocabulary[wordidx-1]
        train_data[docidx-1].append("%s:%d" % (word, count))
    fp.close()
    
    # ファイルに出力
    fp = open(outfile, "w")
    for i in range(num):
        fp.write("%s %s\n" % (train_label[i], " ".join(train_data[i])))
    fp.close()

if __name__ == "__main__":
    # 訓練データを変換
    trans_data("train.label", "train.data", "news20")
    # テストデータを変換
    trans_data("test.label", "test.data", "news20.t")

ナイーブベイズの精度評価

上の形式のファイルを読み込めるように前回作ったナイーブベイズを少しだけ変更します。変更箇所だけ再掲します。前回までは下のように同じ単語が複数あっても直に書いてましたが、今回からは単語:出現頻度というようにまとめている点が違うところです。ナイーブベイズのプログラム内で単語をカウントするときに「:」で分解して頻度分だけ足し込みます

Chinese Chinese Beijing
         ↓
Chinese:2 Beijing:1
class NaiveBayes:
    """Multinomial Naive Bayes"""
    def train(self, data):
        """ナイーブベイズ分類器の訓練"""
        # 文書集合からカテゴリを抽出して辞書を初期化
        ...
        # 文書集合からカテゴリと単語をカウント
        for d in data:
            cat, doc = d[0], d[1:]
            self.catcount[cat] += 1
            # ★変更箇所1
            for wc in doc:
                word, count = wc.split(":")
                count = int(count)
                self.vocabularies.add(word)
                self.wordcount[cat][word] += count
        # 単語の条件付き確率の分母の値をあらかじめ一括計算しておく(高速化のため)
        ...
    def score(self, doc, cat):
        """文書が与えられたときのカテゴリの事後確率の対数 log(P(cat|doc)) を求める"""
        total = sum(self.catcount.values())  # 総文書数
        score = math.log(float(self.catcount[cat]) / total)  # log P(cat)
        # ★変更箇所2
        for wc in doc:
            word, count = wc.split(":")
            count = int(count)
            # logをとるとかけ算は足し算になる
            for i in range(count):
                score += math.log(self.wordProb(word, cat))  # log P(word|cat)
        return score

def sample1():
    # Introduction to Information Retrieval 13.2の例題
    # データ表現を変更
    # 単語:出現回数
    # ★変更箇所3
    data = [["yes", "Chinese:2", "Beijing:1"],
            ["yes", "Chinese:2", "Shanghai:1"],
            ["yes", "Chinese:1", "Macao:1"],
            ["no", "Tokyo:1", "Japan:1", "Chinese:1"]]

では、訓練データnews20でナイーブベイズを訓練して、テストデータnews20.tで分類性能を評価してみます。分類性能にはaccuracyという指標を使います。これは、テストデータ集合全体に対してカテゴリを正しく予測できたデータの割合を表します。テストデータにはあらかじめ正解のラベルが付いているので、ナイーブベイズ分類器の予測カテゴリと一致するか調べればOK。以下、評価プログラムです。今回は英語文書なので不要ですが、あとで日本語文書にも使えるようにcodecsでファイルを開いてます。

#coding:utf-8
import codecs
from naivebayes import NaiveBayes

# eval.py
# ナイーブベイズの性能評価

def evaluate(trainfile, testfile):
    # 訓練データをロード
    trainData = []
    fp = codecs.open(trainfile, "r", "utf-8")
    for line in fp:
        line = line.rstrip()
        temp = line.split()
        trainData.append(temp)
    fp.close()
    
    # ナイーブベイズを訓練
    nb = NaiveBayes()
    nb.train(trainData)
    print nb
    
    # テストデータを評価
    hit = 0
    numTest = 0
    fp = codecs.open(testfile, "r", "utf-8")
    for line in fp:
        line = line.rstrip()
        temp = line.split()
        correct = temp[0]    # 正解カテゴリ
        words = temp[1:]     # 文書:単語の集合
        predict = nb.classify(words)  # ナイーブベイズでカテゴリを予測
        if correct == predict:
            hit += 1  # 予測と正解が一致したらヒット!
        numTest += 1
    print "accuracy:", float(hit) / float(numTest)
    fp.close()

if __name__ == "__main__":
    evaluate("news20", "news20.t")

実行結果は、

documents: 11269, vocabularies: 53975, categories: 20
accuracy: 0.785209860093

となりました。全テストデータのうち約78%のデータは正しく分類できてます。おー、かなりの高精度。この結果は、20 Newsgroupsで評価した論文

の結果ともほぼ一致します。

ストップワードの利用

変換したnews20とnews20.tのデータをよく見るとわかりますが、of, the, in, toのような意味のない単語が大量に含まれています。普通、これらの単語はストップワードを用いてあらかじめ除去します。英語のストップワードはいろいろなところで配布されていますが、Python自然言語処理ライブラリNLTKStopwords Corpusを使います。解凍して出てきたenglishが英語のストップワードリストで127単語が登録されています。これをstopwords.txtと名前を変えて保存します。

ストップワードを使うにはデータ変換プログラム(trans_data.py)を下のように修正します。

    # ボキャブラリをロード
    ...
    # ストップワードをロード
    stopwords = []
    fp = open("stopwords.txt")
    for line in fp:
        line = line.rstrip()
        stopwords.append(line)
    fp.close()
    # ラベルをロード
    ...

    fp = open(datafile)
    for line in fp:
        line = line.strip()
        temp = line.split()
        docidx, wordidx, count = int(temp[0]), int(temp[1]), int(temp[2])
        # 配布データの文書IDと単語IDは1から始まるので1を引いたインデックスを使う
        word = vocabulary[wordidx-1]
        if not word in stopwords:  # ★ストップワードに登録されていない単語だけ追加
            train_data[docidx-1].append("%s:%d" % (word, count))
    fp.close()

ストップワードを除去した訓練データ集合news20とテストデータ集合news20.tで精度を再評価してみます。結果は、

documents: 11269, vocabularies: 53852, categories: 20
accuracy: 0.802265156562

となります。さっきの結果と比べるとボキャブラリ数が53975から53852に減っていてストップワードが除去されているのがわかります。分類精度も78%から80%にわずかに向上しました。実はまだ分類を悪化させる無駄な単語が含まれていて精度を向上させる余地があります。そのような無駄な単語を除去する処理は特徴選択(feature selection)とか次元削減(dimensinality reduction)とか呼ばれていて、これまた多くの方法があります。有名なのは、相互情報量(Mutual Information)やカイ二乗値を用いる手法です。この2つはIIRにも載っているので後で実際に試してみたいと思います。どれくらい上がるんですかね?楽しみです。

まとめ

今回は、20 Newsgroupsというデータセットを用いてナイーブベイズの分類精度を評価してみました。単純なナイーブベイズでも80%というかなりの高精度をたたき出しました。私だったら80%なら十分実用的だねと満足してしまいそうですが、研究者の飽くなき執念はすさまじい。ナイーブベイズの結果が最低ラインでここからさらなる精度向上を目指して血道を上げることになります。

*1:Usenetって読んでたよ。fjとかなつかしいなぁ。スパムが多くてもう崩壊しかけてた。今もまだあるのだろうか。

*2:Mac OS Xmacportsでインストールするとグラフが表示されませんでした。その場合、~/.matplotlib/matplotlibrcにbackend: MacOSXを追加すればOK。

ナイーブベイズを用いたテキスト分類

今までPRMLを読んで実装を続けてきましたが、10章からは難しくて歯が立たなくなってきたのでここらで少し具体的な応用に目を向けてみようと思います。機械学習の応用先としては画像の方が結果を見ていて面白いんですが、当面は自然言語処理を取り上げます。そんなわけで一番始めの応用は機械学習と自然言語処理の接点として非常に重要なテキスト分類(Text Classification, Text Categorization)の技法たちを試していきたいと思います。テキスト分類は文書分類(Document Classification)という呼び方もあります。テキストと文書は同じ意味です。最初なので自分の知識の整理と入門者への紹介のためにちょっと丁寧にまとめてみました。

テキスト分類とは

テキスト分類とは、与えられた文書(Webページとか)をあらかじめ与えられたいくつかのカテゴリ(クラス)に自動分類するタスクです。テキスト分類は対象とするテキストによって幅広い応用が可能です。たとえば、すでに実用化されて身近でお世話になっている機能としては、

  • 電子メールを「スパム」と「それ以外」というカテゴリへ自動分類して「スパム」をゴミ箱へ捨てる(スパムフィルタ)
  • Webページを「政治・経済」「科学・学問」「コンピュータ・IT」「ゲーム・アニメ」などのカテゴリへ自動分類(はてなブックマーク)
  • ニュース記事を「興味あり」「興味なし」というカテゴリへ自動分類して「興味あり」のニュース記事だけおすすめ(情報推薦・情報フィルタリング)

などがあります。それぞれ、電子メール、Webページ、ニュース記事がテキストに当たります。たとえば、私も愛用しているはてなブックマークですが、人間がWebページの内容を読んで、このページは「コンピュータ・IT」だなとか分類しているわけではなく、機械学習の手法を用いた分類プログラム(分類器と呼ぶ)が自動的に分類しています。

大量のWebページが毎日毎日出てくるのにこんなの人手でできるはずないですよねー(Yahoo!は昔これを人手でやってましたが今はどうなんでしょうね?)。

教師あり学習

仕組みはこうです。まず、人間が教師となって分類器を訓練します。こんな感じ。

Webページ1は「IT」
Webページ2は「科学」
Webページ3は「IT」
Webページ4は「政治」
Webページ5は「ゲーム」
・・・

このような(テキスト,人間が与えた正解カテゴリ)を組としたデータを訓練データと呼びます。分類器はこの訓練データをもとに各カテゴリの文書の特徴を自動学習します。たとえば、

  • 「iPhone」「Apple」「Twitter」などの単語が含まれるテキストは「IT」カテゴリである確率が高い
  • 「民主党」「菅直人」などの単語が含まれるテキストは「政治」カテゴリである確率が高い
  • 「研究」「JAXA」「遺伝子」などの単語が含まれるテキストは「科学」カテゴリである確率が高い

などです。このように訓練した分類器を用いて、カテゴリがわからない新しい文書、たとえば、「Apple」「iPhone」が含まれる文書のカテゴリは?と分類器に聞くと「IT」である確率が高いと返してくれます。一般的に訓練データは多ければ多いほど分類器は正確なテキスト分類ができるようになります。このように、人間が正解カテゴリを訓練データとして与える機械学習手法は教師あり学習と呼びます。

Bag-of-words

一般的にテキストは単語の集合として与えます。集合なので並び順は無視されます。つまり、単語が文書内にどこに出てくるかは考慮しません。このようなテキスト表現はbag-of-wordsと呼ばれます。単語をバッグの中にぐちゃぐちゃ詰め込むイメージでしょうか。たとえば、

テキスト分類とは、与えられたテキストをあらかじめ与えられているカテゴリに「自動で」分類するタスクです。

という文書は、bag-of-wordsで表すと

テキスト テキスト カテゴリ タスク 自動 分類 分類 

みたいに単語の集合で表されます。タスクにもよりますが、形態素解析(2009/4/15)で名詞だけ抽出して使うことが多いんじゃないかと思います。話はそれますが、Visual Wordsを用いた類似画像検索(2010/2/27)で取り上げたbag-of-visual wordsはbag-of-wordsの画像版です。bag-of-visual wordsもbag-of-wordsと似ていて画像における単語(局所特徴量のセントロイド)が画像上のどこにあるかは考慮しません。このような単純化のおかげで学習アルゴリズムがシンプルになります。

テキスト分類の技法

テキスト分類は非常に多くの研究があり、そのアルゴリズムも大量にあります。ちょっと思いつくだけでも、ナイーブベイズ、決定木、Rocchio分類法、k-最近傍法、ロジスティック回帰、ニューラルネットワーク、サポートベクトルマシン、ブースティングなどなど。それぞれやり方はだいぶ違っています。また、テキストをベクトルへ変換する手法(TF-IDFとか)や次元削減の方法(LSIとか)もたくさん提案されており、その組み合わせを考えると結局どれ使えばいいの?って感じです。一般的には、サポートベクトルマシンやブースティングが他の手法と比べて高精度な分類ができると言われています。これから実際に試していきます。今回取り上げるのは、よく使われていて実装も簡単、しかも高速というナイーブベイズです。精度評価のベースラインとしてよく使われてます。

ナイーブベイズ

ナイーブベイズの中心となる式はベイズの定理を応用した下の式で表せます。


P(cat|doc) = \frac{P(cat) P(doc|cat)}{P(doc)} \propto P(cat) P(doc|cat)

事後確率P(cat|doc)は文書docが与えられたときカテゴリcatである確率です。カテゴリを予測したい未知の文書は、事後確率がもっとも高いカテゴリへ分類します(MAP推定)。この確率を計算するためには、右辺の事前確率P(cat)と尤度P(doc|cat)が必要になります。P(doc)はどのカテゴリにも共通なので無視できます。事後確率P(cat|doc)と尤度P(doc|cat)はややこしいのですが違うものです。私はこの違いを理解するのにだいぶ苦労した覚えがありますが・・・

まず、P(cat)ですがこれは簡単です。訓練データの各カテゴリの文書数の総文書数に占める割合を計算するだけです。たとえば、

訓練データ100文書中
IT      50文書  →  P(cat=IT)   = 50 / 100 = 0.5
科学    30文書  →  P(cat=科学)= 30 / 100 = 0.3
政治    20文書  →  P(cat=政治) = 20 / 100 = 0.2

のようになります。P(doc|cat)はちょっと複雑です。カテゴリcatが与えられたときに文書docが生成される確率です。ここで、文書docはbag-of-wordsで単語の集合 [word_1,word_2,...,word_k] として表され、単語間の独立性を仮定するとすると下のように計算できます。


P(doc|cat) = P(word_{1} \wedge \; \cdots \; \wedge word_{k} | cat) = \displaystyle \prod_i P(word_i | cat)

上の式で第2式から第3式へは単語の出現確率の間に独立性を仮定しないと成り立ちません。同時確率をそれぞれの確率の積で表せるってのが確率論的独立性の定義です。本来、単語の出現に独立性は成り立ちません。たとえば、「人工」と「知能」は共起しやすいし、「機械」と「学習」は共起しやすいです。これを無視して単語の出現は独立と無理矢理仮定して文書の確率を単語の確率の積で表して単純化するのがナイーブベイズのナイーブたる所以です。単語間の依存関係を仮定したナイーブベイズとしてTAN(Tree-Augmented Naive Bayes)というのも提案されていますが、あまり広まってないところを見ると労多くして功少なしって感じでしょうか?

で、今度はP(word_i|cat)の確率が必要です。これは、単語の条件付き確率と呼びます。カテゴリの中でその単語がどれくらいでてきやすいかを表します。これは簡単で、訓練データのカテゴリcatに単語word_kが出てきた回数をカテゴリcatの全単語数で割ればOKです。T(cat,word)をカテゴリcatに単語wordが出てきた回数、Vを訓練データ中の全単語集合(ボキャブラリ)とすると、


P(word_i | cat) = \frac{T(cat, word_i)}{\displaystyle \sum_{word' \in V} T(cat, word')}

となります。分母はVのすべての単語に関して足し合わせますが、実際は対象カテゴリcatに出てくる単語に絞っても結果は同じです。そのカテゴリに出てこなかった単語はT(cat,word)=0となるからです。

対数

以上の結果をまとめると最終的に分類されるカテゴリcat_mapは


cat_{map} = \arg \max_{cat} P(cat|doc) = \arg \max_{cat} P(cat) \displaystyle \prod_i P(word_i | cat)

となります。argmaxf(x)ってのはf(x)が最大になるようなxを返すっていう意味です。P(word|cat)というのは非常に小さい数な上に文書中にはたくさんの単語が含まれるのでかけ算部分がアンダーフローを起こす可能性があります。そこで、対数をとってかけ算を足し算化します。事後確率の大小関係は対数をとっても変化しない(結果となるcat_mapは変化しない)ので問題ありません。


cat_{map} = \arg \max_{cat} \log P(cat|doc) = \arg \max_{cat} \Bigl( \log P(cat) + \displaystyle \sum_i \log P(word_i | cat) \Bigr)

ゼロ頻度問題

P(doc|cat)は単語の条件付き確率P(word|cat)の積で求まったのですが、アンダーフロー以外にもう1つ大きな問題があります。それは、未知の文書のカテゴリを予測する際、訓練データのボキャブラリに含まれない単語を1つでも含んでいると単語の条件付き確率P(word|cat)は0となり、単語の条件付き確率の積で表されるP(doc|cat)も0となってしまうことです(対数のときはlog 0となり計算できなくなります)。つまり、その新しい文書が生成される確率は0になってしまいます。

たとえば、文書にiPhone、Appleなどの単語が含まれており、「おっ、これはカテゴリITから生成された可能性が高くなってきた」と思っていても、訓練時には含まれなかった新単語iPadが含まれてしまうとP(doc|cat) = 0となり、この文書がカテゴリITから生成された確率は0になってしまいます。iPhoneとAppleが出てるのだからカテゴリはITの可能性が高いだろ!これはおかしい!ってことになります。この問題は、ゼロ頻度問題と呼ばれています。ゼロ頻度問題は、スムージングという方法で緩和できます。よく使われるのが単語の出現回数に1を加えるラプラススムージング(Laplace Smoothing)です。新しい単語が出てくると確率は低くなりますが、0にはなりません。


P(word_i|cat) = \frac{T(cat, word_i) + 1}{\displaystyle \sum_{word' \in V} \Bigl (T(cat, word') + 1 \Bigr} = \frac{T(cat, word_i) + 1}{\displaystyle \sum_{word' \in V} T(cat, word') + |V|

Pythonで実装

上のを素直にPythonで実装すると下のようになります。対数をとり、ラプラススムージングを使っています。単語の条件付き確率P(word|cat)の分母は、引数の単語によらないため訓練時に事前に一括計算しています。これを単語の条件付き確率を求めるたびに計算しようとするとものすごく遅くなります。

#coding:utf-8
import math
import sys
from collections import defaultdict

class NaiveBayes:
    """Multinomial Naive Bayes"""
    def __init__(self):
        self.categories = set()     # カテゴリの集合
        self.vocabularies = set()   # ボキャブラリの集合
        self.wordcount = {}         # wordcount[cat][word] カテゴリでの単語の出現回数
        self.catcount = {}          # catcount[cat] カテゴリの出現回数
        self.denominator = {}       # denominator[cat] P(word|cat)の分母の値
    
    def train(self, data):
        """ナイーブベイズ分類器の訓練"""
        # 文書集合からカテゴリを抽出して辞書を初期化
        for d in data:
            cat = d[0]
            self.categories.add(cat)
        for cat in self.categories:
            self.wordcount[cat] = defaultdict(int)
            self.catcount[cat] = 0
        # 文書集合からカテゴリと単語をカウント
        for d in data:
            cat, doc = d[0], d[1:]
            self.catcount[cat] += 1
            for word in doc:
                self.vocabularies.add(word)
                self.wordcount[cat][word] += 1
        # 単語の条件付き確率の分母の値をあらかじめ一括計算しておく(高速化のため)
        for cat in self.categories:
            self.denominator[cat] = sum(self.wordcount[cat].values()) + len(self.vocabularies)
    
    def classify(self, doc):
        """事後確率の対数 log(P(cat|doc)) がもっとも大きなカテゴリを返す"""
        best = None
        max = -sys.maxint
        for cat in self.catcount.keys():
            p = self.score(doc, cat)
            if p > max:
                max = p
                best = cat
        return best
    
    def wordProb(self, word, cat):
        """単語の条件付き確率 P(word|cat) を求める"""
        # ラプラススムージングを適用
        # wordcount[cat]はdefaultdict(int)なのでカテゴリに存在しなかった単語はデフォルトの0を返す
        # 分母はtrain()の最後で一括計算済み
        return float(self.wordcount[cat][word] + 1) / float(self.denominator[cat])
    
    def score(self, doc, cat):
        """文書が与えられたときのカテゴリの事後確率の対数 log(P(cat|doc)) を求める"""
        total = sum(self.catcount.values())  # 総文書数
        score = math.log(float(self.catcount[cat]) / total)  # log P(cat)
        for word in doc:
            # logをとるとかけ算は足し算になる
            score += math.log(self.wordProb(word, cat))  # log P(word|cat)
        return score
    
    def __str__(self):
        total = sum(self.catcount.values())  # 総文書数
        return "documents: %d, vocabularies: %d, categories: %d" % (total, len(self.vocabularies), len(self.categories))

if __name__ == "__main__":
    # Introduction to Information Retrieval 13.2の例題
    data = [["yes", "Chinese", "Beijing", "Chinese"],
            ["yes", "Chinese", "Chinese", "Shanghai"],
            ["yes", "Chinese", "Macao"],
            ["no", "Tokyo", "Japan", "Chinese"]]
    
    # ナイーブベイズ分類器を訓練
    nb = NaiveBayes()
    nb.train(data)
    print nb
    print "P(Chinese|yes) = ", nb.wordProb("Chinese", "yes")
    print "P(Tokyo|yes) = ", nb.wordProb("Tokyo", "yes")
    print "P(Japan|yes) = ", nb.wordProb("Japan", "yes")
    print "P(Chinese|no) = ", nb.wordProb("Chinese", "no")
    print "P(Tokyo|no) = ", nb.wordProb("Tokyo", "no")
    print "P(Japan|no) = ", nb.wordProb("Japan", "no")
    
    # テストデータのカテゴリを予測
    test = ["Chinese", "Chinese", "Chinese", "Tokyo", "Japan"]
    print "log P(yes|test) =", nb.score(test, "yes")
    print "log P(no|test) =", nb.score(test, "no")
    print nb.classify(test)

上のプログラムでは、Introduction to Information Retrieval(IIR)のTable 13.1の例題を使っています。訓練データは、リストのリストで渡します。内側のリストが1つの訓練データです。リストの0番目の要素がカテゴリになります(よく使われる形式)。たとえば、1つめの訓練データは、bag-of-words表現で[Chinese, Beijing, Chinese]という文書がカテゴリyesであることを意味しています。4つの訓練データを与えてナイーブベイズ分類器を学習し、[Chinese, Chinese, Chinese, Tokyo, Japan]という文書のカテゴリを分類器で予測してます。IIRの結果と同じくyesに分類されます。以下、出力結果です。

documents: 4, vocabularies: 6, categories: 2
P(Chinese|yes) =  0.428571428571
P(Tokyo|yes) =  0.0714285714286
P(Japan|yes) =  0.0714285714286
P(Chinese|no) =  0.222222222222
P(Tokyo|no) =  0.222222222222
P(Japan|no) =  0.222222222222
log P(yes|test) = -8.10769031284
log P(no|test) = -8.906681345
yes

まとめ

ナイーブベイズには2つの代表的なモデルがあります。多項モデル(Multinomial Model)とベルヌーイモデル(Bernoulli Model)です。今回、実装したのは多項モデルです。私の印象では、多項モデルの方がよく使われている気がします。ベルヌーイモデルはあまり見かけません。2つの分類精度を比較した論文(McCallum,1998)によるとボキャブラリ数が多い場合は多項モデルの方が精度が高いことが示されています。ベルヌーイモデルは出現しない単語の確率も考慮するので計算量も大きいです。

今回はもっとも基礎的なテキスト分類のアルゴリズムであるナイーブベイズを実装してみました。用いた例題がすごく単純でありがたみがなかったので、次はスパムメールの分類やこのブログの記事カテゴリ(左にカテゴリーメニューってのがあります)を分類してみようと思います。

参考文献

補足

対数をとって大小を比較することで分類結果を出すことはできますが、分類結果を出すだけでなく、テストデータの各カテゴリへの事後確率 P(cat|doc) を求めたいときは下のようにします。


P(cat|doc) = \frac{P(cat) P(doc|cat)}{P(doc)} \propto P(cat) P(doc|cat)

の式でP(cat|doc)を計算すればいいわけですが、正規化係数(確率の和が1になるように調整するための係数)の分母のp(doc)を求めるのがけっこう大変です。そのため下のようなよく知られた裏技があります。

    def postProb(self, doc, cat):
        """文書が与えられたときのカテゴリの「正規化していない
       (=p(doc)で割らない)」事後確率 P'(cat|doc) を求める"""
        total = sum(self.catcount.values())  # 総文書数
        pp = float(self.catcount[cat]) / total  # 事前確率P(cat)
        # 尤度 P(doc|cat) = P(word1|cat) * p(word2|cat) * ...
        # 対数をとらないので掛け算になる(非常に小さな値!)
        for word in doc:
            pp *= self.wordProb(word, cat)
        return pp

    # テストデータの各カテゴリへの事後確率を求める
    test = ["Chinese", "Chinese", "Chinese", "Tokyo", "Japan"]
    p1 = nb.postProb(test, "yes")  # 正規化されていないので確率ではない!
    p2 = nb.postProb(test, "no")   # 正規化されていないので確率ではない!
    # 下のようにすると足して1になる確率になる
    print "P(yes|test) =", p1 / (p1 + p2)
    print "P(no|test)  =", p2 / (p1 + p2)

結果は、

P(yes|test) = 0.689758611763
P(no|test)  = 0.310241388237

となり、テストデータがyesである確率は69%、noである確率は31%となり、足すと1になる確率になってます。

もちろん、分母のp(doc)をp(cat1)p(doc|cat1) + p(cat2)p(doc|cat2) + ...のように展開して式どおりに計算しても同じ結果になります。