2017年06月19日

 scikit-learnのナイーブベイズ分類器を使ってみる

前回のFamily Out Problemの確率モデル(下図)を題材に、ナイーブベイズ分類器を使ってみる。


Family Out Problem

この問題において、family-out以外の変数の真偽値が与えられた時にfamily-out=TRUEかどうかを判定するよう、ナイーブベイズ分類器を学習させてみる。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.naive_bayes import BernoulliNB
from sklearn.metrics import precision_recall_fscore_support
from sklearn.metrics import roc_curve, auc

def generate_sample():
    "Sample data generator of the Family Out problem"
    fo = np.random.binomial(1, 0.15)
    bp = np.random.binomial(1, 0.01)
    lo = np.random.binomial(1, (0.05, 0.6)[fo])
    do = np.random.binomial(1, ((0.3, 0.97), (0.9, 0.99))[fo][bp])
    hb = np.random.binomial(1, (0.01, 0.7)[do])
    return [fo, bp, lo, do, hb]

# Generate training data and test data
train_data = np.array([generate_sample() for _ in range(1000)])
test_data = np.array([generate_sample() for _ in range(100)])

X_train = train_data[:, 1:5]	# other than fo
y_train = train_data[:, 0]	# only fo
X_test = test_data[:, 1:5]	# other than fo
y_test = test_data[:, 0]	# only fo

# Train a Naive Bayes classifier
clf = BernoulliNB()
clf.fit(X_train, y_train)

# Evaluate with training data
y_pred = clf.predict(X_train)
metrics = precision_recall_fscore_support(y_train, y_pred)
print('Evaluation with training data')
print('Class FamilyOut=False: Precision={:.3f}, Recall={:.3f}, F-measure={:.3f}'.format(metrics[0][0], metrics[1][0], metrics[2][0]))
print('Class FamilyOut=True : Precision={:.3f}, Recall={:.3f}, F-measure={:.3f}'.format(metrics[0][1], metrics[1][1], metrics[2][1]))

# Evaluate with test data
y_pred = clf.predict(X_test)
metrics = precision_recall_fscore_support(y_test, y_pred)
print('Evaluation with test data')
print('Class FamilyOut=False: Precision={:.3f}, Recall={:.3f}, F-measure={:.3f}'.format(metrics[0][0], metrics[1][0], metrics[2][0]))
print('Class FamilyOut=True : Precision={:.3f}, Recall={:.3f}, F-measure={:.3f}'.format(metrics[0][1], metrics[1][1], metrics[2][1]))

# Draw ROC curve
y_train_post = clf.predict_proba(X_train)[:, 0]
y_test_post = clf.predict_proba(X_test)[:, 0]

fpr, tpr, thresholds = roc_curve(y_train, y_train_post, pos_label=0)
roc_auc = auc(fpr, tpr)
plt.plot(fpr, tpr, 'k-', lw=2, label='ROC for training data (area = {:.2f})'.format(roc_auc))

fpr, tpr, thresholds = roc_curve(y_test, y_test_post, pos_label=0)
roc_auc = auc(fpr, tpr)
plt.plot(fpr, tpr, 'k--', lw=2, label='ROC for test data (area = {:.2f})'.format(roc_auc))

plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('ROC curve of BernoulliNB')
plt.legend(loc="lower right")

plt.show()

問題のモデルに従って学習用データとテストデータをランダムに生成し、sklearn.naive_bayes.BernoulliNBの分類器を学習させ、テストデータを分類させ、family-out=FALSEとfamily-out=TRUEのそれぞれについてPrecision, Recall, F値を計算している。
加えて、予測したfamily-out=FALSE/TRUEの確率から、ROC曲線とAUCを出力している。

出力例

Evaluation with training data
Class FamilyOut=False: Precision=0.936, Recall=0.983, F-measure=0.959
Class FamilyOut=True : Precision=0.847, Recall=0.589, F-measure=0.695
Evaluation with test data
Class FamilyOut=False: Precision=0.923, Recall=0.988, F-measure=0.955
Class FamilyOut=True : Precision=0.889, Recall=0.533, F-measure=0.667
ROC曲線
ROC curve
その時のtraining_dataとtest_data(.arff形式)
family_out_train.arff
family_out_test.arff

ナイーブベイズ分類器には、同じクラスのデータでは全ての特徴が独立に出現する、つまりbp,lo,do,hbに依存関係が無く、これらがTRUEになる確率はfoの値だけで決まるという仮定があるが、この問題ではbp,do,hbに依存関係があるので、ベイジアンネットワークの方がうまく学習できると考えられる。
なのでベイジアンネットワークで分類した場合の性能と比較したいが、scikit-learn(0.18)にベイジアンネットワークが無いので、代わりにWeka 3.8.1のBayesNetを用いて学習用データで学習させ、テストデータを分類した結果と比較してみる。

Wekaの操作手順

  1. Weka ExplorerのPreprocessタブでOpen file...ボタンを押し、family_out_train.arffを開く
  2. ClassifyタブでClassifierとしてBayesNetを選択する
  3. "BayesNet -D ..."とあるフィールドをクリックしてパラメーター設定画面を開き、searchAlgorithmのフィールドをクリックして、initAsNaiveBayes=False, maxNrOfParents=2に変更
  4. Test optionsのSupplied test setのボタンを押し、family_out_test.arffを開き、Classとしてfamily-outを選択
  5. クラスをfamily-outを選択し、Startボタンを押す
上記と同じデータを使った出力例
=== Run information ===

Scheme:       weka.classifiers.bayes.BayesNet -D -Q weka.classifiers.bayes.net.search.local.K2 -- -P 2 -N -S BAYES -E weka.classifiers.bayes.net.estimate.SimpleEstimator -- -A 0.5

=== Detailed Accuracy By Class ===

                 TP Rate  FP Rate  Precision  Recall   F-Measure  MCC      ROC Area  PRC Area  Class
                 0.988    0.467    0.923      0.988    0.955      0.651    0.889     0.967     0
                 0.533    0.012    0.889      0.533    0.667      0.651    0.889     0.648     1
Weighted Avg.    0.920    0.398    0.918      0.920    0.911      0.651    0.889     0.919     

=== Confusion Matrix ===

  a  b   <-- classified as
 84  1 |  a = 0
  7  8 |  b = 1

BernoulliNBのテストデータに対する出力と、Precision, Recall, F-measureが完全に一致している。

WekaのBayesNetのデフォルト設定では正しいネットワークが学習されなかったが、上記手順の3.のように設定を変えると、大体次のように正しいネットワークになった。

BernoulliNBの性能には試行毎にばらつきがあったが、1000回の平均を取ってみたのが次の値であり、上記の結果は特別に良い例ではない。

Average with training data
Class FamilyOut=False: Precision=0.925, Recall=0.982, F-measure=0.952
Class FamilyOut=True : Precision=0.844, Recall=0.545, F-measure=0.661
Average with test data
Class FamilyOut=False: Precision=0.925, Recall=0.982, F-measure=0.952
Class FamilyOut=True : Precision=0.842, Recall=0.544, F-measure=0.651
Average of AUC with training data=0.894
Average of AUC with test data=0.895
また、計10回、同じデータでWekaのBayesNetの性能と比較した所、Precision, Recall, F-measureについては、8回は一致していた。
ROC曲線のAUCはWekaのBayesNetの方が高いことが多かったが、BernoulliNBでも大体0.85-0.90の範囲であり、十分に高かった。
従って、この確率モデルに対して、ナイーブベイズ分類器は予測性能が十分に高いと言える。

なお、ここまでの結果では、学習データに対する各種性能値とテストデータに対する値にほぼ差が無いが、これは学習データのサンプル数が1000と十分に多く、偏りが無いからである。これを100にすると、次のように、テストデータに対する成績が少し下がるが、それでも、学習データだけに対して大幅に成績が良い「過学習(overfitting)」の状態と言える程ではない。

学習データ数が100の1000試行の平均

Average with training data
Class FamilyOut=False: Precision=0.931, Recall=0.963, F-measure=0.945
Class FamilyOut=True : Precision=0.793, Recall=0.582, F-measure=0.650
Average with test data
Class FamilyOut=False: Precision=0.926, Recall=0.955, F-measure=0.938
Class FamilyOut=True : Precision=0.767, Recall=0.564, F-measure=0.620
Average of AUC with training data=0.899
Average of AUC with test data=0.890

続きを読む "scikit-learnのナイーブベイズ分類器を使ってみる"

2017年05月29日

 Wekaでベイジアンネットワークの事後確率を計算

ベイジアンネットワークの復習をしていて、確率的グラフィカルモデル‐ベイジアンネットワークとその周辺‐(オペレーションズ・リサーチ 2013年4月号)という記事を見つけた。その中に、Family Out Problemという有名な例題(下図)の紹介があり、p.194に、

表1~3で与えられたCPTの値と式(11)を利用して丹念に計算することにより P(X1=1|X3=1,X5=1)=0.7577···という結論を得る.
と書かれていたので、これを自力で計算できたらベイジアンネットワークを理解できたことにしよう、と思って、丹念に計算してみたら、その値にならなかった。


Family Out Problem

ここでは、
X1: Family Out
X2: Bowel Problem
X3: Light On
X4: Dog Out
X5: Hear Bark
(それぞれ2値の確率変数)であり、それぞれの条件付き確率は次の通りである。
P(X1=1) = 0.15
P(X2=1) = 0.01
P(X3=1 | X1=0) = 0.05
P(X3=1 | X1=1) = 0.6
P(X4=1 | X1=0, X2=0) = 0.3
P(X4=1 | X1=0, X2=1) = 0.97
P(X4=1 | X1=1, X2=0) = 0.9
P(X4=1 | X1=1, X2=1) = 0.99
P(X5=1 | X4=0) = 0.01
P(X5=1 | X4=1) = 0.7

P(X_1=1 | X_3=1, X_5=1) = \frac{\sum_{X_2}\sum_{X_4} P(X_1=1, X_2, X_3=1, X_4, X_5=1)}{\sum_{X_1}\sum_{X_2}\sum_{X_4} P(X_1, X_2, X_3=1, X_4, X_5=1)}
なので、Pythonで

P00101 = 0.85 * 0.99 * 0.05 * 0.7 * 0.01
P00111 = 0.85 * 0.99 * 0.05 * 0.3 * 0.7
P01101 = 0.85 * 0.01 * 0.05 * 0.03 * 0.01
P01111 = 0.85 * 0.01 * 0.05 * 0.97 * 0.7
P10101 = 0.15 * 0.99 * 0.6 * 0.1 * 0.01
P10111 = 0.15 * 0.99 * 0.6 * 0.9 * 0.7
P11101 = 0.15 * 0.01 * 0.6 * 0.01 * 0.01
P11111 = 0.15 * 0.01 * 0.6 * 0.99 * 0.7

P35 = P00101 + P00111 + P01101 + P01111 + P10101 + P10111 + P11101 + P11111
P135 = P10101 + P10111 + P11101 + P11111
P1_35 = P135 / P35
print(P1_35)
とすると、0.8578...という数値になった。どこか読み間違えたかと思って、
def prob(x1, x2, x3, x4, x5):
    p = 1.0
    p *= (0.85, 0.15)[x1]
    p *= (0.99, 0.01)[x2]
    p *= ((0.95, 0.05), (0.4, 0.6))[x1][x3]
    p *= (((0.7, 0.3), (0.03, 0.97)), ((0.1, 0.9), (0.01, 0.99)))[x1][x2][x4]
    p *= ((0.99, 0.01), (0.3, 0.7))[x4][x5]
    return p

P135 = sum([prob(1, x2, 1, x4, 1)
             for x2 in range(2)
             for x4 in range(2)])
P35 =  sum([prob(x1, x2, 1, x4, 1)
             for x1 in range(2)
             for x2 in range(2)
             for x4 in range(2)])
P1_35 = P135 / P35
print(P1_35)
と書いてみたが、やはり0.8578...だった。

上記の記事内の条件付き確率表に誤記があり、0.7577...というのは元の確率表で計算された値か、とも思ったが、どの文書のFamily Out Problemを見ても確率は同じだった。

正解は何なのか、何らかのツールで確認しようと思って、Pythonのベイジアンネットワーク関連のツールを探したが、適当なものがなかなか見つからなかった。

scikit-learnには"Naive Bayes"のAPIはあるがベイジアンネットワークは見つからない。
PyMCBayesPyは、きっとうまく使えばこの計算ができるのだろうが、ネットワークを定義して、CPT(conditional probability table、条件付き確率表)とevidence(観測値)を与えて事後確率を計算する直接的なサンプルコードが見つからなかったので、諦めた。
PBNTにはそういうサンプルコードがあったので、使ってみたが、P(X1=1|X3=1,X5=1)=0.2018...という全然違う値が出力された。P(X5=1)=0.2831と、これは正しい値が出たので、ネットワークとCPTは合ってそうであり、事後確率を計算するにはengine.marginal()でなく別のメソッドを使わないといけないのかとも思ったが、よくわからなかった。

Pythonを諦めてツールを探すと、Wekaでできることがわかった。Wekaは機械学習の勉強をするなら必修らしく、過去にインストールしていたので、やってみた。
Wekaを起動して、Tools->Bayes net editorを開くと、GUIがバグだらけ(Version 3.8.1, WindowsとMacとで確認)で使いにくいが、ノードを追加して右クリックしながらネットワークを作成し、Tools->Layoutでノードの配置を修正し、CPTを設定し、さらに保存したXMLを書き換えて色々修正し、Tools->Show Marginsを選ぶと、次のように結合確率が表示される。

さらに、右クリック->Set evidenceで LightOn=True, HearBark=True と設定すると、次のように、各ノードの事後確率が表示される。

これによると、FamilyOut=Trueの事後確率はやはり0.8578...である。筆者は何か問題を読み間違えているのだろうか?


2017年05月21日

 角交換振り飛車に大ポカ一発で沈む

昨日は毎年参加している地域の将棋大会だった。
今年は2回戦敗退だったが、予選は1勝1敗で抽選で勝ち抜け(3人のブロックで3人とも1勝1敗だった)、抽選で決勝1回戦はシードだったので、昨年の1回戦敗退より悪い内容だった。通常3勝かかる所、たった1勝で2回戦まで進出するとは、何とくじ運の良かったことか。

筆者はここ1年くらい将棋の勉強をしておらず、全て忘れてしまっており、先月からたまにネット将棋を指していたが以前のレートでは全く勝てず、今年は予選突破できないと思った。
昔から記憶力が無い方であるが、たった1年休んだくらいで全て忘れてしまい、レートが200も下がるのでは、筆者は将棋に向いてないのだろうとつくづく思う。

2回戦の相手は昨年優勝のI藤さんだった。筆者が絶好調でもまず勝てない相手であるが、先手の筆者に次の局面のような感触の良い手が出た。

通常は4二の銀が4四に居るので、3四の歩を取っても響かないが、この形だと歩を取った後に3筋の歩を伸ばせる。3四の歩を受けるには△1二角と打つしかないが、形が悪いだろう。
筆者は通常、角交換振り飛車には▲6六歩のようにして角交換を拒否するのだが、1回戦で筆者が100回やっても勝てそうにないY田さんがI藤さんに対して▲6六歩と角交換を拒否して負けたのを見て、直前に何か別の展開をガラケーの自作アプリに仕込んだ棋譜の中から1つ探して、その通りに指してみたものである。そういうことをすると通常は相手の研究にはまるので、ろくなことにならないのだが、今回は功を奏した。
対して、後手のIさんは△3二飛と指した。

これを見た瞬間、チャンスだと思った。2三の地点が空いている。
その誘惑に駆られた上、一手前に長考したこともあり、じっくり考えようとは思わなかった。
▲3四角に△4五桂と両取りに跳ねられても▲同角で両方受かる、と、それ以上考えずに、10秒も使わずに▲3四角と指したら、△2五桂とされて、一発で撃沈してしまった。

あと5秒考えていればこの手に気付いただろう。有段者が何の意味もなく、△3二飛のような隙を作る手を指すはずがない、と、何故一瞬でも思わなかったのか。
△2五桂の後、▲同歩、△3四飛、▲2四歩とすれば、後手は歩切れで、と金ができて大差にはならなさそうだが、この局面では丁度△1五角があって、受かってしまう。

△3二飛の局面は激指14のPro+3でも評価値が+280であり、少しリードしていたのは間違いない。その後、▲3四角ではなく、▲2九飛 △5四歩 ▲4八金と桂馬に紐を付け、△1二角 ▲6六歩 △6四歩 ▲1六歩のように進めれば、まあまあ指しやすそうである。

ここまでの手順は激指14のPro+以上の手を続けたものであり、評価値は+228である。

まあ、相手が相手なだけに、こう進んでもまず勝てなかっただろうが…

続きを読む "角交換振り飛車に大ポカ一発で沈む"

 MacPortsでmaximaをインストールするとエラーになるのを回避

現在、MacPortsを使って

sudo port install maxima
すると、下のようなエラーになって失敗してしまう。(macOS Sierra 10.12.4で確認)
../../doc/info//errormessages.texi:296: `Warning messages' has no Up field (perhaps incorrect sectioning?).
../../doc/info//errormessages.texi:162: Prev field of node `Operators of arguments must all be the same' not pointed to.
../../doc/info//errormessages.texi:152: This node (Only symbols can be bound) has the bad Next.
../../doc/info//errormessages.texi:152: Next field of node `Only symbols can be bound' not pointed to (perhaps incorrect sectioning?).
../../doc/info//errormessages.texi:204: This node (out of memory) has the bad Prev.
../../doc/info//errormessages.texi:8: `Error messages' has no Up field (perhaps incorrect sectioning?).
makeinfo: Removing output file `/opt/local/var/macports/build/_opt_local_var_macports_sources_rsync.macports.org_macports_release_tarballs_ports_math_maxima/maxima/work/maxima-5.39.0/doc/info/maxima.info' due to errors; use --force to preserve.
make[3]: *** [maxima.info] Error 1

原因は、makeinfoのバージョンが古いからのようだ。
http://maxima-discuss.narkive.com/KeLR9qhb/can-t-install-5-39
に、makeinfoのバージョンが4.xだとこうなることが書かれている。 実際、macOS Sierraで makeinfo --version すると、(GNU texinfo) 4.8と表示される。

そこで、

sudo port install texinfo
(GNU texinfo 6.3のmakeinfoがインストールされる)してから再度
sudo port install maxima
すると、上記のエラーは出なくなった。Activationの段階で
--->  Installing maxima @5.39.0_3+xmaxima
--->  Activating maxima @5.39.0_3+xmaxima
Error: Failed to activate maxima: Image error: /opt/local/bin/maxima already exists and does not belong to a registered port.  Unable to activate port maxima. Use 'port -f activate maxima' to force the activation.
というエラーになるが、メッセージの通り、
sudo port -f activate maxima
とすると、一応wxMaximaと併用しても問題なく動いた。

続きを読む "MacPortsでmaximaをインストールするとエラーになるのを回避"

2017年05月06日

 JediにAutoCompleよりCompanyModeを優先させる

2ヶ月前にCarbon EmacsのPython環境を整備したばかりなのだが、MacのOSを10.7.5から10.12.4にバージョンアップすると、Carbon Emacsがまともに動かなくなってしまった。起動はするのだが、すぐに固まってしまう。しかも、その後distnotedというプロセスが全てのCPUを奪い続けて、Macが激重になる。
この現象はOS X 10.9+Emacs 24.3で起こるらしく、Carbon Emacs(Emacs 22)でも起こったという情報は見つけられなかったが、おそらく同じ問題である。どちらかと言うとOSのバグなのだが、現時点ではEmacsを24.4以降にバージョンアップする以外に解決方法が見当たらない。

その為、長年お世話になったCarbon Emacsを手放し、Emacsの最新版である25.2を使うことにした。

EmacsWikiのPython Programming In EmacsのページにはPythonの開発環境を便利にする方法が色々書かれているが、筆者は高機能なPythonの開発環境を使いたい時はSpyderを使っており、EmacsのPython環境はPython.elでPython3が使えて、もう少し便利な補完が効けば十分なので、Jediだけを追加インストールすることにした。

Emacs 25.2にbuilt-inのPython.elはPython3に対応しているのだが、MacPortsでpython35をインストールして、M-x customize-group -> pythonして"Python Shell Interpreter"をpython3にすると、Warningが出まくり、妙に不安定だった。M-x list-packagesしてpython packageを0.25.2にバージョンアップすると少し改善したが、完全には直らなかった。
さらに、Jediをインストールすると、これまた不安定で、時々補完候補が出る前にEmacsが固まった。(C-gで抜けることはできる)

Emacs 24の最新版である24.5.1ではこの問題は起こらなかったので、当面、EmacsでPythonのコードを書く時は24.5.1を使うことにした。

Jediのインストールは、

(require 'package)
(add-to-list 'package-archives '("melpa" . "https://melpa.org/packages/"))
(package-initialize)
としてpackage.elにMELPAのリポジトリを登録し、M-x list-packagesしてjediをインストールし、Jedi.elのドキュメントに従ってセットアップした。

なお、Emacs25でjediをインストールすると、Emacs24で"Symbol's function is void: cl-struct-define"というエラーになった。これはjediに限らず、よくあることだそうで、Emacsの24と25を併用する場合、package.elで何かをインストールするならEmacs24でする方が良さそうだ。

さて、JediはAutoCompleteとCompanyModeの両方に対応しているが、package.elで"jedi"をインストールするとAutoCompleteが使われる。"company-jedi"をインストールするとCompanyModeが使われるのだが、"jedi"と"company-jedi"の両方をインストールすると、AutoCompleteが優先される。特にAutoCompleteに不満は無く、和製なので贔屓したいが、色々読んでいると、世界的には、日本でも現在は、CompanyModeの方が人気があるようなので、何が良いのかを知るために、これからしばらくはCompanyModeを積極的に使うべく、両方をインストールしてCompanyModeをデフォルトにすることに決めた。

しかし、.emacsや.emacs.d/init.elでCompanyModeを優先する適当な方法がわからなかった。Python modeにしてM-x auto-complete-modeとすればAutoCompleteのON/OFFが切り替わり、OFFだとCompanyModeが使われるのだが、Emacs Lispで(auto-complete-mode)としてもON/OFFが切り替わらず、auto-complete-mode関数にはOFFにする引数も無いのである。auto-complete-mode関数のソースコードを見ると、"(if auto-complete-mode ..."とあるので、

(setq auto-complete-mode nil)
(auto-complete-mode)
とすれば良さそうに思ったが、これでもOFFにならない。
但し、(setq auto-complete-mode nil)すると、即座にAutoCompleteがOFFになる。このことを使って、試行錯誤の末、.emacs等に
(add-hook 'python-mode-hook 'jedi:setup)
よりも前のどこかに
(add-hook 'python-mode-hook
(lambda () (setq auto-complete-mode nil)))
と書けば、Python mode開始時にAutoCompleteがOFFになり、CompanyModeが使われることがわかった。

続きを読む "JediにAutoCompleよりCompanyModeを優先させる"