default eye-catch image.

決定木の分割基準,情報ゲイン,エントロピー

[mathjax] 集合に対して再帰的に境界を入れていく操作が決定木の作成。 では、集合のどこに境界を入れれば良いか。 属性をテストすることにより得られる情報量が最も大きくなるように入れる。 汎化能力、みたいな言葉を読んでいくにあたってこの先、結構抽象的な話になるので一度確認。 データが(n)個のクラスに分類できるとして、クラス(i)に属する確率を(P_i)とする。 このとき、あるデータがクラス(i)に属することを知るには(log_2 frac{1}{P_i})の情報量が必要。 その期待値は(I(P_1,P_2,cdots,P_n)=sum_{i=1}^{n} P_i log_2 frac{1}{P_i} = - sum_{n}^{i=1} P_i log_2 P_i)。 情報量の平均をエントロピーとか。 (D_p)個のデータを(D_{left})、(D_{right})に分割するとする。 その時、属性(f)に関する問いを使って分割する。 2分木の左ノード、右ノードだからleft,right。 (I(D_p))は分割前のエントロピー。 (I(D_{left}))、(I(D_{right}))は分割後のエントロピー。 分割前には(I(D_p))ビット必要だったのが、 分割後には(frac{N_{left}}{N_p} I(D_{left}) + frac{N_{right}}{N_p} I(D_{right}))ビット必要になった。と読むらしい。 その差を情報ゲイン(IG(D_p))と呼んで以下のように定義するらしい。 begin{eqnarray} IG(D_p,f) = I(D_P) - frac{N_{left}}{N_p} I(D_{left}) - frac{N_{right}}{N_p} I(D_{right}) end{eqnarray} 分割前よりも分割後の方が(IG(D_p,f))だけエントロピーが低い、という事実に関して、 分割に使った問いにより、(IG(D_p,f))の情報量を獲得した、と考えるらしい。 情報ゲインが最大になるような問いを根にもってきて、 再帰的に情報ゲインが大きいものから順に問うことで決定木を作っていく。

default eye-catch image.

交差検証(CrossValidation)

同じ出処から取ってきたデータを全て訓練データとして使わずに、 訓練データとテストデータに分割して、訓練データで作ったモデルに対するテストデータの精度を返す、 みたいなことをやるらしい。交差検証(CrossValidation)という名前が付いている、 sklearn.model_selection.cross_val_score( estimator, X, y=None, groups=None, scoring=None, cv=’warn’, n_jobs=None, verbose=0, fit_params=None, pre_dispatch=‘2*n_jobs’, error_score=’raise-deprecating’) estimatorとして、例えば決定木分類であればDecisionTreeClassifierのインスタンスを渡す。 Xは説明変数、yは目的変数。交差検証自体のアルゴリズムを選択できてcvに渡す値で制御できる。 学習済みのモデルを渡してスコアが戻るのではなく、 単なるモデルのインスタンスと学習前のデータを別々に渡している作りなのを見ると、 モデル毎のスコアを並べてみたくなる..。そういう使い方するのか? \"スコア\"と言ってるこの値は具体的には何の値なのか..? K-分割交差検証(K-fold cross-validation)の場合cvは以下のようになる。 cvを省略するとK=3が使われる。使ったデータは例のあやめ。 clf = DecisionTreeClassifier(max_depth = 3) # integer, to specify the number of folds in a (Stratified)KFold, score1 = cross_val_score(estimator = clf, X = iris.data, y = iris.target, cv = 10) # None, to use the default 3-fold cross validation, score2 = cross_val_score(estimator = clf, X = iris.data, y = iris.target) 決定木の深さを交差検証で求める(失敗) suumoから引いてきた家賃~占有面積,築年数,階数,バストイレ別データについて、 目的変数を家賃、説明変数を占有面積、築年数、階数、バストイレ別として決定木を作ってみた。 決定木の深さは交差検証で求める、と書かれているので以下のようにしてみた。 import sys import pandas as pd import matplotlib.pyplot as plt from sklearn import tree from sklearn.model_selection import train_test_split from sklearn.model_selection import cross_val_score path = \"./fuchu.csv\" train = pd.read_csv(filepath_or_buffer=path) train = train.drop_duplicates() feature_name = [\"area\",\"age\",\"bt_separated\"] feature_name2 = [\"area\"] train_x = train[feature_name2].values train_y = train[\"value\"].values sort_idx = train_x.flatten().argsort() max_depth_list = [1,3,5,7,9,11,13,15,17] for max_depth in max_depth_list: regressor = tree.DecisionTreeRegressor(max_depth = max_depth) regressor.fit(train_x,train_y) score = cross_val_score(estimator=regressor,X=train_x,y=train_y) print([max_depth,score.mean()]) 結果、おかしい...。マイナスってなんだ。 [1, -0.020456885057637583] [3, 0.012376893473610817] [5, -0.014519973352621932] [7, -0.03332646095737912] [9, -0.06862181989491711] [11, -0.09590667631145984] [13, -0.1345799052512994] [15, -0.1529266053016934] [17, -0.17054978124438266] 説明変数を占有面積だけに絞って散布図書いて、その上に回帰結果をプロットすると、 異常値の0に引っ張られて乱れてた...。以下は深さ=17の時。むぅ. 真ん中の塊の中を左下から右上に向かってギザギザに進んで欲しいのだが、 物凄い引っ張られよう。異常値が原因というよりは、引っ張られ過ぎなんだな。

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.

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

[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で詳細を見る