最尤推定・交差エントロピー誤差（実装編）
今回は3回に分けて学んだ分類問題の実装編です。
以前の記事をまだ読んでいない方はあらかじめ目を通しておくことをおすすめします。
まずは手順の再確認です。
実装に躓いたら上記の図や記事を参考にして、基本に立ち返ります。
1.モジュールのimport
以下のモジュールを使って実装します。
import matplotlib.pyplot as plt import numpy as np import math import copy import random %matplotlib inline
2.データの作成
2次元配列でデータを作成。[体重データ、性別データ]のセットで作成します。
#データの作成 x=np.array([]) y=np.array([]) for i in range(20): if i<=9: sex=0#random.randint(0,1) y=np.insert(y,len(y),sex,axis=0) weight=np.random.normal(60.0 , 5.0) x=np.insert(x,len(x),weight,axis=0) else: sex=1#random.randint(0,1) y=np.insert(y,len(y),sex,axis=0) weight=np.random.normal(70.0 , 5.0) x=np.insert(x,len(x),weight,axis=0) data=np.c_[x,y] print(data)
3.データの可視化
今回は女性と男性の体重データがそれぞれ10個ずつの想定で作ってます。
#データの可視化 plt.scatter(data[:9,0], data[:9,1],color="red") plt.scatter(data[10:19,0], data[10:19,1],color="blue") plt.show()
4.パラメータと関数の定義
4.1.学習パラメータ
param = np.array([1.0 for i in range(2)])
4.2.シグモイド関数
def sigmoid(x): e = math.e s = 1.0 / (1.0 + e**-x) return s
4.3.尤度関数
戻り値は尤度の符号を変えているので交差エントロピーに変換しています。
def likelihood(param,data): mean = sum(data[:,0])/len(data) data_copy=copy.deepcopy(data) data_copy[:,0]-=mean y=np.array([0.0 for i in range(len(data_copy))]) #直線の方程式 y = param[0]*data_copy[:,0]+param[1] #シグモイド関数へ代入 y = sigmoid(y[:]) prob = 1.0 for i in range(len(data_copy)): prob *= y[i]**data_copy[i,1]*(1- y[i])**(1-data_copy[i,1]) return prob*(-1)
4.4.微分関数
下記の関数は線基底関数モデル（実装編）でもご紹介しました。
def DER(param,func,data,index): epsilon=0.0001 p1=[0.0 for i in range(len(param))] p1+=param p1[index]+=epsilon p2=[0.0 for i in range(len(param))] p2+=param p2[index]-=epsilon mse=func(p1,data)-func(p2,data) der=mse/(2*epsilon) return der
4.5.学習部分
微分関数と同様以前ご紹介しました。
def fitting(param,data): alpha=0.5 for i in range(len(param)): der=DER(param,likelihood,data,i) param[i]-=der*alpha
5.学習結果の可視化
mean = sum(data[:,0])/len(data) x = np.arange(mean-20,mean+20,0.01) #もっと見やすく改良、正規化されたデータを元に戻して表示↓ plt.plot(x[:], sigmoid(param[0]*(x[:]-mean)+param[1])) plt.scatter(data[:9,0], data[:9,1],color="red") plt.scatter(data[10:19,0], data[10:19,1],color="blue") plt.show()
実際の学習結果の例↓
縦に赤と青の点が並んでいる場所では１と０の間に曲線があり、そうでない場所では１か０に曲線が振り切れているのがわかります。
このように確率のアイデアを取り入れることで複数の種類のデータを曖昧さも加味しながら分類することができました。
最尤推定・交差エントロピー誤差(2)
前回：最尤推定・交差エントロピー誤差(1)に引き続き進めていきます。
まだご覧になっていない方は、前回記事から見ることをオススメします。
今回は今までやってきたことの総まとめです。分類問題の最後ですので頑張りましょう。
上記の図は今まで学んできた分類の処理を簡単にまとめたものです。
順を追って見直していきます。
(1)シグモイド関数の入力にする。(2)シグモイド関数の決定
ここでは直線の式（１次関数）をシグモイド関数の入力として用いることで、まるで直線の両端を押しつぶしたようなグラフが得られることが分かりました。
直線のパラメータ（a,b）を求めることで、尤もらしい確率のグラフが完成します。
(3)尤度の計算
前提としてyは「与えられたデータが男性の体重である確率」を表しており、tは「実際に体重データが男性のものだったのかどうか」を表します。
（※yは(1)(2)で作成したシグモイド関数のこと）
プログラムの実装を考えると与えられたデータが「男性だった場合」を「女性だった場合」で条件分岐式（if文）を書きたくなるところですが、数学の力でうまく表現してやりましょう。
L= y ^ t ( 1 - y ) ^ ( 1 - t )…上記の図にも同じ数式があります。
このような表現をすることで条件分岐式を使わずともうまく計算できます。
男性データ（識別子：１）を扱う時は( 1 - y ) ^ ( 1 - t )の部分が１になります。
女性データ（識別子：０）を扱う時はy ^ t の部分が１になり掛け算されます。
ここで分類問題とは別に学べることは、わざわざ条件分岐式のようなプログラミング記法を頼らなくても、数学の力を使うとうまく処理できることがあるということです。
地道で単調なコーディングをすることは時として重要ですが、物事の構造に迫り活用することもプログラミングの醍醐味だなと感じます。
(4)交差エントロピー誤差に直してパラメータの導出
パラメータの導出のために最尤推定法としてここまで話を進めてきました。
つまり尤度の最大化問題を考えてきたわけです。
しかし、以前学んだ勾配降下法を使うためには値の最小化問題を考えたいので尤度に-1をかけてしまいます。
すると、パラメータ（a,b）の推定を行うことができます。
さて、今回は2つのものに分類する問題について考えてきました。
分類（実装編）の記事は今まで3回に分けてお伝えしてきたことを1つの記事にまとめて掲載しようと考えています。
最尤推定・交差エントロピー誤差(1)
前回：確率の導入では分類の問題に確率を用いるメリットと、シグモイド関数について解説してきました。
今回からはそのシグモイド関数の形の決定について扱っていきましょう！
少し長いので2回に分けて学んでいきます！
問題：体重データから性別を推定する。
解決したい問題の例は前回同様です。
今回は前回勉強したシグモイド関数を色々な形に変形してみるところからスタートしてみましょう！
シグモイド関数の変形
以前も紹介しましたが、シグモイド関数の式は以下の通りです。
（※wikipediaから引用）
この関数の引数に様々なパターンの値を代入して形がどう変化するのか見てみましょう。
左上：ε(x) 標準シグモイド関数
左下：ε(x-20) 右に20平行移動
右上：ε(-0.3*x) 形を滑らかに変形・左右を反転
右下：ε(0.15x+5) 上記３つの応用
ただデータ(x)を用いるのでなく、１次関数(ax+b)になるよう加工してからシグモイド関数へ代入すると、変形できることがわかります。
「形が自由自在に操れる」ということは、「どんなデータセットに対しても適切な形のシグモイド関数を求められる」ということになります。
すると、どんなデータセットが与えられたとしても確率が正確に表せるシステム構築へと繋がるわけです。
最尤推定
さて、ここで確率の求め方の紹介を挟みます。
上図の緑色の線で囲まれた部分に注目してください。
青点が２個、赤点が４個入っているのが確認できると思います。
では、問題です。
この緑色の線の中から一つデータを選ぶとき、男性のデータが選ばれる確率はいくらでしょう？
もちろん全部で６個あるデータのうち、２個が男性のデータなので確率は2/6（33％）が正解ですよね。
この方法は小学生の時から使っている分かりやすい解析的な確率導出の手段です。
しかし、以下２点の理由からオススメの方法とは言えません。
- 機械学習では逐次教師データが更新される場合があるため、その都度数え直す手間が大きい。
- 男性データの個数を数えるのに条件分岐式を使うことが予想される。（if文などの条件分岐式は多用すると実行に時間が多くかかる恐れがある。）
そこで出てくるのが最尤推定法です！
先ほど緑色の線で囲まれたデータが、男性のものである確率は33%と求めました。
その確率が尤もらしい（もっともらしい）、つまり「それっぽさ」の度合いを表す言葉が尤度です。
求め方は男性のデータ２個（各33%）、女性のデータ４個（各66%）の場合なので…
0.33 × 0.33 × 0.66 × 0.66 × 0.66 × 0.66 = 0.0206… となります。
0.33を2回、0.66を4回掛け算すれば良いです。
実はこの尤度という値は、確率（今回でいうと男性のデータが選ばれる確率の値）が正しく導かれているときに最大値をとります。
（よって男性33％、女性66％のとき最大値）
本当に他の値の組み合わせで最大値を取らないか検証します。
男性のデータが得られる確率が10%、女性のデータが得られる確率が90%じゃないかと予想して尤度の計算をしてみます。
0.1 × 0.1 × 0.9 × 0.9 × 0.9 × 0.9 = 0.0065… となりました。
0.9という数字を使っていることから数値が大きくなるかと思いきや、結果はやはり先ほどの男性33%、女性66%の組み合わせの方が圧倒的に尤度が大きいことが分かります。
このように確率の「それっぽさ」は尤度の最大化計算によってなされます。
さて、以前：勾配降下法（理論編）にて、ある関数が最小値をとるようなパラメータの導出法を学びました。
今回の記事で学んだ尤度は最大化することを目標としていますが、-1をかけてやることで最小化問題に変化します。このように尤度の符号を変化させたものを交差エントロピー誤差と呼びます。
交差エントロピー誤差を使うことで以前学んだ勾配降下法の学習プロセスをそのまま使うことができます。
まとめ
次回はシグモイド関数で表された確率の尤度について学びます。
確率の導入
前回まではデータの散らばりを「それっぽい直線」を考えることによって学習する回帰を扱ってきました。
今回からはあるデータ群を複数のグループ（クラスタとも呼ぶ）に分類する問題に取り組みましょう。そのためには確率の考え方が必要不可欠になってきます。
まずは問題を参考にしながらなぜ確率の考え方が必要なのかという点から考えていきましょう。
問題：体重データから性別を推定する。
例えば、下図のような20個のデータがあったとしましょう。
横軸は体重を表していますが、縦軸には今は意味がありません。
実はこのデータは男女10人ずつから収集したものと仮定します。
しかし、このままではどのデータが女性のもので、どのデータが男性のものか分かりませんね。なので分かりやすいよう男女を色別に分け、下図のようにしてみます。
これら20個のデータは教師データとして扱い、21個目のデータが男性なのか女性なのか推論していきます。
（※男性を青、女性を赤としてプロットした図）
さて、このデータを元に男性・女性を分類するための「判断基準」を手に入れることが学習段階では必要です。では、今までのように線形回帰で問題を解くことはできるでしょうか？
答えはNOです。
線形回帰で扱っていた問題は体重のデータ（連続数）をもとに身長の値（連続数）の推定を行ってきました。
しかし、今回推定する問題は体重のデータ（連続数）をもとに性別（離散数）の推定を行います。
取り扱う数字の性質が異なるので、同じ方法で推論しようとしてもうまくいきません。
そこで、方針としては体重のデータから男性と女性との変わり目を見つけ、境界線を引くのが分かりやすそうです。下図でいうと、体重60~65の間でだいたい性別が分かれているように見えます。
この境界線を決定境界と呼びます。
このままでも良いのですが、ハッキリと2分したことで左側に青点、右側に赤点が紛れ込んでいますね。つまり決定境界周辺は、精度があまり高くないことが分かります。
この多分合ってるんじゃないかな？という不確定さ（もしくは曖昧さ）を数学的に表現できるとシステムの精度を客観的に測ることができますね。
そこで使うのが確率の知識です。
ここで、もっと図を見やすくするために性別のデータも交え、しっかりと教師データを見つめてみます。今回は男性を表す値として１、女性を表す値として０を与えています。
ちなみに性別データは男性を１女性を０として値を割り振ることはできますが、割り振った値の数学的な大小関係や性質を使って処理することはできません。
（無論、割り振る数値は男性と女性で逆になっていても構いません。）
※今回使う教師データをプロットしたもの。体重と性別の2次元データになっている。
シグモイド関数
ここで確率を表現するのに便利な関数を紹介します。
確率は基本的に0から1までの実数で表されるため、特殊な形の関数が必要です。
それがシグモイド関数です。これは漸近線をy=1とy=0の２本とる関数になっているため、確率を表すのに非常に都合の良い関数になっています。
式で表すと以下のようになります。
（※wikipediaから引用）
実際にpythonで描画したもの↓
ちなみに、式の内部にあるαはゲインと呼ばれ、ゲインが１のシグモイド関数を標準シグモイド関数と呼びます。
今後シグモイド関数を用いる際は、基本的にこの標準シグモイド関数を使っていきます。
このシグモイド関数をデータ点と重ねて表示すると下図のようになります↓
これは「与えられた体重のデータが男性の物である」確率を表していると思ってみてください。
すると、確実に男性の体重しか得られていない65kg以上の範囲ではシグモイド関数はほぼ１に近づいています。
それに対して60kg未満の範囲では女性の体重のみ得られているため、シグモイド関数は０に近づいています。
まとめ
- 分類・グループ化の問題は確率を導入して考える
- シグモイド関数の知識について
今回は確率の概念の必要性と確率を表現するための関数について記載しました。
次回はデータからシグモイド関数の形を評価する手法をご紹介していきたいと思います。
線形基底関数モデル（実装編）
前回：線形基底関数モデル（理論編）で紹介した内容を実際にプログラミングしていきます。
以前学んだ内容もバージョンアップ（より一般化）してもう一度組み直してみます。
では始めて行きましょう↓
線形基底関数モデルの実装
1.データの準備
今回採集された（と仮定する）データを適当に作ります。
データは[x座標,y座標]となるような行列で表します。
def data_log(x): y=np.log(x*10) return y #データの作成 x=np.array([]) y=np.array([]) for i in range(20): tall=random.uniform(0,15) x=np.insert(x,len(x),tall,axis=0) norm=np.random.normal(0 , 0.05) weight=data_log(x[i])+norm y=np.insert(y,len(y),weight,axis=0) data=np.c_[x,y]#データの結合 print(data)
可視化もします↓
#データの可視化 x_axis=np.arange(0,15.0,0.1) plt.plot(x_axis[:],data_log(x_axis[:])) plt.scatter(data[:,0], data[:,1],color="red") plt.show()
2.正規関数の用意
上で準備したデータに正規関数の重ね合わせで学習して行きます。
各正規関数には重みづけをしてから重ね合わせます。この時に使う重みの値が推定するパラメータになります。
#初期パラメータ param=np.array([1.0 for i in range(6)]) #正規関数 def gauss(x,mean,sigma): y=np.exp(-(x-mean)**2/(2*sigma**2)) return y #正規関数の重ね合わせ def mixgauss(x,param): y=0 for i in range(len(param)): y+=param[i]*gauss(x,i*3,5) return y
好きな範囲でプロットして確認↓
x_axis=np.arange(-15.0,15.0,0.1) for i in range(len(param)): plt.plot(x_axis[:],gauss(x_axis[:],i*3,5)) plt.plot(x_axis[:],mixgauss(x_axis[:],param)) plt.show()
3.学習に必要なメソッドの定義
根本の考え方は最小二乗法なので誤差を求める↓
（下記のプログラムはyとzの内積計算をしているが、分かりにくければfor文を使って和を取れば良い。結果は同じ。）
def min2error(param,data): y=0.0 z=[1.0 for i in range(len(data))] y+=(mixgauss(data[:,0],param)-data[:,1])**2 error2=y.dot(z)/len(data) return error2
与えられたインデックス（index）の値を収束させた極限を求める↓
（つまりは偏微分。以下のコードはパラメータの個数が可変長でもうまく動作するはず。）
def DER(param,func,data,index): epsilon=0.0001 p1=[0.0 for i in range(len(param))] p1+=param p1[index]+=epsilon p2=[0.0 for i in range(len(param))] p2+=param p2[index]-=epsilon mse=func(p1,data)-func(p2,data) der=mse/(2*epsilon) return der
全てのパラメータを1回ずつ更新する↓
def fitting(param,data): alpha=0.5 for i in range(len(param)): der=DER(param,min2error,data,i) param[i]-=der*alpha
4.実際に学習してみる
今回は1000回学習する。
for i in range(10000): fitting(param,data) print('学習結果') print(param)
最終確認↓
#データの可視化 x_axis=np.arange(0,15.0,0.1) plt.plot(x_axis[:],mixgauss(x_axis[:],param)) plt.scatter(data[:,0], data[:,1],color="red") plt.show()
なかなか良い形になったのではないでしょうか？
学習回数を増やすとさらに値が変化し、データに沿うような関数が得られます。
ただし、正規関数の平均値と分散を固定、または任意な値で固定していたので学習には限度があるかもしれません。
線形基底関数モデル（理論編）
前回：勾配降下法（理論編）まではデータ点がある直線的な分布で得られると仮定してフィッティングを行ってきました。
しかし、データの中には明らかに直線的な分布にならないデータもあるはずです。
そこで、複雑に曲がったデータの分布を線形回帰で表現するための方法をご紹介します。
線形基底関数モデル
下図のようにデータの分布が曲線を描いているようなものについて考えていきます。
線形回帰ではこのデータの分布にそれっぽい線を引くことによって、未知のデータに対して推論することができるのでした。このアイデアについては最小二乗法（理論編）でも解説しました。
ここで引きたい線は次の図のような線ですね。
この関数（曲線）の式が分かれば推論ができるのですが、実験するたびにデータの分布がどんな関数で表せるか調べるのは手間ですし、もっと複雑でたくさんの要素が絡み合った結果得られるデータは、単純な関数で表せないことがほとんどです。
そのため、とても単純な関数を用意し、その重ね合わせで表現してやります。
この時用意した単純な形の関数を基底関数と呼び、基底関数を用いたフィッティング方法を基底関数モデルと呼びます。
正規関数
ここで関数の紹介です。
次回：線形基底関数モデル（実装編）でも基底関数として使うのが正規関数です。
機械学習の勉強をしていると必ず出会うことになる重要な関数の一つです。
この関数は確率密度を表す目的でも使うことができ、今後記事にするであろう教師なし学習ではこの正規関数（確率の世界では正規分布と呼ぶ。）なしでは話が進みません。
ちなみに、正規関数のことをガウス関数と呼ぶ場合もありますが、全く同じもののことを指します。
※wikipediaより引用
実際にプロットしてみると以下の通りです。
山のような形をしており、-∞から∞までの積分を行うと１になります。（つまり面積が１。確率を表すのに都合が良いことが分かります。）
正規関数は一般にN[µ,σ]と表します。
μは正規関数の平均値で、山形の頂点に当たる部分がどこにあるのかを表しています。上図だと頂点が横軸の０の地点にありますね。（μ＝０）
σは正規分布のなだらかさを表します。値が大きくなればなるほど平たい正規分布になっていきます。
今回のフィッティング法ではあまりにも複雑になってしまうのでσは固定し、μは任意の値を我々で指定していきます。
さて、この正規分布をいくつか用意し重ね合わせます。
下図は赤、青、緑の３つの正規関数を合成し、黄色の曲線を作った図です。
N[μ,σ] → N[-6,3] N[0,3] N[5,5]
このようにいくつかの正規関数を合成し、目的のデータ点とぴったり合うような曲線を作っていきます。
複数の正規関数を合成してフィッティングしていきます。そこでは各正規関数に重みをつけて合成しますが、この重みをパラメータとして学習していきます。
つまり、今回調べたい「それっぽい線」とは以下のように表せるはずです。
aN[μ,σ]+bN[μ,σ]+cN[μ,σ]+...
この正規関数にかけられている[a,b,c,...]を勾配降下法で推定していきます。
次回：線形基底関数モデル（実装編）では勾配降下法で分かりやすいように記述していたところを、行列式を使って一般化していきますので違いもお伝えしていきます。
勾配降下法（実装編）
今回の記事では勾配降下法の実装を行います。
理論の勉強をしたい方は前回：勾配降下法（理論編）を参考にしてください。
また、今回用いるデータセットは最小二乗法（実装編）で作成したものをそのまま使っていきます。
では、始めていきましょう。
勾配降下法の実装について
微分法の実装についての注意
勾配降下法を実装するうえで避けては通れないことは偏微分を実装することです。「偏微分が実装できる」ということは「微分が実装できる」ことにもなりますし、「極限を理解している」ということにもなります。
理論編でもそういった高校数学の知識は端折って説明していたので実装編でも簡単な説明のみで進めさせていただきます。
注意事項としては微分は極限をとります。我々人間はある値を０に近づけたり、無限に近づけたりすることで大体の値を得ます。
しかし、コンピュータは厳密に数値計算をしようとするため、桁が溢れて正しい数値が出なかったり、収束せずに解が得られなかったりします。
なので、「真の値に近似するように仕組みを作り、ある程度の妥協をする必要がある」ということを念頭において作業を進めてください。
１変数関数でプログラムを理解
1.)初期パラメータと関数の準備
仕組みを理解するためにまずは１変数関数を適当に用意。
初期パラメータも適当な値をxに代入します。
x=-100#適当な初期パラメータ def g(x): y=2*(x-3)**2#任意の関数（極値持っていればどんな関数でも良いと思います。） return y
今回の関数はx＝３で極値を持っています。
2.)微分法の実装
微分の定義は平均変化率の増加量を０に近づけることでした。
よってプログラム中のepsilonに対し、０に近い数値を与えることで近似します。
def DER(func,x): epsilon=0.001 mse=func(x+epsilon)-func(x-epsilon)#yの増加量 der=mse/(2*epsilon)#yの増加量/xの増加量 return der
3.)微分した値を使ってパラメータの更新
微分した数値はそのままパラメータの更新に用いると発散して解が求まらない可能性があるため、学習率（プログラム中のalpha）をかけます。
for i in range(1000):#学習回数1000回 alpha=0.1#学習率 der=DER(g,x) x=x-der*alpha print(x)
xの値が３に近づいているのが分かると思います。
良い感じに出来上がりましたね。
最小二乗法を勾配降下法で解く
さて、本題です。
最低でも２変数関数の偏微分が必要です。しかし、偏微分といっても特別なことをするわけではありません。
ある一つの変数以外を固定してパラメータの更新をかけてやれば良いだけです。
4.)初期パラメータ・データセットの準備
実際に見てみます↓
#2.2 勾配降下法 import matplotlib.pyplot as plt import numpy as np import math import random %matplotlib inline #初期パラメータ a=1000 b=1000 #身長から算出する理想体重 def BMI(x): y=x*x*22 return y #データの作成（x：身長、y：体重） x=np.array([]) y=np.array([]) for i in range(20): tall=np.random.normal(1.75 , 0.1) x=np.insert(x,len(x),tall,axis=0) norm=np.random.normal(0 , 2) weight=BMI(x[i])+norm y=np.insert(y,len(y),weight,axis=0) print(x) print(y)
※初期パラメータ以外は前回：最小二乗法（実装編）のコードから引用しただけです。
5.)偏微分の実装
今回はパラメータがaとbの２つあるので分かりやすいようあえて２つの関数に分けて定義しました。
本来であれば行列式を使い、関数も一つにまとめた方がスマートなのですが、分かりやすさを重視して実装していきたいと思います。
（「スマートではない」ということは３変数、４変数関数…など、「変数が増えるごとに実装が面倒なプログラムの書き方である」ということになります。）
１つ目はaの更新用に偏微分したものです。
def DER1(func,a,b,x,y): epsilon=0.0001 mse=func(a+epsilon,b,x,y)-func(a-epsilon,b,x,y) der=mse/(2*epsilon) return der
２つ目はbの更新用に偏微分したものです。
def DER2(func,a,b,x,y): epsilon=0.0001 mse=func(a,b+epsilon,x,y)-func(a,b-epsilon,x,y) der=mse/(2*epsilon) return der
3.)学習させる。
aの更新とbの更新は交互に行うことで勾配を降っていきます。。
for i in range(1000): alpha=0.5 der=DER1(min2error,a,b,x,y) a=a-der*alpha#aの更新 print(a) der=DER2(min2error,a,b,x,y) b=b-der*alpha#bの更新 print(b)
学習して得られたパラメータ（a,b）が、前回：最小二乗法（実装編）で得られた値と近い数値が得られればそれで学習成功となります。
print("学習結果") print(a) print(b)
（今回私が実行したデータセットは、前回：最小二乗法（実装編）で実行したデータセットと異なるため、少し違う数値が現れています。）
