default eye-catch image.

決定木回帰と決定木の作り方1

[mathjax] アンサンブル学習とかランダムフォレストに入門する前に決定木に入門する。 決定木はやっていることが直感的でわかりやすい。 決定木回帰と決定木分類。 ここよりはドメインとの連結部分が大変なんだろうと思った。 あと、Pythonは練習しないとな。 CART(Classification and Regression Tree)法(単に分類と回帰を英語にしただけだ!) 木を作るのだけれども、それが面白かったので今回と次回で書いてみる。 決定木回帰 (y=(x-0.5)^2)という2次曲線に従う事象があるとして、(y)を観測するとする。 観測による誤差が平均(mu=0)、分散(sigma^2=0.1)の正規分布に従うとして(y=(x-0.5)^2+epsilon)。 区間([0.0,1.0])の間に等間隔に存在する観測値(x,y)。 試行毎に(epsilon)が変わってくるので、毎回異なる。 この区間に(16)個のデータがあるとして、それを訓練データとして使ってモデルを作る。 (x_{train},y_{train})とする。 モデルの作成(学習)はfit()。 Scikit-learnに全て用意されていてデータを放り込むだけでモデルが出来る。 決定木の葉の最大値を(5)としている。他のパラメタは全部デフォルト値。 import numpy as np # 区間[0,1]上に16個の点を等間隔に生成する X_train = np.linspace(start=0,stop=1,num=16) y_train = (X_train - 0.5) ** 2 + np.random.normal(loc=0.0,scale=0.1,size=16) # 16x1配列を1x16に整形 X_train = X_train.reshape(16,1) print(X_train) # 決定木回帰 from sklearn.tree import DecisionTreeRegressor DTR = DecisionTreeRegressor(max_leaf_nodes=5) DTR.fit(X_train,y_train) 出来上がったモデルにテスト用データを流し込んでみる。 区間([0.0,1.0])に100個のデータを発生させてpredict()を呼ぶ。 最後、訓練データと回帰結果を同じグラフを書いてみて終了。 # 区間[0,1]上に100個の点を等間隔に生成する X_test = np.linspace(0,1,100) X_test = X_test.reshape(100,1) # 回帰! y_predict = DTR.predict(X_test) X_train = X_train.reshape(16) X_test = X_test.reshape(100) # 描画 import matplotlib.pyplot as plt plt.scatter(X_train,y_train) plt.plot(X_test,y_predict) plt.savefig(\'img.png\') plt.show() 葉の最大値を(5)としたので、木の深さが規定されて、 階段の個数が決まっている。 以下、6回分のモデルと予測の同時プロット。 (epsilon)の変化により訓練データが微妙に変わるだけで、 決定木の構造がむちゃくちゃ変化するのが特徴。 それっぽく言うとロバスト性が無いとか。 訓練データによって決定木がかなり違うことを利用して、 複数の決定木から多数決で結果を得ようというのがアンサンブル学習の試み。 むー。順々に詰めていこう。 max_leaf_nodesをデータの個数-1と一緒にすると、 全ての訓練データを通るモデルを作ることができる。 訓練データに対しては100%の精度が出るが、未学習のデータに対して答えられなくなる(ほんとに?)。 これが過学習(overfitting)。 Rによる 統計的学習入門posted with amazlet at 19.04.22Gareth James Daniela Witten Trevor Hastie Robert Tibshirani 朝倉書店 売り上げランキング: 118,792Amazon.co.jpで詳細を見る

default eye-catch image.

賃貸物件データをスクレイピングして散布図を出力してみた話

Hello CRISP-DM 1日目。 統計の基本は一応学びなおしたつもりなのと、 機は熟した感があり、Hotなところを突いてみる。 多変量解析と機械学習 多変量解析は、統計学の1分野で、人間が説明できるモデルを作ることが目的。 モデルができれば予測できるんじゃなかろうか、とも思うけども、 対比されちゃってるけど、統計学なんて凄まじい歴史がある訳で、 予測のための新しいツールとして機械学習があるという話。 モデルを説明することを放棄して、高性能な予測なり分類なりをしましょうと。 説明できないモデルを使って責任ある決定とか出来んのか、というのもあり。 完全に理解することを放棄するという位置付けでもない様子。 ルールベースと機械学習 予測器としての側面で、ルールベースとの比較がある。 密度関数が100%理解可能(予測値が出現する規則が100%理解可能)なのであれば、 いっそのこと全てif文で書けば良いわけで、わざわざ難しくする必要はなく..。 ルールベースの予測器を作ってみて結構上手く行った経験があって、 これはこれで何が起きても説明できるし、これで行けるならこれで良いかも。 Windows95の時代からある。 家賃予測してみる 家賃の予測モデル?を作るのを目的にCRISP-DMというプロセスを回す話。 ドメイン領域から目的変数と説明変数を取り出して、モデルを作って評価して、 正しくなさそうなら元に戻ってやり直す、という普通っぽいやり方にCRISP-DMという名前が付いてる。 データマイニングのプロセスで、特別に何か特別なツールを使わなければならない、ということはなし。 kaggleのTitanic..だと確かにリアル感がないのと、 結構、家賃だけでも知らないデータが埋もれてるものだな、と関心したり。 スクレイピングしてみる suumoから地元の賃貸物件情報を取得して表示してみた。 スクレイピングするためのコードはこちらを参考にしてみた。 同時接続数が1なので4200件くらいの取得に8時間くらいかかってしまうけど、 待てば良い話なので、寝る前にスタートして起きたら確認。 ログだとかセンサーだとかWebだとか、いろいろなデータソースがあるけれども、 取ってくるためにはかなり泥臭い感じになるし、他に方法がないのかも。 なるべくコードを書かずに済ませたい、けども、致し方ないのかも。 Octparseで取ろうとしたらWindows専用でMacのは無いと言われ..。 散布図を作る pandas(panel datas)というPython製のライブラリと、 seabornというPython製の描画ライブラリだけで簡単に散布図が出来る。 散布図出しただけで気分が良いな。 データが規則正しく並んでるのが気持ちが良いのか。 現実世界に存在する賃貸物件の情報が綺麗に2軸のグラフに収まってる。 もうこれだけで目的の半分くらいは達成したんじゃなかろうか..。 データの次数は5なので5x5の散布図ができる。 散布図作成時にある列に関して色分けできる機能があって、 試しに階数で色分けしてみた。微妙..。 suumoで期待した通りの書式になってなってないデータは0を入れてるので、 下限の0に張り付いてるデータが存在する。 1Kとか2LDKとか、間取りの区別はしていないから同時に出る。 家賃と占有面積は関係があるということがわかる様子。 広い方が高いよ、は言える様子。 他の因子に引きずられず、築浅であろうが駅近であろうが、広さは正義らしい。 家賃と築年数は、全体的に右肩下がりだけども、 占有面積ほど関係があるとは言えない様子。 現実的には、築浅だとしても駅から遠ければお安い、という話があると思い..。 全体的に右肩下がりなので、そうであっても古ければ概ね安いという。 もっと意外なのは、家賃と最寄駅からの時間。 駅近であることが一番影響してそうに思ってたけどもそんなことはなかった。 確かに、駅近であっても古くて狭いのは安いからな。 占有面積と築年数ほど、家賃に影響を与えないらしい。 思い込みというのは怖いもの..。 多変量解析なら、どれの係数が高いかはここで当たりがつくけれども、 ランダムフォレストで家賃予測するとすると、どうなるか、次回に続く。

default eye-catch image.

共役事前確率分布と逆ガンマ分布再考

[mathjax] 今回から、atomにmathjax-wrapperを入れてatomで数式入り文章を書いてみる。 数式を書くと瞬時にプレビューができて、10倍は速くなった気がする。1時間くらいで書いてみる。 母集団の(sigma)、(mu)が既知なのなら、共役事前分布は正規分布で良いのだから、 一体何をやりたいのか...という記事になってることに1日経って気づいた...。 チラシの裏程度なので気にしない。 何故、事前分布と事後分布を同じ確率分布に揃えるかというと、計算済みの事後分布を次の事前分布として 使うことができるから。事後分布の計算を繰り返していくと理論上は精度が上がっていく。 事前分布と事後分布を同じ分布にするために、逆ガンマ分布を事前分布として採用する。 そうすることで、左辺の事後分布と、右辺の尤度x共役事前分布が共に逆ガンマ分布となる。 (母平均が既知で分散が未知、という条件がつく。) 正規分布の共役事前確率分布として逆ガンマ分布を使う、というところまでは良かったのだけども、 逆ガンマ分布のパラメタが巷には2種類なのと4種類なのがあってもやもや。 ガンマ分布をさらに一般化した一般化ガンマ分布が4パラメータなのだそうな。 そもそもガンマ分布はカイ二乗分布と指数分布の一般化なので、どこまで一般化するのか。 ひとまず触れてはいけないところに触れてしまったようなので、2パラメータで出直してみる。 合ってるんだか間違ってるんだかも不明なのだけども、結論としては、細かいことはどうでもよくて、 事後分布と尤度x共役事前分布の関係を脳裏に焼き付けるためのプロセス。 何周かしないと真理にはたどり着けない...。 ガンマ分布と逆ガンマ分布 ガンマ分布の確率密度関数は以下。(alpha)は形状パラメータ,(beta)はスケールパラメータ。 (alpha) を固定して (beta)を動かすとカイ二乗分布。 (beta) を固定して (alpha)を動かすと指数分布。 Excelで (alpha) と (beta)をグリグリ動かしてみると意味がわかる。面白い。 良くできてるなー。 $$ begin{eqnarray} f(x|alpha,beta) = frac{1}{beta^alpha Gamma(alpha)} x^{alpha-1} expleft(-frac{x}{beta}right) end{eqnarray} $$ そして、(alpha) を大きくしていくと形状が正規分布に近づいていく。 累積確率分布関数は以下。 $$ begin{eqnarray} F(x|alpha,beta) = frac{1}{beta^alphaGamma(alpha)} int_0^xt^{alpha-1} exp (e^{-frac{t}{beta}})dt end{eqnarray} $$ 対して逆ガンマ分布の確率密度関数は以下。ガンマ分布で(x\'=frac{1}{x}) とすると出てくる。 $$ begin{eqnarray} f(x|alpha,beta)=frac{beta^{alpha}}{ Gamma(alpha)} x^{-(alpha+1)} expleft(-frac{beta}{x}right) end{eqnarray} $$ The shorthand for the distribution, X~inverted gamma(α,β), or IG(α, β), means that a random variable X has this distribution with positive parameters α and β. 再びベイズ統計へ ベイズの定理を使って出てくる事後確率分布と事前確率分布の関係は以下の通りだった。 $$ begin{eqnarray} p(mu,sigma^2|boldsymbol{x})&=&frac{p(boldsymbol{x}|mu,sigma^2)p(mu,sigma^2)}{p(boldsymbol{x})} \\ &propto& p(boldsymbol{x}|mu,sigma^2)p(mu,sigma^2) end{eqnarray} $$ 尤度は周辺確率の積として表されるけれども、母集団の平均が既知で分散が未知の場合の話をするので、 確率変数(X)を分散にもつ正規分布を考える。(y)は標本平均。(mu)は既知の母平均。 $$ begin{eqnarray} p(x | mu, X) &=& N(x|mu,X) \\ &=& frac{1}{sqrt{2pi x}}expleft(-frac{(y-mu)^2}{2x}right) end{eqnarray} $$ 正規分布である尤度に逆ガンマ分布をかけて式変形していく。 比例の関係を見たいので定数をのぞいてシンプルにするところがポイント。 $$ begin{eqnarray} p(mu,sigma^2|boldsymbol{x}) &propto& frac{1}{sqrt{2pi x}}expleft(-frac{(y-mu)^2}{2x}right) frac{beta_0^alpha}{Gamma(alpha_0)} left(frac{1}{x}right)^{alpha_0+1} expleft(frac{-beta_0}{x}right) \\ &propto & frac{1}{x} exp left( - frac{(y-mu)^2}{2x}right) x^{-alpha_0-2} exp left( -frac{beta_0}{x} right) \\ &propto & x^{-alpha_0-2} exp left( -frac{1}{x}(beta_0 + frac{1}{2}(y-mu)^2 right) end{eqnarray} $$ 右辺の確率分布が逆ガンマ分布であるためには、定数のぞいた箇所、つまり(x)の肩と(exp)の肩が 逆ガンマ分布のそれと同じでなければならないから、 $$ begin{eqnarray} -(alpha_0-2) &=& alpha +1 \\ alpha &=& alpha_0+frac{1}{2} \\ end{eqnarray} $$ また、 $$ begin{eqnarray} beta &=& beta_0 + frac{1}{2} (y-mu)^2 end{eqnarray} $$ これで、左辺の逆ガンマ分布の (alpha)、(beta)が作れる。 更新するたびに (alpha)、(beta)が増加していくのがわかる。

default eye-catch image.

MAP推定

[mathjax] 真の(mu)、(sigma)があるものとしてデータ(boldsymbol{x}=(x_1,x_2,cdots,x_n))が最も起こりやすくなるように (mu)、(sigma)を計算するのが最尤推定だった。 最尤推定においては、(mu_{ML})、(sigma_{ML})が先にあって、その結果として(boldsymbol{x})があると考えている。 逆に、データ(boldsymbol{x}=(x_1,x_2,cdots,x_n))があって、 その結果として(mu)、(sigma)が確率的な分布をもって存在すると考える。 (mu)、(sigma)はデータ(boldsymbol{x}=(x_1,x_2,cdots,x_n))次第で確率的に決まる。 この分布を事後確率分布と呼んでいて、 事後確率分布を最大にする(mu_{MAP})、(sigma_{MAP})を求めるのが最大事後確率推定(MAP推定)。 まずは、事後確率分布と事前確率分布の関係をまとめてみて、 最大化すべき事後確率分布とは何なのか、詳細に追ってみる。 最尤推定 (N)個のデータ(boldsymbol{x}=(x_1,x_2,cdots,x_n))が観測されたとして、 観測値にはガウスノイズが含まれる(=観測値は正規分布に従う)。 (boldsymbol{x})が互いに独立であるならば、同時確率は周辺確率、つまり各々の確率の積で表せた。 (boldsymbol{x})は既知のデータ、統計量である(mu)、(sigma)は変数であるから、 同時確率は(mu)、(sigma)を変数とする関数として表せた。 begin{eqnarray} p(x | mu, sigma^2) = prod_{n=1}^{N}N(x_n|mu,sigma^2) end{eqnarray} 既にある観測値(boldsymbol{x})を最も高い確率で実現するための統計量(mu)、(sigma)が 一つ存在するという仮定のもと(boldsymbol{x})を推定した。 (p(x|mu,sigma^2))は、(mu)、(sigma)がとある値であるときに、 (boldsymbol{x})が発生する確率として捉えることができる(条件付き確率)。 事後確率分布 逆に、(boldsymbol{x})が得られたときに(mu,sigma^2)が得られる確率(p(mu,sigma^2|x))を考える。 (boldsymbol{x})自体が確率的に生起するデータであり、(mu,sigma^2)が(boldsymbol{x})次第でブレる。 そんな状況。 ベイズの定理を使って因果を入れ替えると以下の通り。 begin{eqnarray} p(mu, sigma^2 | x) = frac{p(boldsymbol{x|}mu,sigma^2)p(mu,sigma^2)}{p(boldsymbol{x})} end{eqnarray} 分子の(p(boldsymbol{x|}mu,sigma^2))は尤度。(p(mu,sigma^2))は事前分布。 データ(boldsymbol{x})が事前分布(p(mu,sigma^2))に従うと仮定したもので、 最初からそれがわかっていれば苦労しない訳だから、(p(mu,sigma^2))はわからない。 (p(boldsymbol{x}))はデータ(boldsymbol{x})が発生する確率。 観測することから得られる確率。(mu,sigma^2)に対して定数。 この(p(mu, sigma^2 | x))を最大化する(mu, sigma^2)を求めようというのが、 最大事後確率推定。 (p(mu, sigma^2 | x))の式において、分母の(p(boldsymbol{x}))は(mu, sigma^2)に対して定数だから、 (p(mu, sigma^2 | x))を最大化する(mu, sigma^2)は、分母を払った(p(boldsymbol{x|}mu,sigma^2)p(mu,sigma^2))も最大化する。 だから、代わりに(p(boldsymbol{x|}mu,sigma^2)p(mu,sigma^2))の最大化を考える。 begin{eqnarray} p(mu, sigma^2) p(x | mu, sigma^2) &=& N(x_n|mu,sigma^2) prod_{n=1}^{N}N(x_n|mu,sigma^2) \\ &=& N(x_n|mu,sigma^2) p(x | mu_{ML}, sigma^2_{ML}) end{eqnarray} 尤度関数(p(boldsymbol{x|}mu,sigma^2))の最尤解(mu_{ML})と(sigma^2_{ML})の求め方は前回記事。 [clink url=\"https://ikuty.com/2019/01/23/maximum_likelihood_estimation-2/\"] データさえあれば標本平均、標本分散を計算して終了。 begin{eqnarray} mu_{ML} &=& frac{1}{N}sum_{i=1}^{N}x_i \\ sigma^2_{ML} &=& frac{1}{N} sum_{i=1}^N (x_i -mu_{ML})^2 end{eqnarray} 尤度(prod_{n=1}^{N}N(x_n|mu,sigma^2))は計算により値が決まる。 begin{eqnarray} prod_{n=1}^{N}N(x_n|mu,sigma^2) &=& prod_{j=1}^{n} frac{1}{sqrt{2pi} sigma_{ML}} exp left( - frac{1}{2} left( frac{x_j -mu_{ML}}{sigma_{ML}} right)^2 right) end{eqnarray} 共役事前確率分布 さて... 事後確率と事前確率の関係をわかりやすく書くと以下の通り。 だが、事前確率(p(mu,sigma^2|x))は分からない。 begin{eqnarray} p(mu,sigma^2|x) = frac{p(x|mu,sigma^2)}{p(x)} p(mu,sigma^2) end{eqnarray} 事後確率分布(p(mu,sigma^2|x))が事前確率分布(p(mu,sigma^2))と同じ確率分布であれば、 事前確率分布と尤度の計算から事後確率分布を求められる。 そういう操作をするために、事後確率分布と同じ形の事前確率分布を仮に作っておいて、 まず事後確率分布と仮の事前確率分布を結ぶ式を立てる。 仮に立てた事前確率分布を、共役事前確率分布と言うらしい。 いきなり事前確率分布を(mu)、(sigma)だけをもつ別の確率分布に置き換えることはできないから、 確率分布を制御する超パラメータ(alpha,beta,gamma,cdots)を使って表現する。 正規分布の共役事前分布は逆ガンマ分布を使う。 逆ガンマ分布は4個の超パラメータ(alpha,beta,gamma,delta)を持つ。 begin{eqnarray} p(mu,sigma^2) = frac{sqrt{gamma}}{sigma sqrt{2pi}} frac{beta^alpha}{Gamma[alpha]} left( frac{1}{sigma^2} right)^{alpha+1} exp left( - frac{2beta + gamma (delta - mu)^2 }{2sigma^2} right) end{eqnarray} なので、ベイズの式は超パラメータ(alpha,beta,gamma,delta)を使って以下のように書ける。 左辺の事後確率分布と右辺の共役事前確率分布は同じ分布であって、 右辺には定数(1/p(x))と尤度( prod_{n=1}^{N}N(x_n|mu,sigma^2))をかけた形になっている。 begin{eqnarray} p(mu,sigma^2|x) &=& frac{p(x|mu,sigma^2)}{p(x)} p(mu,sigma^2) \\ &=& frac{1}{p(x)} prod_{n=1}^{N}N(x_n|mu,sigma^2) frac{sqrt{gamma}}{sigma sqrt{2pi}} frac{beta^alpha}{Gamma[alpha]} left( frac{1}{sigma^2} right)^{alpha+1} exp left( - frac{2beta + gamma (delta - mu)^2 }{2sigma^2} right) end{eqnarray} 左辺の事後確率を最大化する(mu,sigma)を見つけるわけだけども、 それは、尤度と共役事前確率で表された右辺を最大化する(mu,sigma)を見つけるのと同じ。 対数をとると積の部分が和になって計算しやすいから対数をとる。 やり方は(mu)、(sigma)で偏微分したものが0とする。 begin{eqnarray} frac{partial}{partial mu} p(mu,sigma^2|x) &=& frac{partial}{partial mu} frac{p(x|mu,sigma^2)}{p(x)} p(mu,sigma^2) \\ &=& frac{partial}{partial mu} frac{1}{p(x)} prod_{n=1}^{N}N(x_n|mu,sigma^2) frac{sqrt{gamma}}{sigma sqrt{2pi}} frac{beta^alpha}{Gamma[alpha]} left( frac{1}{sigma^2} right)^{alpha+1} exp left( - frac{2beta + gamma (delta - mu)^2 }{2sigma^2} right) \\ &=& 0 end{eqnarray} ちょっと計算しようとすると目眩がするけど、超パラメータ(alpha,beta,gamma,delta)はそのまま残して、 (mu,sigma)を(alpha,beta,gamma,delta)を含む式で表す。 begin{eqnarray} mu_{MAP} &=& frac{sum_{i=1}^{N}x_i + gamma delta}{N+gamma} \\ &=& frac{Nbar{x} + gamma delta}{N+gamma} end{eqnarray} 同様に、 begin{eqnarray} frac{partial}{partial sigma^2} p(mu,sigma^2|x) &=& frac{partial}{partial sigma^2} frac{p(x|mu,sigma^2)}{p(x)} p(mu,sigma^2) \\ &=& frac{partial}{partial sigma^2} frac{1}{p(x)} prod_{n=1}^{N}N(x_n|mu,sigma^2) frac{sqrt{gamma}}{sigma sqrt{2pi}} frac{beta^alpha}{Gamma[alpha]} left( frac{1}{sigma^2} right)^{alpha+1} exp left( - frac{2beta + gamma (delta - mu)^2 }{2sigma^2} right) \\ &=& 0 end{eqnarray} 計算できないけど begin{eqnarray} sigma^2_{MAP} &=& frac{sum_{i=1}^{N} (x_i-mu)^2 + 2beta + gamma (delta-mu)^2}{N+3+2alpha} \\ &=& frac{Nsigma^2 + 2beta + gamma (delta-mu)^2}{N+3+2alpha} end{eqnarray} 最尤推定による解と並べてみる。 平均、分散ともに、最尤解の分母、分子をパラメータ(alpha,beta,gamma,delta)で調整したものになっている。 begin{eqnarray} mu_{ML} &=& bar{x} \\ mu_{MAP} &=& frac{Nbar{x}+gamma delta}{N+delta} \\ sigma_{ML}^2 &=& bar{sigma^2} \\ sigma_{MAP}^2 &=& frac{Nsigma^2 + 2beta + gamma (delta-mu)^2}{N+3+2alpha} end{eqnarray} パターン認識と機械学習 上posted with amazlet at 19.01.22C.M. ビショップ 丸善出版 売り上げランキング: 21,190Amazon.co.jpで詳細を見る 統計学入門 (基礎統計学Ⅰ)posted with amazlet at 19.01.22東京大学出版会 売り上げランキング: 3,133Amazon.co.jpで詳細を見る [youtube id=QnqlP4O0TXs]

default eye-catch image.

最尤推定

[mathjax] これは単なる趣味の記事です。 たいして数学が得意でもなかった偽理系出身のおっさんの懐古趣味。 PRMLは省略が多くて難しいです。 省略を紐解いても出てくるのが統計学入門レベルです。 鬼ですね。 ただ、厳密さを求めない簡単な説明をしてくれているところを先に読むと、 割と簡単に紐解けます。 \"互いに独立であること\"の復習 少し前に、\"違いに独立\"の定義を詳しめに読んだことがあった。 2つの確率変数(X)、(Y)があったとき、 (X=x)であり同時に(Y=y)である確率(P(X=x,Y=y)=f(x,y))について、 確率を定義できる空間(2次元ユークリッド空間(S)...)において、 (X,Y)が同時に起こる特別な事象(A)((S)の部分集合)を定義した。 begin{eqnarray} iint_S f(x,y)dydx = 1, f(x,y) geq 0 \\ P((X,Y) in A) = iint_A f(x,y)dydx end{eqnarray} 2変数の片方について積分した分布を周辺確率分布(marginal probability distribution)とよんだ。 begin{eqnarray} g(x) &=& sum_y f(x,y) \\ h(x) &=& sum_x f(x,y) end{eqnarray} (g(x),h(x))は、(f(x,y))において、片方の変数を固定して足し算したもの。 では、(g(x),h(x))から(f(x,y))を作れるかというと、 (g(x),h(x))が\"違いに独立\"である場合に限り、(g(x)とh(x))の積と(f(x,y))が等しくなった。 [clink url=\"https://ikuty.com/2018/10/17/statistical-independence/\"] 最尤推定 (N)個のデータ(boldsymbol{x}=begin{pmatrix}x_1 \\ x_2 \\ cdots \\ x_n end{pmatrix})が、平均(mu)、標準偏差(sigma)の正規分布(N(mu,sigma))から独立に生成されるとする。 また、(boldsymbol{x})は互いに独立であるとする。 各々のデータが発生する確率(P(x_i|mu,sigma)=Nleft(x_n | mu, sigma^2right))を 同時確率(P(boldsymbol{x}|mu,sigma))の周辺確率と考えることができる。 つまり以下の通り。 begin{eqnarray} p(boldsymbol{x} | mu, sigma^2 ) = prod_{n=1}^N Nleft(x_n | mu, sigma^2right) end{eqnarray} ここで、(boldsymbol{x})は既に与えられている訳なので、 (p(boldsymbol{x} | mu, sigma^2 ) )の変数部分は(mu)、(sigma)となる。 つまり、(p(boldsymbol{x} | mu, sigma^2 ) )は(mu)、(sigma)の関数。 (p(boldsymbol{x} | mu, sigma^2 ))を最大化する(mu)、(sigma)を求めようというのが、 最尤推定の考え方。 (p(boldsymbol{x} | mu, sigma^2 ) )の最大値を与える(mu)、(sigma)は、 両辺に対数をとった後でも最大値を与えるから、 計算しやすくするために両辺に対数をとる。 対数をとることで、右辺の積が対数の和になる。 高校数学で暗記したやつの連発。 begin{eqnarray} log p(boldsymbol{x} | mu, sigma^2 ) &=& log prod_{n=1}^N Nleft(x_n | mu, sigma^2right) \\ &=& -frac{1}{2sigma^2} sum_{n=1}^N left( x_n -mu right)^2 - frac{N}{2} log sigma^2 - frac{N}{2} log left( 2pi right) end{eqnarray} (mu)に関して最大化するためには、 (sigma)を定数として固定して(mu)で偏微分して、それが0とする。 右辺は第2項,第3項がゼロとなり、第1項の微分が残る。 begin{eqnarray} frac{partial }{partial mu} log p(boldsymbol{x} | mu, sigma^2 ) &=& frac{1}{sigma^2} sum_{n=1}^{N} left(x_n-muright) \\ &=& 0 end{eqnarray} これ、すなわち以下。 begin{eqnarray} frac{1}{sigma^2} left( sum_{n=1}^{N} x_n - Nmu right) &=& 0 \\ sum_{n=1}^{N} x_n &=& Nmu \\ mu = frac{1}{N} sum_{n=1}^{N} x_n end{eqnarray} なんと、(mu)は標本平均のときに最大化する!。 (sigma)に関して最大化するには、 (mu)を定数として固定して(sigma)で偏微分して、それが0とする。 これも高校数学で暗記したやつを連射する。 begin{eqnarray} frac{partial }{partial sigma} log p(boldsymbol{x} | mu, sigma^2 ) &=& left( -frac{1}{2sigma^2} sum_{n=1}^N left( x_n -mu right)^2 right)\' - left( frac{N}{2} log sigma^2 right)\' \\ &=& frac{1}{sigma^3} sum_{n=1}^N left( x_n -mu right)^2 - frac{N}{sigma} \\ &=& 0 end{eqnarray} これ、すなわち以下。 begin{eqnarray} frac{N}{sigma} = frac{1}{sigma^3} sum_{n=1}^N left( x_n -mu right)^2 \\ sigma^2 = frac{1}{N} sum_{n=1}^N left( x_n -mu right)^2 \\ end{eqnarray} なんと、(sigma^2)は標本分散のときに最大化する。 なぜ、(mu)、(sigma)それぞれ独立に偏微分して求めた値を使って、 もう片方を計算して良いのか、については時間がないから省略。 標本分散は不偏分散ではないことに起因して、 最尤推定により求められた分散は、不偏分散から過小に評価される。 [clink url=\"https://ikuty.com/2018/10/27/unbiased_variance/\"] PRMLでは、多項式曲線のフィッティングに最尤推定を使う内容になっていて、 最尤解が不偏でないことがわかりやすくなっている。 次回は、最大事後確率推定(maximum posterior)について。 以下を最大化する。 begin{eqnarray} pi(mu) L(mu) = frac{1}{sqrt{2pi}sigma_m} exp left( -frac{1}{2} left(frac{mu}{sigma_m}^2 right)right) prod_{j=1}^{n} frac{1}{sqrt{2pi} sigma_v} exp left( - frac{1}{2} left( frac{x_j -mu}{sigma_v} right)^2 right) end{eqnarray} パターン認識と機械学習 上posted with amazlet at 19.01.22C.M. ビショップ 丸善出版 売り上げランキング: 21,190Amazon.co.jpで詳細を見る 統計学入門 (基礎統計学Ⅰ)posted with amazlet at 19.01.22東京大学出版会 売り上げランキング: 3,133Amazon.co.jpで詳細を見る

default eye-catch image.

コサイン類似度 Cosine Simirality

[mathjax] BoW化した文章同士の類似度を求める一番メジャーなやり方。 内積の定義を式変形しただけ。 $$ begin{eqnarray} cos(vec{q},vec{d}) &=& frac{vec{q}cdot vec{d}}{|vec{q}| cdot |vec{d}|} \\ &=& frac{vec{q}}{|vec{q}|} cdot frac{vec{d}}{|vec{d}|} \\ &=& frac{sum_{i=1}^{|V|}q_i d_i}{sqrt{sum_{i=1}^{|V|}q_i^2} cdot sqrt{sum_{i=1}^{|V|}d_i^2}} end{eqnarray} $$ BoWが正規化された単位ベクトルの場合は分母が1になるから、 $$ cos(vec{q},vec{d}) = vec{q} cdot vec{d} = sum_{i=1}^{|V|} q_i d_i $$ 次元が3以下の場合はベクトルがどれだけ同じ向きを向いているか、という説明で十分だけど、4以上の場合は結局「類似」のような言葉になってしまう。 範囲 シュワルツの不等式というのがあって、 $$ sum_{i=1}^{n} p_i^2 cdot sum_{i=1}^{n} d_i^2 geq sum_{i=1}^{n} p_i d_i $$ コサイン類似度の分子は分母以下であることが決まる。 シュワルツの不等式の等号が成り立つとき、 コサイン類似度の分子と分母が同じとなり1となる。 pi, di のいずれかが 0のとき、qidiは0となる。 全部0となるとき、コサイン類似度の分子は0となる。 だから、定義の上ではコサイン類似度は0以上1以下だと思うんだ。 類似度? BoW化してタームの頻度を要素とした時点で、ベクトルは意味を持たないことになるから、ベクトル同士のコサイン類似度を取ったところでベクトル同士が似た意味を持つとは言えない。単にタームの出現頻度が似ているというだけ。 BoWではなく、特徴量を要素とするベクトル同士のコサイン類似度をとったのであればベクトル同士の特徴量が似ていることを表す。 文章をなんとか直交する特徴量からなるベクトルで表現することができれば、意味の類似度を求めることができる。 tf-idfに触れずに書くとこんな感じ。 文章が短いと 文章が短いと、文章がだいぶ異なるのに出現頻度は近い、という状況が起きやすくなり、あてにならなくなる。

default eye-catch image.

ロジスティック回帰とROC曲線

[mathjax] ロジスティック回帰 2値の特徴量tを持つトレーニングデータ(x,y,t)を正解/不正解に分類する問題を考えたとき、最尤推定によりモデルf(x,y)を決めるのがロジスティック回帰。モデルf(x,y)がトレーニングデータ(x,y,t)を正解と分類する確率P(x,y)を立てて、P(x,y)の最大化問題を解くことでf(x,y)のパラメタを決定する。その際、P(x,y)はロジスティック関数を使って以下の通り表される。 $$ begin{eqnarray} z_n &=& sigma(w^Tphi_n) \\ P(x,y) &=& z_n (1-z_n) end{eqnarray} $$ P(x,y)において、トレーニングデータを定数とすることでP(x,y)をパラメタwに関する関数とした。P(x,y)を最大化するwを求めたいのだが、P(x,y)を解析的に最大化することはできず、パラメタwを数値計算(反復再重み付最小二乗法、ニュートンラフソン法)することによりf(x,y)のパラメタwを更新していく。wが一定値に収束した場合、f(x,y)=0はトレーニングデータ内の正解データ(t=1)を1/2の確率で正解(t=1)として分類できることを示す。 正解分類確率の制御 トレーニングデータ(x,y,t)に対して、f(x,y)t>0の場合は「正解」、f(x,y)t0の方向あるt=0のデータ、f(x,y)0の方向に距離α分平行移動していくと、f(x,y)>0のときは正解であったもののf(x,y)>αのときは不正解となるデータが発生する。 逆にf(x,y)<0の方向に距離β分平行移動していくと、f(x,y)<0のときは正解であったもののf(x,y)0方向にα移動することで、トレーニングデータを「正解」として分類する確率を50%から75%に変更する等の操作にあたる。 ロジスティック回帰とROC曲線 分類問題を整理すると、分類対象とするトレーニングデータの特徴量tを2値のいずれであるか識別することであり、分類とその結果の組み合わせとして以下の組み合わせが考えられる。それぞれ統計の分野で名前が付いている。 t=1であるデータを1として分類する TP, True Positive t=1であるデータを0として分類する FN, False Negative t=0であるデータを0として分類する TN, True Negative t=0であるデータを1として分類する FP, False Positive かなり混乱するので図を張り付けておく。観測データが1なのか0なのかは\"Positive\"/\"Negative\"、それを1と分類したのか0と分類したのかは\"True\"/\"False\"。t=1/0は符号に過ぎないので\"Positive\"/\"Negative\"はその時次第。True/Falseは分類結果が特徴量と一致するかどうかを表す。 ロジスティック関数P(x,y)=σ(w0+w1x+w2y)は、トレーニングデータ(x,y)を1として分類する確率を返す関数である。 最尤推定によりwを求めたのだから、f(x,y)=w0+w1x+w2yは、確率的に最もトレーニングデータを正しく分類する多項式となっている。 \"確率的に最も\"の度合いを定量的に表すのがROC曲線とAUC。 全トレーニングデータ(x,y)に対してP(x,y)を計算し降順ソートすると、P(x,y)が高ければ(つまり上位であれば)t=1のデータが多いはずで、P(x,y)が低ければ(つまり下位であれば)t=0のデータが多いはずだ。P(x,y)に相当するTP,FPが期待される。実際には上位にt=0のデータ、下位にt=1のデータが観察されることがあり、単純にP(x,y)からTP,FPが決まるものではない。 TPがP(x,y)よりも相対的に高く、(1-P(x,y))よりも相対的に低ければ、多項式の性能が良いということになる。 定数c=P(x,y)を定めたときP(x,y)<c<1の範囲でt=1のものは True Positive, P(x,y)<c<1の範囲でt=0のものは False Positive。0<c<P(x,y)の範囲でt=0のものは False Negative, t=1のものは False Positive。cを1から0に下げていくと True Positive, False Positiveの数がそれぞれ独立に増加していく。各々のcを真陽性率、偽陽性率という。 横軸に偽陽性率FP、縦軸に真陽性率TPをとり、トレーニングデータ(x,y)をプロットした図がROC(Receiver Operating Characteristic)曲線。 P(x,y)を1から0に減らしていったときTP/FPそれぞれ増加するが、ROC曲線によりTPが大きくFPが小さいプロットが性能の良い分類であることを視覚的に表現できる。 ROC曲線の左下と右上を結ぶ直線は、0<P(x,y)<1の全域に渡りFP=TPであったことを示し、ランダムにTrueかFalseかを答えたときの視覚的表現にあたる。 左上に凸となっているROC曲線が性能の良さを表し、ROC曲線の下の面積AUC(Area Under the Curve)が分類器の性能を表す数値。

default eye-catch image.

ロジスティック回帰 教師有り確率的分類モデル

[mathjax] 2種類の値tを持つトレーニングデータ{(xn,yn,tn)}を正解/不正解に分類するパラメトリックモデルを立てて、パラメトリックモデルと不正解に分類されるトレーニングデータの距離E(w)が最小になるwを決定しモデルを作成するのがパーセプトロンであった。 E(w)においてトレーニングデータ(xn,tt)を定数として使用することでE(w)をwから決まる関数とし、wの決定をE(w)の最小化問題に帰着させる考え方は最小二乗法と同様であった。パーセプトロンは最小二乗法と異なり、E(w)の最小化問題を解析的に解けないため、数値計算によりE(w)を小さくする方向にwを更新する手順を取った。 最小二乗法はトレーニングデータから決定した多項式が未知のサンプルデータにおいても二乗誤差が最小になることを期待していたのと同様に、パーセプトロンも誤分類する未知のサンプルとの距離が最小になることを期待した分類モデルであることを意味する。最小二乗法もパーセプトロンも、サンプルデータのバラツキを考慮していない。 サンプルデータがモデルからσのバラつきをもって存在することを加味することで、パラメタwとバラつきσから決まる尤度関数を導入し、wの決定を尤度関数の最大化問題に帰着する最尤推定法に拡張した。最尤推定法により、パラメトリックモデルが変数とσの関数となり、確率的な回帰を検討できるようになった。 サンプルデータがモデルにより正解であると分類される確率を考慮することで、回帰モデルとして未来を予測する要素を増やすことができる。そこで、本エントリにおいてその手法(ロジスティック回帰)をまとめる。 ITエンジニアのための機械学習理論入門posted with amazlet at 17.03.10中井 悦司 技術評論社 売り上げランキング: 8,130Amazon.co.jpで詳細を見る ロジスティック回帰 モデルf(x,y)を仮定したとき、トレーニングデータ(x,y)について|f(x,y)|はモデルとの距離を表す。また、f(x,y)>0の方向にt=1のデータ、f(x,y)<0の方向にt=-1のデータが存在するモデルであることを示す。 f(xn,yn)tn > 0 であれば分類は正解であるから、f(x,y)が大きくなるにつれてt=1であるデータの発生確率が上がり、f(x,y)が小さくなるにつれてt=1であるデータの発生確率が下がることになる。逆に、f(x,y)が大きくなるにつれてt=-1であるデータの発生確率が下がり、小さくなるにつれてt=-1であるデータの発生確率が上がる。 横軸にf(x,y)、縦軸にt=1であるトレーニングデータがf(x,y)から得られる確率を取ると以下のような曲線が得られる。一般にロジスティック曲線と呼ばれる。 ロジスティック関数は以下のように表現される。 $$ sigma(a)=frac{1}{1+e^{-a}} $$ トレーニングデータ(x,y)で得られたデータの特徴量tが1である確率は以下の式で表される。 $$ P(x,y) = sigma(w_0 + w_1 x + w_2y) $$ また、t=1の場合とt=-1の場合で確率が逆になるため $$ begin{eqnarray} t_n=1の場合 &:& P(x_n,y_n) \\ P_n &=& P(x_n,y_n)^1{1-P(x_n,y_n)}^0 \\ &=& P(x_n,y_n) \\ t_n=-1の場合 &:& 1-P(x_n,y_n) \\ P_n &=& P(x_n,y_n)^0{1-P(x_n,y_n)}^1 \\ &=& 1-P(x_n,y_n) \\ end{eqnarray} $$ まとめて書くと以下となるが、これは全てのトレーニングデータがモデルから得られる確率であり尤度関数である。 $$ P_n = P(x_n,y_n)^{t_n}{ 1-P(x_n,y_n)}^{1-t_n} $$ ここで、 $$ P_n = sigma(w_0+w_1x+w_2y) $$ だから、 $$ P_n = sigma(w_0+w_1x_n+w_2y_n)^t{t_n}{1-sigma(w_0+w_1x_n+w_2y_n)}^{1-t_n} $$ 多項式のベクトル表現をロジスティック関数に入れた結果をznとする。 $$ begin{eqnarray} w &=& left ( begin{array}{c} w_0 \\ w_1 \\ w_2 end{array} right ) \\ phi_n &=& left ( begin{array}{c} 1 \\ x_n \\ y_n end{array} right ) \\ z_n &=& sigma (w^Tphi_n) end{eqnarray} $$ 尤度関数Pnをznを使って表すと以下となる。 $$ P_n = z_n^{t_n}(1-z_n)^{1-t_n} $$ 尤度関数を最大とするwを決定する最尤推定法によりwを決定する方法をロジスティック回帰という。 ロジスティック回帰における尤度関数の最大化 パーセプトロンにおいて数値計算(確率的勾配降下法)によりwの修正値を求めwを更新していったのと同様に、ロジスティック回帰の尤度関数最大化問題も解析的に解くことができず、数値計算によりwを修正していく。ニュートンラフソン法を使った反復再重み付最小二乗法により修正値を求めていく。導出方法は別途記述。wの修正式は以下の通りである。 $$ w_{new} = w_{old}-(Phi^TRPhi)^{-1}Phi^T(z-t) $$ なお、各行列は以下の通りである。Rはzn(1-zn)を対角成分とする対角行列。 $$ begin{eqnarray} t &=& left ( begin{array}{c} t_1 \\ t_2 \\ vdots \\ t_n end{array} right ) \\ Phi &=& left ( begin{array}{c c c} 1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ vdots & vdots & vdots \\ 1 & x_n & y_n \\ end{array} right ) \\ z_n &=& left ( begin{array}{c} z_1 \\ z_2 \\ vdots \\ z_n \\ end{array} right) \\ end{eqnarray} $$ 収束条件は以下の通り。ニュートン法の収束条件ぽい。 $$ frac{||w_{new}-w{old}||^2}{||w_{old}||^2} < 0.001 $$ なお、wの更新が進み、ロジスティック関数znが0または1に近づくか、0または1そのものになる場合がある。 $$ z_n = sigma (w^T phi_n) $$ zn=0または1となったとき、Rの対角成分zn(1-zn)が0となる。すると $$ Phi^TRPhi $$ の逆行列が存在しなくなりwを更新できなくなる。数値計算上の精度の話でもあるし、もともと全てのトレーニングデータが正解/不正解に分類できるのであれば、f(x,y)を大きくしていったときにf(x,y)と誤分類データとの距離が0になる、すなわちモデルから正解データを100%の確率で取得できる時がくるので、収束停止条件として考慮すべき内容である。

default eye-catch image.

パーセプトロン

[mathjax] (x,y)平面上のサンプルデータ(xn,yn)が、(x,y)平面を分割する1次多項式f(x,y)=w0+w1x+w2y を基準に f(xn,yn)>0 であれば正解、f(xn,yn)<0であれば不正解 のいずれかに分類する問題を考える。 正解、不正解を表す変数tについて、正解の場合に+1、不正解の場合に-1としたとき、トレーニングデータ(xn,yn,tn)を与えて多項式の係数wを求める分類アルゴリズムをパーセプトロンという。 ITエンジニアのための機械学習理論入門posted with amazlet at 17.03.10中井 悦司 技術評論社 売り上げランキング: 8,130Amazon.co.jpで詳細を見る パラメトリックモデルの決定 (x,y)平面を分割する1次多項式f(x,y)を以下の通り定義する。 $$ f(x,y)=w_0+w_1x+w_2y $$ トレーニングデータ(xn,yn)に対してパラメタtを以下の通り決定する。 $$ begin{eqnarray} f(x_n,y_n) > 0 Rightarrow t=+1 \\ f(x_n,y_n) < 0 Rightarrow t=-1 end{eqnarray} $$ tの符号とf(x,y)の符号の関係性から以下が言える。正解/不正解それぞれで場合分けして考えなくても良い。 f(xn,yn)>0なのにt=-1であるデータや、 0 Rightarrow 正解 \\ f(x_n,y_n)t < 0 Rightarrow 不正解 $$ モデルの評価式 トレーニングデータ(xn,yn,tn)を与えたとき、最良のパラメータw=(w0,w1,w2)を決めるため、多項式の評価式を決定する。正解、不正解を分類する問題であるから、全ての不正解データをf(x,y)に与えたときにその評価値が最小になるようにする。 これは、不正解データを無くすか、なるべく不正解データと距離が近いところを通るf(x,y)を見つけることに他ならない。ということで不正解データとの距離を定式化する。 $$ begin{eqnarray} E_n = | f(x_n,y_n) | \\ end{eqnarray} $$ 不正解データについてのみEnを加算する。ただし (xn,yn)は不正解データのみ $$ begin{eqnarray} E &=& sum_n E_n = sum_n |f(x_n,y_n)| \\ &=& -sum_n (w_0 + w_1x_n+w_2y_n)t_n end{eqnarray} $$ ここで、無理くりベクトル(行列))計算に持ち込む。 $$ begin{eqnarray} w &=& left ( begin{array}{c} w_0 \\ w_1 \\ w_2 end{array} right) \\ phi &=& left ( begin{array}{c} 1 \\ x_n \\ y_n end{array} right) \\ E &=& -sum_n t_nw^Tphi_n end{eqnarray} $$ モデルの評価式Eが決まった。Eはwにより決まるためE(w)とも書ける。トレーニングデータからE(w)を最小化するパラメタを求めていく。 確率的勾配降下法 E(w)をwについて偏微分するとE(w)の勾配ベクトルを求められる。(あぁ線形代数...。E(w)の勾配ベクトル∇E(w)とは、ベクトルwにおいてE(w)を最大化する向きと大きさを持つベクトルのことです...) $$ begin{eqnarray} nabla E(w) &=& left ( begin{array}{c} frac{partial E}{w_0} \\ frac{partial E}{w_1} \\ frac{partial E}{w_2} end{array} right ) \\ end{eqnarray} $$ 一方、上のベクトルを微分すると以下の式となる。 $$ nabla E(w) = - sum_{n}t_nphi_n $$ E(w)の最小化を考えると以下が成り立つはず。 $$ nabla E(w) = - sum_{n}t_nphi_n = 0 $$ ∇E(w)はtとxを含むがwを含まず、∇E(w)を式変形してwの式に展開することができない。wの式に展開できないと∇E(w)を0にするwを求める式を立てられない。 そこで∇E(w)の幾何学的な特徴を利用した数値計算によりE(w)を更新する。勾配ベクトル∇E(w)はwにおいてE(w)を最も大きくするベクトルだから、wと-∇E(w)を加算することでE(w)を小さくできる。 これを繰り返し実行することでE(w)を最小化できる。 実際のトレーニングデータ数Nは大きく、∇E(w)の行列計算は大変な作業となる。そこで、トレーニングデータから無作為にデータを1個取り出し、その1個分のデータについてのみwを更新する。N=1のとき∇E(w)は以下の通りとなる。 $$ - nabla E(w) = sum_n t_n phi_n = t_n phi_n $$ そのため、wの修正は以下の通りとなる。 $$ begin{eqnarray} w\' &=& w - ( -nabla E(w) ) \\ &=& w + nabla (Ew) \\ &=& w + t_n phi_n end{eqnarray} $$ E(w)が最小になるwが決まったとき、そのwを採用した多項式は、トレーニングデータについて誤りが少ない分類を行う式になっている。 トレーニングデータから無作為にサンプルを選んで勾配ベクトルの逆ベクトルを加算して降下していくから「確率的勾配降下法」という。へぇ。 パーセプトロンの収束速度 数値計算により値を更新していくアルゴリズムは、その収束速度が話題に上がる。パーセプトロンの確率的勾配降下法においてはトレーニングデータφのcにより調整する。 $$ begin{eqnarray} phi &=& left ( begin{array}{c} c \\ x_n \\ y_n end{array} right) end{eqnarray} $$ 確率的勾配降下法によりwを更新するときw0は1回の繰り返しで±cの範囲で増減する。cを大きくするとw0の増減範囲が大きくなり収束速度が速くなる。ではこのパラメータcはどのように決めるのか..はまた別のエントリで..。 $$ begin{eqnarray} w\' = w + t_n phi_n end{eqnarray} $$ 分類アルゴリズムとして、パラメトリックモデルが線形でない場合や、そもそも教師無しの場合など、他にもいろいろあるので、別エントリに書いていく。

default eye-catch image.

最尤推定 確率分布をもったトレーニングデータの学習

[mathjax] 前回までのエントリでトレーニングデータからパラメトリック曲線とそのパラメタの導出方法を書いた。 観測点(xn,tn)の背景にあるパラメトリック曲線を求める パラメタの評価式を決定する 評価式の値が最良になるようにパラメタを決定する 実際にはトレーニングデータはパラメトリック曲線から発生したデータではなく、パラメトリック曲線はあくまでトレーニングデータとの最小二乗誤差が最小になる曲線に過ぎない。 未知のサンプルデータの特徴をよりうまく含む方法を考える。 実際のサンプルデータが、トレーニングデータから導出したパラメトリック曲線から±σの誤差を持つデータであると仮定することで、サンプルデータの特徴を確率的に捉えることができる。これによりサンプルデータが「どれ位の範囲で当てはまるか」を推定できるようになる。 ITエンジニアのための機械学習理論入門posted with amazlet at 17.03.10中井 悦司 技術評論社 売り上げランキング: 8,130Amazon.co.jpで詳細を見る 正規分布の確率密度とパラメトリック曲線 統計・確率の基本中の基本、正規分布の確率密度。データxの発生頻度の平均がµ、分散がσ2である場合の確率密度は次の関数で与えられる。 $$ N(x|mu,sigma^2)=frac{1}{sqrt{2pisigma^2}}e^{-frac{1}{2sigma^2}(x-mu)^2} $$ で、今、M次多項式を考えている。 $$ begin{eqnarray} f(x) &=& w_0 + w_1x + cdots + w_Mx^M &=& sum_{m=0}^M w_mx^m end{eqnarray} $$ サンプルデータxnに対し観測値tがf(xm)を中心として±σ散らばることを考えると、tはfn,σ2から決まる式となる。 $$ begin{eqnarray} N(t|f(x_n),sigma^2)= frac{1}{sqrt{2pisigma^2}}e^{-frac{1}{2sigma^2}{(t-f(x_n))}^2} end{eqnarray} $$ トレーニングデータからパラメータを決定する 未知のサンプルデータ(xn,t)ではなく、トレーニングデータ(xn,tn)を代入すると、 $$ begin{eqnarray} N(t_n|f(x_n),sigma^2)= frac{1}{sqrt{2pisigma^2}}e^{-frac{1}{2sigma^2}{(t_n-f(x_n))}^2} end{eqnarray} $$ である。このモデルからトレーニングデータ(x0,t0)が得られる確率P0は以下の通り求まる。 $$ begin{eqnarray} P = frac{1}{sqrt{2pisigma^2}}e^{-frac{1}{2sigma^2}{t_0-f(x_0)}^2} end{eqnarray} $$ ここで、このモデルから全てのトレーニングデータ{(xn,tn)}n=1Nのデータが得られる確率Pは条件付き確率として以下の通り求まる。このPを尤度関数という。 $$ begin{eqnarray} P &=& N(t_1|f(x_1),sigma^2) N(t_2|f(x_2),sigma^2) cdots N(t_n|f(x_n),sigma^2) \\ &=& prod_{n=1}^N N(t_n|f(x_n),sigma^2) end{eqnarray} $$ f(xn)は係数行列wが未決定であり、分散σ2も未決定である。Pを最大化するようにf(xn)=wとσ2を決定する。積を分解して整理すると、 $$ begin{eqnarray} P &=& prod_{n=1}^N frac{1}{sqrt{2pisigma^2}}e^{-frac{1}{2sigma^2}{t_n-f(x_n)}^2} \\ &=& bigl(frac{1}{2pisigma^2}bigr)^{frac{2}{N}}expBigl[-frac{1}{2sigma^2}sum_{n=1}^N { t_n-f(x_n) }^2 Bigr] end{eqnarray} $$ ここで、尤度関数に二乗誤差EDが現れるので、EDで置き換える。 $$ begin{eqnarray} E_D &=& frac{1}{2} sum_{n=1}^N{f(x_n)-t_n}^2 P &=& Bigl( frac{1}{2pisigma^2} Bigr)^{frac{N}{n}} e^{-frac{1}{sigma^2}E_D} end{eqnarray} $$ 尤度関数がσとwから決まる2変数関数であることを明確にするため、σ^2の逆数をβとし、二乗誤差EDが係数wから決まることを明確とするためED=ED(w)とする。 $$ P(beta,w) = bigl( frac{beta}{2pi} bigr)^{frac{N}{2}} e^{-beta E_D(w)} $$ 両辺に自然対数を取る。Pが最大になることとPの自然対数が最大になることは同じ。 $$ begin{eqnarray} ln{P(beta,w)} &=& ln{bigl( frac{beta}{2pi} bigr)^{frac{N}{2}} e^{-beta E_D(w)}} \\ &=& frac{N}{2}ln{beta}-frac{N}{2}ln{2pi}-beta E_D(w) end{eqnarray} $$ 単調増加関数が最大になる条件は微分した値が0になること。2変数の単調増加関数が最大になる条件は、それぞれの変数で偏微分した値が0になること。つまり、 $$ begin{eqnarray} frac{partial (ln{P})}{partial w_m} &=& 0 \\ frac{partial (ln{P})}{partial beta} &=& 0 end{eqnarray} $$ まず、wの偏微分。EDw項以外は定数になるから、 $$ begin{eqnarray} frac{partial E_D}{partial w_m} = 0 \\ end{eqnarray} $$ パラメトリックモデルに確率密度関数の誤差を含まない最小二乗法における二乗誤差を最小にする条件と同じとなる。つまり、トレーニングデータ{xn}n=0N},{t1,...,tN}から求まる。 $$ begin{eqnarray} w = bigl(Phi^TPhibigr)^-1Phi^T t end{eqnarray} $$ 次に、βの偏微分。β項以外は定数になる。対数の微分は、(logx)\'=1/x だから、 $$ begin{eqnarray} frac{partial P(beta,w)}{partial beta} &=& 0 \\ frac{N}{2} frac{1}{beta} - E_D(w) &=& 0 \\ frac{1}{beta} &=& frac{2E_D}{N} end{eqnarray} $$ ここで、以下としてしたから、 $$ beta = frac{1}{sigma^2} $$ 式変形すると、なんとσと、最小二乗法でパラメトリックモデルの次数Mを評価したときの、モデルとサンプルとの差(平方根平均二乗誤差)ERMSと同じになる。 $$ sigma = sqrt{frac{1}{beta}} = sqrt{frac{2E_D}{N}} = E_{RMS} $$ つまり... 最初に戻ると、トレーニングデータ{xn}が与えられたとき、トレーニングデータがM次多項式f(x)から±σの範囲にバラついていると仮定した場合、モデルから{tn}が得られる確率は、2つのパラメタw,σの関数として以下のように決まる。 $$ begin{eqnarray} N(t_n|f(x_n),sigma^2)= frac{1}{sqrt{2pisigma^2}}e^{-frac{1}{2sigma^2}{(t_n-f(x_n))}^2} end{eqnarray} $$ 全てのトレーニングデータがこのモデルから得られる条件付き確率(尤度関数)は条件付き確率により定まり、尤度関数の最大化問題を解くことで2つのパラメタw,σを得られる。 wは、モデルに含まれる正規分布に従う誤差を含まないで作成したモデルの場合と同じ方法で求められる。また、σはそのモデルとの誤差評価値ERMSと同じである。 最尤推定法は最小二乗法を一般化した話だった、というか最小二乗法がサンプルデータの生起確率が正規分布を前提としていた、という話。 ともかく、トレーニングデータから確率変数付きパラメトリックモデルを立てて、パラメタを確率の最大化問題として解くことを最尤推定という。 サンプルデータの生起確率が正規分布に従う例がわかりやすい。より一般的に、観測した事象から確率的モデルを立てることは統計的なモデル化の王道らしく、そこから派生している機械学習アルゴリズムでも頻出であるようだ。