Pythonでパーセプトロンのふるまいを可視化する
普段分析で機械学習は使いますが、実際に学習器がどう動いて回帰や分類をするのか理屈では理解していた（つもり）ですが、体感としてはよくわかってませんでした。
それに機械学習パッケージを使うと、恐ろしいことに仕組みを全く知らなくてもそれらしい結果が出せてしまうので、本質を理解できていなくても使いこなしている気分になってしまいます。
まあ使えて結果が出れば、それでいいじゃんという話もありますが、なんとなくもやもやしてたんですよね。
で、そんなときに以下の記事を見つけて、すごく肚に落ちたというかもやもやがかなりスッキリしました。
ロジスティック回帰 (勾配降下法 / 確率的勾配降下法) を可視化する - StatsFragments
やっぱり可視化されると納得感が違いますねー。
ということで、自分でもやってみたいと思いこの記事を書きました。
（sinhrksさんの許可を得て一部コードは流用させていただきました。ありがとうございます。）
前置きはこのくらいで、本題に入っていきます。
パーセプトロンについて
パーセプトロンが何かについては、検索すればいくらでも記事や書籍がありますのでここで詳細な説明は省きます。
例えば私は以下の記事などを見て勉強させていただきました。
TokyoNLP#5で「パーセプトロンで楽しい仲間がぽぽぽぽ〜ん」を発表しました - 睡眠不足？！
機械学習超入門III 〜機械学習の基礎、パーセプトロンを30分で作って学ぶ〜 - EchizenBlog-Zwei
テキストマイニングのための機械学習超入門 二夜目 パーセプトロン - あんちべ！
機械学習モデルの分類の１つとして、事前分布から仮定する生成モデル・事後分布のみ仮定する識別モデル・分布を仮定しない識別関数という見方*1がありますが、パーセプトロンは識別関数の最も基本的な学習器と言えます。
最近の機械学習モデルの学習ルールと言えば、損失関数を最小化するようにパラメータを更新していきますが、パーセプトロンは間違えた時だけ更新する（誤り訂正法）ところが男らしいといえば男らしい気がしています(?)。
データを準備
まずはデータを作っていきます。
irisデータを使って、デモ分析用に少し加工します。
def load_data(): # データ準備 iris = pandas.read_csv("http://aima.cs.berkeley.edu/data/iris.csv", header=None) names = ['Sepal.Length', 'Sepal.Width', 'Petal.Length', 'Petal.Width', 'Species'] iris.columns = names # ２説明変数、目的変数（２クラス）に加工 data = iris[:100] x = data[['Petal.Width', 'Petal.Length']] y = data['Species'] # ラベルを-1,1の列に変換 y[ y=='setosa' ] = 1 y[ y=='versicolor'] = -1 # ランダムシャッフル indexer = numpy.arange(x.shape[0]) numpy.random.shuffle(indexer) # x,y並び替え y = y.iloc[indexer, ] x = x.iloc[indexer, ] # indexを初期化 x = x.set_index([range(100)]) y.index = [range(100)] return x,y
合わせて、可視化用の関数も作ります。
def plot_xy(x, y, colors, ax=None, figsize=(5,3.5)): if ax is None : # 描画領域作成 fig = pyplot.figure(figsize=figsize) # 描画領域にaxesを追加、マージン調整 ax = fig.add_subplot(1, 1, 1) fig.subplots_adjust(bottom=0.15) x1 = x.columns[0] x2 = x.columns[1] for (species, group), c in zip(x.groupby(y), colors): ax = group.plot(kind='scatter', x=x1, y=x2, color=c, ax=ax, figsize=figsize) return ax
デモ分析用のデータをplot_xyで可視化すると、こんな感じです。普通な感じですね。
探索法
決定境界を探索するため、xを元にyの予測値を定義していきます。
まずはxをwx+bという形でyと同じ次元に写像します。その後、さらに[-1,1]という離散値を持つ値に変換します。
なぜこんなシンプルな変換で良いのかというと、パーセプトロンは与えられた問題が線形分離可能であることを前提として構築されているからですね。なので「実際のyと予測値のyが違った時だけ、パラメータを少し更新する」という戦略で更新と探索を繰り返していけば、いつかは必ず決定境界が見つかることになります。よってyの予測値は「当たっているか外れているか」だけわかればよく、どの程度当たったか外したかという度合い（距離）は不要なので、すごくシンプルな写像になります。
これを表現してyの予測値を求めるのが以下のコードです。
def yhat(x, w, b): out = numpy.dot(x, w) + b def link(a): if a >= 0: res = 1 else : res = -1 return res return link(out)
学習ルール
すでにくどくどと書いてしまっていますが、パーセプトロンでは予測が違った時のみ更新する「誤り訂正法」でパラメータを更新していきます。どの程度更新するか（学習率）はケースバイケースですが、あまり大きく更新すると振動して収束しなくなることもあるようです。今回は学習率として定数0.2を設定しました。更新の方法としては勾配降下法を用いています。
def learning_rule(x, y, w, b, eta=0.2): #パーセプトロン学習規則 res = yhat(x, w, b) if ( res*y )<0: w = w + eta*x*y b = b + eta*y return w, b
上記学習を全データに対して実行する関数を作ります。
def update(x, y, w, b, num=3): for i in range(1, num): for index in range(x.shape[0]): # 1行ずつ逐次処理 _x = numpy.array(x.ix[index, :]) _y = numpy.array(y.ix[index, :]) w,b = learning_rule(_x, _y, w, b) yield i, w, b, _x
さらに、実行結果を逐次グラフとして保存するように上記関数をラップします。
def run_perceptron(x, generator, fig): figs = [] def db_2d(x, w, b): x2 = -w[0]/w[1]*x - b/w[1] #決定境界 return x2 for i, w, b, _x in generator: by = db_2d(x, w, b) _fig = fig.plot(x, by, color='gray', linestyle='dashed') wt = """Iteration = {0} times w = [{1[0]:.2f}, {1[1]:.2f}] b = {2:.2f}""".format(i, w, b) # axes上の相対座標(0.1, 0.9)にtextの上部を合わせて描画 tx = fig.text(0.1, 0.9, wt, va='top', transform=fig.transAxes) # プロトタイプを描画 _p = fig.plot(_x[0], _x[1], 'ks') figs.append(tuple(_fig) + tuple(_p) + (tx, )) return figs
これで必要な道具は揃いました。
これらを元に実行スクリプトを書くと、以下のようになります。
# 描画準備 ########### figsize =(5, 3.5) fig = pyplot.figure(figsize = figsize) ax = fig.add_subplot(1, 1, 1) fig.subplots_adjust(bottom=0.15) # 学習準備 ########### # データ作成 x, y = load_data() # w, bの初期値を作成 w, b = [0,0], 0.1 gen = update(x, y, w, b) # 学習開始 ########### bx = numpy.arange(x.ix[:, 0].min(), x.ix[:, 0].max(), 0.1) figs = run_perceptron(bx, gen, ax) # 描画開始 ########### ax = plot_xy(x, y, colors=['red', 'blue'], ax=ax) ani = animation.ArtistAnimation(fig, figs, interval=3, repeat=False) pyplot.show()
上記コードを実行することで、以下のようにパーセプトロンの動きが可視化されました！
元データがきれいなので、比較的簡単に決定境界を見つけてしまいますね。。
一応、その時点で評価されているxの座標を黒い四角で表しています。
探索中に、この評価点の決定境界から得られる推定値が異なっている（誤りがある）時点で、パラメータが更新されて決定境界が動く様子が見て取れます。
今回の記事は以上です！
もしなにか間違いがあれば、ご指摘いただけると助かります。
せば
Pythonでジニ係数を計算してみる
データ分析で予測分類モデルとか作るときに、目的変数に対するカテゴリカルな特徴量の寄与をただ知りたいときってありますよね。
Rを使ってたときは、パッケージでジニ係数とか情報利得とかを計算する関数があったんですが、Pythonでは決定木とかに組み込まれているものしか見つけられなかったので作ってみました。
以下の流れでご紹介
ジニ係数って
wikiによれば
ジニ係数(ジニけいすう、Gini coefficient)とは、主に社会における所得分配の不平等さを測る指標。ローレンツ曲線をもとに、1936年、イタリアの統計学者コッラド・ジニによって考案された。所得分配の不平等さ以外にも、富の偏在性やエネルギー消費における不平等さなどに応用される。
とのことですが、データ分析の文脈では分布の不均衡さを表す指標として使います。
Pythonで実装
まずは単一の状態でのジニ係数を計算する関数を実装します。
こいつをある変数に対して分割前の状態と分割後の各状態で計算して、その差分を見るのがデータ分析でのジニ係数の扱いですね。
# -*- encoding=utf-8 -*- import numpy import pandas from collections import Counter def gini(vec): prob2s = list() count = Counter(vec) countall = float(numpy.sum(count.values())) for item in count.items(): counteach = item[1] prob = counteach/countall prob2 = prob**2 prob2s.append(prob2) gini = 1 - numpy.sum(prob2s) print "%s : %f" % (vec.name, gini) return gini
続いて、このginiを呼び出して分割前の目的変数yのジニ係数と、あるカテゴリカル変数xで分割した後のジニ係数を比較します。この差分が正の方向に大きいほど、「xによる分割で各サブセット内のyが均質化した」ことになるので、xは影響力のありそうな特徴量かな、となります。
def giniIndex(x,y): # x,yともにカテゴリカル変数 root_gini = gini(y) grouped = y.groupby(x) leaf_gini_0 = grouped.apply(gini) leaf_weight = grouped.apply(len) / float(len(y)) leaf_gini = leaf_gini_0 * leaf_weight # ０に近いほど均質 = 正例・負例の純度が高いカテゴリが存在 giniindex = root_gini - leaf_gini.sum() print "gini index : %f" %giniindex return giniindex
ついでにxが連続量の場合についても作ってみました。
厳密にはxについて積分が必要な気がしますが、簡易版として100分割して
離散化しています。
def giniIndexNum(x,y, bins=100): # xが連続変数、yがカテゴリカル変数 root_gini = gini(y) xname = x.name yname = y.name xy = pandas.concat([x,y], axis=1) xbins = pandas.cut(xy[xname], bins) grouped = xy[yname].groupby(xbins) leaf_gini_0 = grouped.apply(gini) leaf_weight = grouped.apply(len) / float(len(y)) leaf_gini = leaf_gini_0 * leaf_weight giniindex = root_gini - leaf_gini.sum() print "gini index : %f" %giniindex return giniindex
では実際に使ってみます。*1
x = numpy.random.choice(["A","B","C"], 100) x = pandas.Series(x) y = numpy.random.choice([0,1], 100) y = pandas.Series(y) y.name = "testy" giniIndex(x,y)
これを実行すると、以下のようになります。
testy : 0.499800 A : 0.482422 B : 0.495317 C : 0.499635 gini index : 0.007012
せば
*1:12/26 コードに一部誤りがあったので修正しました。m(_ _)m
