Python を使った実験データ処理

権藤研講習会資料 2019/09/14, 09/16, 10/07
新山

0. ウォームアップ問題
1. UNIXスクリプト処理の基本
2. 実験データの管理
3. 大量のデータを蓄積・処理する場合のTips
4. データのSerializationについて
5. Python から CSV/JSON/SQLite を使う
6. Vector Space Model を使った類似度計算
7. 類似度をもちいたクラスタリングの実装
7.99. ウォームアップ問題 (その2)
8. Python による決定木の実装
9. Python による Naive Bayes の実装

0. ウォームアップ問題

( ) 内はおおまかな目安時間です。

Python 3.7 をインストールせよ。(5分)
以下のようなCSVファイルから単語と出現回数を読み込み、回数が多い順にソートして表示するプログラムを書け。なお、同じ単語が複数回出てきた場合は、それらの合計を使用するものとする: (15分)
```
foo,3
baa,1
john,5
foo,1
```
出力例:
```
john,5
foo,4
baa,1
```
SQLite3 をインストールせよ。(5分)
以下の表を作成し値を入力する SQL を書け。(10分)

Id Name Score

1 Alcie 100

2 Bob 50

3 Carol 75
上の表から、Scoreが50点であるような行を表示する SQL を書け。(5分)
a001.txt, a002.txt, ... という名前のついた複数のファイルがあるとする。これらの名前を一括して a001.html, a002.html, ... に変更するようなシェルスクリプトを書け。 (5分)
上のスクリプトで、ファイル名の一覧があらかじめ決まっておらず、あるテキストファイル files.txt 中に一行ずつ書かれている場合はスクリプトをどう変更すべきか。(5分)

Id	Name	Score
1	Alcie	100
2	Bob	50
3	Carol	75

1. UNIXスクリプト処理の基本

1.1. シェルの基本

特殊な記号 ({, }, $, *, ;) を引数に渡す場合は '〜' で囲む。

演習. カレントディレクトリの中にある「-t」というファイル中から * という文字列を検索する grep コマンドを書け。

1.2. 変数の活用

a="foo"
echo "$a"
b="$a $a"
c='$a $a'
d='$a '"$a"

演習. シェル変数と環境変数の違いは何か?

1.3. Historyを活用する

「こんなコマンドを実行したはずなんだけど、なんだっけ」

$ history | grep なんか

注意: historyファイルはときどき消されることがある。

過去のコマンドラインをすべて記録しておく

function _prompt_cmd {
    local s=$?
    echo "`date '+%Y-%m-%d %H:%M:%S'` $HOSTNAME:$$ $PWD ($s) " \
         "`history 1`" >> $MYHISTFILE
}
PROMPT_COMMAND=_prompt_cmd

コマンドラインの記録は、そのまま実験ノートにもなる。あとで実験手順を再現したいときに参考になる。

1.4. パイプの使い方

演習. 以下のコマンドラインを順に実行し、つぎの表現が何をするか予測せよ。

$ ls -l
$ ls -l | wc
$ ls -l | sort
$ ls -l | sort -k4
$ ls -l | sort -k4 -r
$ ls -l | sort -k4 -r -n
$ ls -l | awk '{print $5;}'
$ ls -l | awk '{a+=$5;}'
$ ls -l | awk 'BEGIN{a=0;} {a+=$5;} END{print(a);}'
$ ls -l | awk '/euske/ {print $4;}' | uniq
$ ls -l | awk '/euske/ {print $4;}' | uniq -c
$ ls -l | awk '/euske/ {print $4;}' | uniq -c | wc

$ find ~
$ find ~ -type f
$ find ~ -type f | grep test
$ find ~ -type f | grep -i test
$ find ~ -type f | grep -i test | wc
$ find ~ -type f -name '*test*'
$ find ~ -type d -ctime -3
$ find ~ -type d -mtime +3
$ find ~ -type d -mtime +3

1.5. shスクリプトの基本

shスクリプトはこの行から始める。

#!/bin/sh
  または
#!/bin/bash

$ chmod 755 foo.sh

コマンドライン引数の扱い

echo "$0"
a=$1
shift
b=$1
c="$@"

あるプロセスの標準出力を値として使う

`コマンド`
  または
$(コマンド)

よく使う制御構造

if 式; then ...; else ...; fi

if 式; then
    ...
else
    ...
fi

for 変数 in 式; do
    ...
done

演習. 与えられた引数を1行ずつ表示するシェルスクリプトを書け。ただし、その行が foo であるときのみ、 bar と表示すること。

while read 変数; do
    ...
done

case 式 in
パターン1)
    ...
    ;;
パターン2)
    ...
    ;;
*)
    ...
    ;;
esac

1.6. xargs を使う

$ cat files.txt
a.txt
b.txt
c.txt
$ cat files.txt | xargs echo
$ cat files.txt | xargs cat
$ find -type d | xargs ls

1.7. 実験パイプラインの設計

UNIXプログラムのお約束:

出力は、標準出力に。
ログは、標準エラー出力に。

できるだけ各コマンドが ひとつだけのことをするように設計する。
できるだけ標準入力から入力を受けとり、標準出力に結果を表示する:
```
$ cmd1 input.txt | cmd2 | cmd3 > output.txt
```
途中結果を確認したい場合:
```
$ cmd1 input.txt
```

途中結果を保存しておき、後で続きをやりたい場合:

$ cmd1 input.txt > temp.txt
$ cmd2 < temp.txt | cmd3 > output.txt

途中結果をモニタしたい場合:

$ cmd1 input.txt | tee temp.txt | cmd2 | cmd3 > output.txt

a. フィルタとして設計する場合

$ コマンド [オプション] < 入力ファイル > 出力ファイル

b. 1つのファイルを入力する場合

$ コマンド [オプション] 入力ファイル

c. 可変個のファイルを入力する場合 (理想形)

$ コマンド [オプション] 入力ファイル1 入力ファイル2 ...

ファイル名が与えられない場合は標準入力を使用する。

こうしておくと

$ find ... | xargs コマンド [オプション]

のようにできる。

また、実験パラメータの変更はコードをじかに変更するのではなく、コマンドラインオプションとして処理すること。

例

$ exp1 -a1 -p2 -k0 input.txt > output_a1_p2_k0.txt
$ exp1 -a2 -p3 -k5 input.txt > output_a2_p3_k5.txt
...

ログに関する注意

ログは、たとえ人間が読む場合でも、なるべく機械的に処理できるようにしておくこと。 (grep/awk 等でのおおまかな統計が簡単にとれる。)

×: print("kokodayo", x, y)
○: print(f"DEBUG: data read complete: x={x}, y={y}.")

長く走らせるスクリプト

#!/bin/sh
exec </dev/null
exec >log
exec 2>&1
renice +20 -p $$

echo "*** START `date` ***"
/usr/bin/time -v 実際のコマンド
echo "*** END `date` ***"

1.8. Pythonスクリプトの定石

Pythonスクリプトは慣例によりこの行から始める。

#!/usr/bin/env python

だいたい以下のようなパターンで書くと、上に示した「お約束」に沿ったコマンドになる。

import sys
import fileinput

def doit(args):
    for line in fileinput.input(args):
        print(line)
    return

def main(argv):
    import getopt
    def usage():
        print('usage: %s [-d] [-o output] [file ...]' % argv[0])
        return 100
    try:
        (opts, args) = getopt.getopt(argv[1:], 'do:')
    except getopt.GetoptError:
        return usage()
    debug = 0
    output = None
    for (k, v) in opts:
        if k == '-d': debug += 1
        elif k == '-o': output = v
    return doit(args)

if __name__ == '__main__': sys.exit(main(sys.argv))

sys.argv でコマンドライン引数を取得。
getopt でオプションを取得。
main() 関数はコマンドラインの解析のみをおこない実際の処理は doit() 関数に正式なパラメータを渡してやらせる。 (これにより、他モジュールから利用することが可能)

$ python test.py input.txt
  または
$ test.py input.txt
  または
$ cat input.txt | test.py

演習.

上のプログラムを書き換え、変数 debug により doit() 内のなんらかの挙動が変わるようにせよ。
引数をとる -p オプションを追加せよ:
```
$ test.py -p 4 input.txt
```

2. 実験データの管理

2.1. ファイル名のつけ方

基本戦略は、シェルのワイルドカード (*) である条件をもったファイルだけを簡単に指定できるようにすることである。

名前は重要である。できるだけ意味がわかる名前にすること。 (さもないと、あとでワケわかんなくなる)
規則的にする。(日付・条件・用途によって予測可能にする)
- src_exp1_ver2.csv
- result-v4-20190803-temp.txt
- out.3cddd4d137ad4f794a8ccf0763b4d5a6450934b5
一定の prefix をつける。
- ×: jikken, jiken2, ..., kekka3
- ○: jikken_1, jikken_2, ..., jikken_3
データの種類・用途によって拡張子を変える。
- jikken_2.input
- jikken_2.output
- jikken_2.graph
- jikken_2.graph.svg
数値が入る場合は、桁数を揃える。
- ×: f1, f2, ..., f443
- ○: f0001, f0002, ..., f0443

2.2. ディレクトリ構造

基本的にUNIXのファイル名は逐次探索である。したがって、あまり1個のディレクトリに沢山のファイルを置くと遅くなる (せいぜい1000個程度)。

それからパス名が長くなりすぎると見にくいし、入力も大変。

データの種類・用途ごとにまとめる (input/, output/)
日付ごとにまとめる (s201909121012/, ...)
実験条件ごとにまとめる (data_seg01_p3_q4/, ...)
1., 2., 3. の混合 (data_201909121012_seg01_p3_input, ...)

3. 大量のデータを蓄積・処理する場合のTips

可能なかぎり stream可能にする (データ形式が重要)。

なるべく高速に parse できる形式にする。 (しかし自己流バイナリ形式はおすすめしない)

変更頻度が少ないものはディスク上に置いてもよい。

参照頻度が多くても、シーク可能ならディスク上に置けるかもしれない。 (OSがキャッシュするので)

高速化のためのよくある手段

計算をサボる。(事前条件の強化)
計算をごまかす。(近似値の使用)
計算を後回しにする/先にやっておく。(キャッシュ・lookup tableの使用)

4. データのSerializationについて

4.1. 考慮する要素

複雑さ (数値、テキスト、構造体)
データ量 (オンメモリ, 10GBytes, 1TBytes)
変更する頻度 (0回 / 1回以上)
参照パターン (シークの要・不要)
共有度 (複数のプロセス、タスク、マシン間で?)
キャッシュ可能性
堅牢性 (データの一部が崩れても全体に影響を与えないか?)

重要: できるだけ既存のツール・ライブラリで処理できるようにする。

4.2. テキストファイル (自分フォーマット)

おすすめしない。もしやるとしたら、parseが簡単にできるようにすること。

新山がときどき使っているフォーマット

# コメント
+キー1 バリュー1
+キー2 バリュー2
(空行がレコード区切り)

利点: 特定のフィールドが grep で簡単に見れる。データに注釈が入れられる。 Pythonで簡単に解析可能。エディタで修正が簡単。
欠点: 自己流。 1行 = 1レコードではない。

rec = {}
for line in fp:
    line = line.strip()
    if line.startswith('#'): continue
    if line.startswith('+'):
        (k,_,v) = line.partition(' ')
        rec[k] = v
    elif not line:
        yield rec
        rec = {}

4.3. バイナリファイル (自分フォーマット)

超おすすめしない。

簡単なデータだけならいいかも (たとえば 32ビット列の羅列ひたすら1G個とか)。

「SQLite は fopen() に対抗するために作られた」

4.4. よく知られている形式 (おすすめ)

CSV

利点: 簡単。テキスト。 Excelで編集できる。
欠点: 微妙に仕様が標準化されていない。 (とくに " や , が入ったセルの場合) 巨大なデータには向かない。シーク不可。書き換え不可。

JSON

利点: 簡単。テキスト。言語非依存。ストリーム可。
欠点: 巨大になるとエディタでは見にくい。シーク不可。書き換え不可。

XML

利点: 複雑な階層構造をもったデータ向け。テキスト。言語非依存。
欠点: 巨大になるとエディタでは見にくい。シーク不可。書き換え不可。

SQLite

利点: 型が決まっているデータ向け。言語非依存。 4Gぐらいのデータまでなら余裕。シーク可。書き換え可。堅牢。
欠点: あらかじめテーブルの設計が必要。ツールを使わないと見れない。

SQLite + JSON

複雑な構造 × 膨大な数があるときに使う方法。

4.5. もっと高度な方法

ProtocolBuffer, HDF, MongoDB, ...

導入に手間がかかりすぎて、個人でやる実験には向かない。

5. Python から CSV/JSON/SQLite を使う

5.1. CSV

書き込み

import csv
with open('output.csv', 'w') as fp:
    writer = csv.writer(fp)
    writer.writerow(['a', 'b', 'ccc'])

読み込み

import csv
with open('input.csv') as fp:
    for row in csv.reader(fp):
        print(row)

5.2. JSON

書き込み

import json
with open(output.json') as fp:
    data = {'a':123, b:['x','y']}
    fp.write(json.dumps(data))

読み込み

import json
with open('input.json') as fp:
    for line in fp:
        data = json.loads(line)

5.3. SQLite

書き込み

import sqlite3
db = sqlite3.connect('data.db')
cur = db.cursor()
cur.executescript('''
CREATE TABLE Student (
    Id INTEGER PRIMARY KEY,
    Name TEXT,
    Score INTEGER);
''')
for (name,id,score) in scores:
    cur.execute('INSERT INTO Student VALUES (?, ?, ?);', name, id, score)

読み込み

import sqlite3
db = sqlite3.connect('data.db')
cur = db.cursor()
for row in cur.execute('SELECT Name,Id FROM Student;'):
    (name,id) = row

6. Vector Space Model を使った類似度計算

Vector Space Model (VSM) とは離散的な素性からなるデータの類似度を計算する簡単な方法である。

A = {a_i}, B = {b_i}

というデータがあるとき、簡単な方法では:

Jaccard係数 = (A ∩ B) / (A ∪ B)

VSM はこれをもうちょっと高度にしたもので、各素性の「重み (頻度)」を考慮する。これはようするに {a_i}, {b_i} の各要素を並べ、これをベクトルとしてみたときの cosine距離である:

類似度 = Σ (a_i × b_i) / √ (Σ |a_i|²) × √ (Σ |b_i|²)

これがなぜうまくいくのかは、次の例を考えてみるとわかる。

素性が揃っているとき	素性が揃ってないとき

類似度: 高い	類似度: 低い

さらに、次のようなケースもカバーできる:

高さが違うが形は似ているとき

類似度: 高い

6.1. Pythonによる実装

各素性の集合からなるベクトルは、Python では辞書 (dict) として表すのがもっとも自然である。そこで、2つの与えられた辞書オブジェクトから cosine距離を求めるような関数 calcdot() を作る。

def calcdot(a, b):
    ...
    return

v1 = {'x':1, 'y':2}
v2 = {'x':2, 'y':1, 'z':3}
print(calcdot(v1, v2))  # 0.47809144373375745

演習. 上の関数を完成させよ。

6.2. 実際のデータへの適用

Wikipedia英語版の記事データ (10000記事)

演習a. 上の記事データをダウンロードし、各記事中の単語を単純な正規表現で切り出してカウントする Python プログラムを書け。このデータは以下のような構造になっている。

# 2428190 Melbourne Shuffle
The Melbourne Shuffle (also known as Rocking or simply The Shuffle) is
a rave and club dance that originated in the late 1980s in the
underground rave music scene in Melbourne, Australia. The basic
...
(空行)
# 442370 List of prime numbers
By Euclid's theorem, there is an infinite number of prime
numbers. Subsets of the prime numbers may be generated with various
formulas for primes. The first 500 primes are listed below, followed
...
(空行)

ヒント

まず文字列を単語のリストに変換する関数 splitwords を考える。とりあえず、正規表現を使って簡単にやる:

import re
def splitwords(text):
    return [ w.lower() for w in re.findall(r'\w+', text) ]

つぎに上の doit() を改良して、読み込みんだデータファイルを解析する:

def doit(args):
    for line in fileinput.input(args):
        line = line.strip()
        if line.startswith('#'):
            (artid, _, title) = line[2:].partition(' ')
            artid = int(artid)
        elif line:
            # 単語に区切る。
            words = splitwords(line)
        else:
            # この時点で artid, title, words が設定されているはず。
            print(artid, title, words)
            # 各単語の頻度情報を格納した辞書 wordcount を求める。
            wordcount = countwords(words)

gzip圧縮されたデータを読み込むには、以下のようにする手もあるが:

$ gzip -dc enwiki-20140102-10000.txt.gz | python doit.py

Pythonのfileinputにオプションを与えると、gzipをそのまま読み込むことができる。 (ただし、UTF-8をデコードする必要があるので注意!)

    for line in fileinput.input(args, openhook=fileinput.hook_compressed):
        ...

演習b. 各記事ID, 記事タイトルおよび a. でカウントした頻度情報 (JSON形式) を入れる SQLテーブルを設計し、SQLiteを使って10000記事ぶんの情報を格納したデータベースを作成せよ。

ヒント

まずSQLiteのテーブルを作成しよう。

db = sqlite3.connect('articles.db')
cur = db.cursor()
cur.executescript('''
CREATE TABLE Article (
    ArtId INTEGER PRIMARY KEY,
    Title TEXT,
    Words TEXT);
''')

ここに、上で求めた (artid, title, wordcount) を INSERT していく。ただし、wordcount は辞書なので、JSONに変換する:

cur.execute('INSERT INTO Article VALUES(?,?,?);',
            (artid, title, json.dumps(wordcount)))
# 注意: commitしないとデータは書き込まれない。
db.commit()

演習c. 上で設計した関数 calcdot() を 10000×10000のベクトル対に適用し、もっとも高い類似度をもつ記事ペアを表示せよ。

ヒント

これは別の Python スクリプトにする。

基本戦略は、すべての記事対 |A| × |B| に対して、 calcdot(a,b) を計算し、最高となる a, b を求めればよい:

articles = [ ... ]
maxsim = 0
maxpair = None
for a in articles:
    for b in articles:
        sim = calcdot(a, b)
        if maxsim < sim:
            maxsim = sim
            maxpair = (a,b)

実際には calcdot(a,b) == calcdot(b,a) であることを利用して、計算時間を半分にする:

for (i,a) in enumerate(articles):
    for b in articles[i+1:]:
        sim = calcdot(a, b)
        ...

問題は、articles は巨大で一度にメモリに読み込みたくないということである。そこで、まず artid の一覧を取得し毎回 SQLite から記事のベクトルを読み込むことにする。

def getarticle(cur, artid):
    for (wordcount,) in cur.execute(
        'SELECT Words FROM Article WHERE ArtId = ?;', (artid,)):
        return json.loads(wordcount)

# すべての Artid のリストを取得。(これは小さい)
artids = [ artid for (artid,) in cur.execute('SELECT ArtId FROM Article;') ]

(注意: なお、ここで紹介した方法は完璧ではない。一般的には、自然言語文の類似度計算には各単語の出現頻度だけでなく、単語の重み (IDF) も考慮している)

新山による実装: https://github.com/euske/python3-toys/blob/master/vsm.py

7. 類似度をもちいたクラスタリングの実装

上で求めた類似度の計算を使って、類似記事をグループ化 (クラスタリング) してみる。ここではもっとも簡単に、クラスタ間の最短距離 (半径) でまとめることを考える:

# しきい値
threshold = 0.9

for ...:
    a = getarticle(artid1)
    b = getarticle(artid2)
    sim = calcdot(a, b)
    if threshold < sim:
       # まとめる
       connect(artid1, artid2)

関数 connect(p, q) は、すこし複雑である:

# 各グループの所属を保持する辞書:
groups = {}

def connect(p, q):
    if p in groups and q in groups:
        # p, q 両方が別々のグループに含まれている場合: g[p] と g[q] を統合する。
        g = groups[p]
        g.merge(groups[q])
        # 各要素の所属をつけかえる。
        for x in g.artids:
            groups[x] = g
    elif p in groups:
        # p のみがグループに含まれている場合: q を g[p] に追加する。
        groups[p].add(q)
    elif q in groups:
        # q のみがグループに含まれている場合: p を g[q] に追加する。
        groups[q].add(p)
    else:
        # p, q どちらもグループに含まれていない場合: 新規グループを作成する。
        g = Group()
        g.add(p)
        g.add(q)
        groups[p] = groups[q] = g

Group はクラスを使ってこのように定義する。

class Group:

    def __init__(self):
        self.artids = []
        return

    def add(self, artid):
        self.artids.append(artid)
        return

    def merge(self, group):
        self.artids.extend(group.artids)
        return

7.99. ウォームアップ問題 (その2)

CSV ファイルを以下のような JSON 形式に変換するスクリプト csv2json.py を書け。

8. Python による決定木の実装

決定木は、少数の素性 (〜1000個程度) のみからなるデータを学習したいときに有用である。

利点

誰でも実装できる。
一度学習してしまえばメチャクチャ高速。

欠点

過学習しやすい。
確からしさを出したいときには使えない。

基本的な手順

使う素性をいっこ選び、その素性を使った「split」を列挙する。
もっとも相対エントロピーが高い split を選択する。
再帰的に split を繰り返す。

8.1. 準備

まず、学習する対象のモノ (オブジェクト) を Python でどのように表すかを決定する。オブジェクトは基本的に素性の集合なので、ここでは単に Python の辞書を使う。

obj = {'Temperature': 'High', 'Wind': 'Strong' }

つぎに Feature (素性) クラスを設計する:

class Feature:

    def __init__(self, attr):
        self.attr = attr
        return

    # 与えられたオブジェクトから素性を取り出す。
    def ident(self, obj):
        return obj[self.attr]

    # 与えられたオブジェクト列をこの素性で分割する。
    def split(self, objs):
        raise NotImplementedError  # 未実装

こうすると何が嬉しいかというと、「素性」というものを抽象的に定義することで、あとあと別のタイプの素性 (連続的な量を表す素性など) に拡張できるためである。

>>> f = Feature('Temperature')  # 素性 Temperature を定義する。
>>> f.ident(obj)                # obj からその素性を取り出す。
'High'

つぎに、作成する決定木をあらわす抽象クラス Node を設計する:

class Node:

    # 与えられたオブジェクトに対する判定結果を返す。
    def test(self, obj):
        raise NotImplementedError  # 未実装

実際には、決定木は再帰的な構造をもっていて、これは Branch と Leaf に分けられる。それぞれを実装する:

class Leaf(Node):

    def __init__(self, answer):
        self.answer = answer
        return

    # 与えられたオブジェクトに対する判定結果を返す。
    def test(self, obj):
        # これは木の終端なので、すでに答えは決まっている。
        return self.answer

class Branch(Node):

    def __init__(self, feature, children):
        self.feature = feature
        self.children = children
        return

    # 与えられたオブジェクトに対する判定結果を返す。
    def test(self, obj):
        # 決められた素性を使って、子ノードのいずれかを選ぶ。
        v = self.feature.ident(obj)
        branch = self.children[v]
        # 子ノードにオブジェクトを渡して予測させる。
        return branch.test(obj)

ひとつのオブジェクト obj はルート Branch から始まって、次々に test(obj) に渡されていく。木を下るにしたがって、いずれ Leaf ノードに到達する。 Leaf ノードはつねに決まった値を返すので、この時点で判定結果が決まることになる。

ここで先ほどの抽象的な Feature クラスを継承し、離散的な素性を扱うための DiscreteFeature クラスを実装する:

class DiscreteFeature(Feature):

    # 与えられたオブジェクト列をこの素性で分割する。
    def split(self, objs):
        assert 2 <= len(objs)
        # 自分の素性に従って objs を分割する。
        d = {}
        for obj in objs:
            v = self.ident(obj)
            if v in d:
                a = d[v]
            else:
                a = d[v] = []
            a.append(obj)
        # この時点で
        # d = { 値1: [オブジェクト, オブジェクト, ...],
        #       値2: [オブジェクト, オブジェクト, ...],
        #       ... }
        # 各項目の平均エントロピーを計算する。
        n = len(objs)
        ent = sum( len(a) * calcent(a) for a in d.values() ) / n
        # 平均エントロピーと分割結果を返す。
        return (ent, d)

8.2. エントロピーの計算

ある分布をもつ確率変数 P のエントロピーは次の式で表される:

E(P) = Σ -p × log(p)

ここで、与えられたオブジェクト列に対するエントロピーを計算する関数 calcent(objs) を作成する。ここでいう「オブジェクト列に対するエントロピー」とは、実際には判定結果のエントロピーのことであるので、まず学習データの各オブジェクトのどこに判定結果が入っているかをあらかじめ決めておく。ここでは obj['answer'] に入っているとしよう。 (なお log の底は自然対数でもよいが、2 を使うほうがわかりやすい)

# 与えられた学習データ objs のエントロピーを計算する。
def calcent(objs):
    answers = [ obj['answer'] for obj in objs ]
    probs = []
    ...
    return sum( -p*log2(p) for p in probs )

演習. 上の関数 calcent() を完成させよ。

8.3. 決定木の構築

以上の準備ができたら、実際に与えられた素性 feats と学習データ objs から決定木を構築する関数 buildtree を作成する。

# buildtree:
#   素性の集合 feats を使って、オブジェクト列 objs を分類する決定木を返す。
def buildtree(feats, objs):
    # これ以上分岐しても意味がない場合は Leaf を返す。
    if calcent(objs) < EPSILON:
        return Leaf(getbest(objs))
    # objs をもっともよく分割するような素性を探す。
    minent = 9999
    minfeat = None
    minsplit = None
    for feat in feats:
        (ent, split) = feat.split(objs)
        if ent < minent:
            # もっともエントロピーが少ない時の feat と split を記録。
            minent = ent
            minfeat = feat
            minsplit = split
    assert minsplit is not None
    # split の各結果に対してさらに部分木を構築する。
    children = {}
    for (v,a) in minsplit.items():
        children[v] = buildtree(feats, a)
    return Branch(minfeat, children)

変数 feats にはあらかじめ定義した Feature の列を与える。 picnic.csv の場合、これは "Outlook", "Temp", "Humidity", "Wind" の 4種類を与えればよい。

Day,Outlook,Temp,Humidity,Wind,Decision
1,Sunny,Hot,High,Weak,No
2,Sunny,Hot,High,Strong,No
3,Overcast,Hot,High,Weak,Yes
...

feats = [
    DiscreteFeature('Outlook'),
    DiscreteFeature('Temp'),
    DiscreteFeature('Humidity'),
    DiscreteFeature('Wind')
]
tree = buildtree(feats, objs)

8.4. 決定木を使う

完成した決定木を使うには、トップの Node.test() メソッドに判定したいオブジェクト (素性の集合) を渡せばよい。子ノードの test() が再帰的に呼ばれ、解が決定される。

obj = {'Outlook': 'Sunny', 'Temp': 'Hot', 'Humidity': 'High', 'Wind': 'Weak'}
print(tree.test(obj))

演習. picnic.csv のデータを使って決定木を学習し、実際に適用せよ。

実際には、学習した Leaf / Branch オブジェクトをファイルに保存しておき、あとで使うといった工夫が必要になる。

新山による実装: https://github.com/euske/python3-toys/blob/master/dtree.py

9. Python による Naive Bayes の実装

Naive Bayes 法は、離散的な少数の素性 (〜1000個程度) のみからなるデータを学習したいときに有用である。

利点

学習が簡単 (数えるだけ)。
確率の異なる複数の結果を返せる。

欠点

複雑な相関関係をもつデータには不向き。
一度も学習データに表れなかった素性があるとつねに「確率ゼロ」になってしまう。

9.1. 原理

学習: 素性の集合 F = {f_i} とクラス k に対して、条件つき確率 P(k | f₁, f₂, ..., f_n) を計算する。
予測: F が与えられたら、 argmax_k P(k | F) となるような k を求める。

で、どうやって P(k | F) を計算するのか? 以下の式を使う。 (Naive Bayes といわれるゆえんである。)

P(k | F) = P(k)・P(F|k) / P(F)

ここで F はあらかじめ与えられているので P(F) は無視できて

argmax_k P(k | F) = argmax_k P(k)・P(F | k)

さらに、各素性 f_i と k は相関しているが、各素性 f₁, f₂, ... f_n どうしはそれぞれ独立 (independent) して現れると仮定する。つまり

P(f₁ | f₂) = P(f₁ | f₃) = ... P(f₁ | f_n) = P(f₁)

実際には、この仮定は正しくない。 (Naive Bayes といわれるゆえんである。) しかしこの仮定により、

P(F | k) = P(f₁, f₂, ..., f_n | k) = P(f₁ | k) × P(f₂ | k) × ... × P(f_n | k)

と近似することができる。 P(f_i | k) を求めるのは簡単である。学習データを見て、各素性 f_i とそれに対応するクラス k が同時に現れる回数をただ数えればよい。このように、 Naive Bayes では素性と応答の数をただかぞえるだけで学習が可能である。

9.2. Python における実装

Naive Bayes を実装するには

各クラス k が学習データ中に何回現れたか。
各クラスと素性の対 (f_i, k) が学習データ中に何回現れたか。

を記録しておく必要がある。 a. は簡単である。いっぽう b. は、以下のように格納しておくと便利である:

fcount = {
  クラス1: { 素性a: 回数, 素性b: 回数, ... }
  クラス2: { 素性a: 回数, 素性b: 回数, ... }
  ...
}

さらに (素性と関係なく) k が現れた回数は、それ自体をひとつの特殊な素性 ALL とみなせるので

fcount = {
  クラス1: { ALL: 回数, 素性a: 回数, 素性b: 回数, ... }
  クラス2: { ALL: 回数, 素性a: 回数, 素性b: 回数, ... }
  ...
}

のようにできる。

これをふまえて、 NaiveBayes クラスを定義する:

class NaiveBayes:

    def __init__(self):
        self.fcount = {}  # 素性 (k,f) の出現回数。
        return

    # 素性とクラスの相関をひとつ学習する。
    def learn(self, k, feats):
        # クラス k と同時に現れた素性一覧をとりだす。
        if k in self.fcount:
            fc = self.fcount[k]
        else:
            fc = self.fcount[k] = {}
        # k の数を数える。
        if 'ALL' not in fc:
            fc['ALL'] = 0
        fc['ALL'] += 1
        # (f,k) の数を数える。
        for f in feats:
            if f not in fc:
                fc[f] = 0
            fc[f] += 1
        return

9.3. 予測する

モデルが学習できたら、予測である。素性の集合 feats が与えられたら、各 k に対して

P(k)・Π P(f_i | k) = P(k)・Π {P(f_i, k) / P(k)}

を計算すればよいのであるが、実際には母数 N が同じなのでこれは確率である必要がない。

fcount[k]['ALL'] * Π (fcount[k][f] / fconut[k]['ALL'])

だけで済んでしまう。さらに、あらかじめ kcount と fcount の log を記録しておき

log(fcount[k]['ALL']) + Σ (log(fcount[k][f]) - log(fconut[k]['ALL']))

のようにすれば加減算だけでよくなる。これをふまえて、メソッド predict() を設計する:

class NaiveBayes:
    ...

    # 与えられた素性から推定される各クラスの確率を返す。
    def predict(self, feats):
        klass = []
        for (k,fc) in self.fcount.items():
            # P(k)・P(f_i | k) を計算する。
            pk = log(fc['ALL'])
            p = (pk + sum( (log(fc[f]) - pk) for f in feats ))
            klass.append((p, k))
        # クラスの一覧を確率の大きい順にソートする。
        klass.sort(reverse=True)
        return klass

この方法がよいのは、結果が確率 (のlog) つきで返されるということである。もっとも確実な予想だけを知りたければ klass[0] を使えばよいし、第2候補も欲しければ klass[1] も見ればよい。複数の候補が返されるのは Naive Bayes の大きな利点である。

注意点

決定木における「素性」とは、何がしかの値を持つものであったが、 Naive Bayes における「素性」は、実際には「存在するか否か」という二値的なものであることに注意。したがって、 picnic.csv のようなデータを使うには、各素性を "Outlook=Sunny" のようにまるごと文字列として表してやる必要がある。つまり「素性 Outlook の値が "Sunny" / "Overcast" / "Rain" のどれかだ」と考えるのではなく、「"Outlook=Sunny"、 "Outlook=Overcast"、 "Outlook=Rain" という別々の素性が存在する」と考えるのである。当然、Outlook の値が排他的だという情報は Naive Bayes にはわからない。したがって Naive Bayes は「"Outlook=Sunny" かつ "Outlook=Overcast"」というありえない状況も排除しない。これは独立性の仮定を置いたことによる帰結で、 Naive Bayes 法の限界である。

演習. 同じく、上で示した picnic.csv のデータを今度は NaiveBayes クラスに適用し、結果を観察せよ。

nb = NaiveBayes()
FEATS = ( 'Outlook', 'Temp', 'Humidity', 'Wind' )
for obj in objs:
    # オブジェクトの各素性の値を、二値的な素性に変換する。
    feats = [ f'{k}={obj[k]}' for k in FEATS ]
    # Decision の値と各素性との相関を学習する。
    nb.learn(obj['Decision'], feats)

新山による実装: https://github.com/euske/python3-toys/blob/master/naivebayes.py

9.4. スムージング

実際に上の例を実行してみると、 3番目のオブジェクトのあたりで「"Outlook=Overcast" という素性が存在しない」という KeyError例外が発生してしまう。これは P(No | Outlook=Overcast) の確率を計算しようとしたことによる。 (fcount['No'] の fc 中には Outlook=Overcast というキーが存在しない。) これは Naive Bayes を使うさいによく現れる問題で、このような学習データが存在しなかったのであるから、そもそも確率を計算できないのである。

このような場合、逃げの一種として「素性の各出現回数に 1 を出す」という方法がある。いわゆる "Laplace smoothing" である。これは Python のコード上では、 fc[f] でキーが存在しなかったときに 1 を返すように実装するだけである。

# 使用前
p = (pk + sum( (log(fc[f]) - pk) for f in feats ))
# 使用後
p = (pk + sum( (log(fc.get(f,0)+1) - pk) for f in feats ))

Yusuke Shinyama