カテゴリー: 教養

[mathjax] 機械学習の分類問題の中心にある決定境界の決定方法について かなり要領を得た説明を聞いて理解が2段階くらい先に進んだのでまとめてみます。 データが与えられただけの状態から決定境界を決める問題はNP困難ですが 別の問題に帰着させることで解を得る、というのが基本的なアイデアです。 分類の正誤とその度合いを一度に表現できるマージンを定義し、 マージンを使って与えた代理損失を最小にする問題にします。 分類問題を代理損失の最小化に帰着させるのですね。 任意の決定境界を決める問題は線形分類であってもNP困難 2値のラベルA,B付きの2次のデータポイントが与えられたとして、 入力空間(X1-X2)におけるA,Bの分離境界(decision boundary)を求める問題が\"分類\"。 直線で分離境界を書くとして、それを求めるための最も愚直な方法は以下のようなもの。 その分離境界によりデータポイントが正しく分類出来ていれば1をカウントする。 正しく分類出来ていなければ0をカウントする。 全データポイントにおける正答率を求める。 正答率が最大になるような決定境界を求める。 そもそも分離境界は直線でなくても良いのに、あえて直線ですよ、と仮定をしたとしても、 分離境界が完全に自由で、全データに対して正答率を求めないといけない。 上記の問題の計算量は(mathcal{O}(n^3))では済まない。NP困難。 計算できるように改善 分離境界の初期値を決めて、そこから正答率が良くなる方向に少しずつずらしていこうにも、 \"正しく分類されている\"=1,\"分類されていない\" =0 は、少しの変化に影響されない。 正しい=1/正しくない=0、という損失とは別の損失を作って、 その損失を使った別の問題を解くことを、上記を問題を解くことに帰着させる。 決定境界の変化に敏感な損失を作る サンプルサイズが十分大きいとき、1.で作った損失による学習結果が、「正しく分類」「正しくない分類」という損失の学習結果と一致する margin 線形分類において、分離境界(f(x_1,x_2,cdots,x_n)=w_0+w_1x_1+w_2x_2+cdots+w_nx_n)とする。 この多項式と分離の正誤、正誤の度合いは以下のように決まる。 分類の正誤は(f(x_1,x_2,cdots,x_n))の符号が決める。 分類の正誤の度合いは(f(x_1,x_2,cdots,x_n))の絶対値が決める。 (f(x_1,x_2,cdots,x_n))が正の場合、決定境界から近い場所にあるデータポイントは もしかしたら誤って分類してしまったものかもしれない。 決定境界から遠い場所にあるデータポイントは近いものよりは正しく分類しているかもしれない。 同様に(f(x_1,x_2,cdots,x_n))が負の場合、決定境界から近い場所にあるデータポイントは もしかしたら正しい分類かもしれないし、遠いデータポイントはより近いものより間違っている可能性が高い。 この事実を1つの式で表す。 データポイントには出力ラベル(y=pm 1)が付いているものとする。 判別関数を(f(x_1,x_2,cdots,x_n))とする。決定境界は(f(x_1,x_2,cdots,x_n)=0) begin{eqnarray} m = yf(x_1,x_2,cdots,x_n) end{eqnarray} ラベル1を-1と分類した場合(f(x_1,x_2,cdots,x_n)<0)。 同様に-1を1と分類した場合も(f(x_1,x_2,cdots,x_n)<0)。 つまり、誤分類したときにラベルと判別関数の符号が異なり(m0)となる。 ということで、(m)をマージン(margin)と呼ぶ。 サポートベクトル marginが最大になるように各データポイントの中にある決定境界を決めていく。 全てのデータポイントについて距離を計算する必要はなく、決定境界と距離が一番近いデータポイントとの 距離を最大化すれば良いらしい。（それが一番近いかどうかはいずれにせよ距離を求める必要がありそうだけど..) marginが最大になるように決めた決定境界と距離が最も近いデータポイントをサポートベクトルと言うらしい。 マージンを使った損失 最初に戻ると、決定境界の変化に敏感な損失を作ることが目的だった。 マージンが正の方向に大きいほど正しい分類であると言えるし、 マージンが負の方向に大きいほど誤った分類であると言えるけれども、 正しい度合いが高ければ小、誤りの度合いが高ければ大、となる損失を考えることで、 誤った方向に決定境界を修正すれば敏感に値が上昇する損失にすることができる。 (正しい方向に移動しても変わらない。) 横軸にマージン、縦軸に損失を取ったとして、以下のような損失(h(m))を考える。 もちろん、(m = yf(x_1,x_2,cdots,x_n))。 (m=1)より大きいマージンについては損失が0。(m=1)より小さいマージンについて線形に増加する。 (m=1)を境にヒンジの形をしているのでhinge損失という名前が付いてる。 begin{eqnarray} h(m) = max(0,1-m) end{eqnarray}

[mathjax] 前の記事で線形単回帰において訓練データから回帰係数を求める方法を書いてみた。 標本平均を使って母平均を推測する話とリンクさせることで、 回帰係数の95%信頼区間を求めることができた。 回帰係数(hat{beta_0},hat{beta_1})と真の回帰係数(beta_0,beta_1)の関係がこれ。 [clink url=\"https://ikuty.com/2019/05/15/linear_regression_evaluate/\"] RSE,真の回帰直線と観測データがどれくらい離れているか 真の回帰直線がわかったとしても、全てのデータが回帰直線の上に乗っているのでなければ、 回帰直線を使って値を予測したときに誤差が出てくる。 残差平方和(Residual sum of square)。WikipediaにもRSS。 (hat{y_i})は訓練データを使って得られた回帰係数で作った回帰直線で予測した値。 だから、RSS自体も訓練データに対応して変動する。 begin{eqnarray} RSS=sum_{i=1}^n (y_i-hat{y_i})^2 end{eqnarray} で、知りたいのはRSSが訓練データに対してどの程度変動するかだから標準偏差。 標本分散は不偏推定量ではなくて分布の自由度で割る必要がある...という話があって、 不偏推定量を求める段取りが必要。(n-1)ではなく(n-2)で割る!。詳しくは以下。 カイ2乗分布になりそうだけれども、自由度が何故(n-2)なのだろうか...。 begin{eqnarray} RSE= sqrt{frac{1}{n-2}sum_{i=1}^{n}(y_i-hat{y_i}^2)} end{eqnarray} [clink implicit=\"false\" url=\"https://stats.stackexchange.com/questions/204238/why-divide-rss-by-n-2-to-get-rse/377869\" imgurl=\"https://cdn.sstatic.net/Sites/stats/img/logo.svg?v=60d6be2c448d\" title=\"Why divide RSS by n-2 to get RSE?\" excerpt=\"The reason is based on trying to get an unbiased estimator of the underlying error variance in the regression. In a simple linear regression with normal error terms it can be shown that:That is, under the standard assumption of normally distributed errors, the residual sum-of-squares has a chi-squared distribution with ?−2 degrees of freedom. \"] 決定係数(R^2) RSS,RSEは(Y)の単位で値が決まる。(y_i)が無茶苦茶大きいとRSEは大きくなる。 RSEだけ見て回帰直線がどれだけ当てはまっているか言えない様子。 当てはまりの良さを(0)から(1)の範囲におさめる別の指標もある。 TSS (Total sum of square)として以下。 begin{eqnarray} TSS = sum_{i-1}^{n}(y_i-bar{y_i})^2 end{eqnarray} (R^2)として以下。 begin{eqnarray} R^2 &=& frac{TSS-RSS}{TSS} \\ &=& 1-frac{RSS}{TSS} \\ &=& 1 - frac{sum_{i=1}^n (y_i-hat{y_i})^2}{sum_{i-1}^{n}(y_i-bar{y_i})^2} end{eqnarray}

[mathjax] 線形単回帰で推定する回帰係数の精度を評価する方法を読んだのでまとめてみる。 当然、真の直線はわからないのだけれども、真の直線があると仮定した上で 推定した回帰係数との関係を考えることで、回帰係数の精度について話せるようになる。 回帰係数の導出 データポイントが(n)個ある状況。 ( (x_1,y_1),(x_2,y_2),cdots,(x_n,y_n) ) 回帰係数(hat{beta_0})と(hat{beta_1})を使って線形回帰したい。 begin{eqnarray} hat{y} = hat{beta_0} + hat{beta_1} x end{eqnarray} データポイントと回帰直線の差を残差平方和(RSS,redisual sum of square)で表す。 データポイントは既に与えられているデータなので、(hat{beta_0},hat{beta_1})の関数。 begin{eqnarray} f(hat{beta_0},hat{beta_1}) = (y_1 -hat{beta_0}-hat{beta_1}x_1)^2 + (y_2 - hat{beta_0}-hat{beta_1}x_2)^2 + cdots + (y_n - hat{beta_0}-hat{beta_1}x_n)^2 end{eqnarray} RSSを最小にする(hat{beta_0})と(hat{beta_1})を求めるために、(hat{beta_0})、(hat{beta_1})それぞれで偏微分して(0)として解く。 なんでそれぞれ個別に偏微分して0と置いて良いかは、 RPML読もうとして力尽きたときに理解したので省略。 参考にした本に( hat{beta_0}),(hat{beta_1}),RSSの3次元の図があって、確かにそれで良さそうな予感。 begin{eqnarray} frac{partial}{partial hat{beta_0}} f(hat{beta_0},hat{beta_1}) = 0 \\ frac{partial}{partial hat{beta_1}} f(hat{beta_0},hat{beta_1}) = 0 \\ end{eqnarray} 以下のようになるらしい。(bar{x})、(bar{y})はデータポイントの標本平均。 なので、データポイントがわかれば計算で求まる。 begin{eqnarray} hat{beta_1} &=& frac{sum_{i=1}^n (x_i-bar{x}) (y_i-bar{y}) }{sum_{i=1}^n (x_i-bar{x})^2 }\\ hat{beta_0} &=& bar{y}-hat{beta_1}bar{x} end{eqnarray} 母回帰直線の推定 データポイントが同じであれば(hat{beta_0}),(hat{beta_1})は同じになるけれども、 データポイントを取り直して異なるデータセットにすると、(hat{beta_0}),(hat{beta_1})は微妙に違う値になる。 じゃあ、データセットを大量に用意したとして、(hat{beta_0}),(hat{beta_1})を計算しまくると、 どこかに収束するんじゃなかろうか。 標本が大量にあると標本平均は母平均に収束する。標準偏差はより小さくなる。 つまりデータが大量にあると、母平均からのズレが小さくなっていく。 大数の弱法則、中心極限定理、ルートnの法則。 begin{eqnarray} hat{sigma} &=& frac{sigma}{sqrt{n}} \\ hat{sigma}^2 &=& frac{sigma^2}{n} end{eqnarray} begin{eqnarray} lim_{n rightarrow infty} hat{sigma}^2 = lim_{n rightarrow infty} frac{sigma^2}{n} = 0 end{eqnarray} [clink url=\"https://ikuty.com/2018/07/17/sample_sigma/\"] (hat{beta_0}),(hat{beta_1})は母回帰直線からどれくらいばらついているのか。 (hat{beta_0}),(hat{beta_1})の分散は以下を使うらしい。 両方に出てくる(sigma^2)は、母回帰直線と回帰直線の差となる項の散らばり度合い。 つまり、(Y=beta_0 + beta_1 X + epsilon )としたときの(epsilon)の分散。 begin{eqnarray} sigma_{hat{beta_0}}^2 &=& sigma^2 Bigl[frac{1}{n} + frac{bar{x}^2}{sum_{i=1}^n (x_i-bar{x})^2} Bigr] \\ sigma_{hat{beta_1}}^2 &=& frac{sigma^2}{sum_{i=1}^n (x_i -bar{x})^2} end{eqnarray} (x_i)が散らばれば散らばるほど、(sigma_{hat{beta_1}}^2)は小さくなる。 データポイントの(x)成分が小さい方から大きい方まで含まれれば、傾き(beta_1)を推定しやすくなる。 そして、(bar{x}=0)であるならば、(hat{beta_0})の散らばりは、(hat{mu})の散らばりと等しくなる。 最終的に求めたいのは不明な(sigma^2)だが、(sigma^2)はデータから計算できる。 (sigma)の推定値(RSE,Resual Standard Error)はRSSから推定する。 begin{eqnarray} sqrt{frac{f(hat{beta_0},hat{beta_1})}{(n-2)}} end{eqnarray} (hat{beta_1})の標準偏差がわかったので、95%信頼区間を求めることができる。 線形回帰における(hat{beta_1})の95%信頼区間は、 begin{eqnarray} Bigl[ hat{beta_1} - 1.96 sigma_{hat{beta_1}},hat{beta_1} + 1.96 sigma_{hat{beta_1}} Bigr] end{eqnarray} 同様に(hat{beta_0})の95%信頼区間は、 begin{eqnarray} Bigl[ hat{beta_0} - 1.96 sigma_{hat{beta_0}},hat{beta_0} + 1.96 sigma_{hat{beta_0}} Bigr] end{eqnarray}

[mathjax] 損失関数をバイアス項、バリアンス項、削減不能誤差の和に分解できることと、 損失は削減不能誤差より下回らないこと、バイアス項、バアリアンス項のトレードオフが起こること、 を読んだ。過学習っていうのはこういうことなのか、と腑に落ちたので記念に記事を書いてみる。 （式変形は細かいところで間違ってるのと、おっさんのチラシの裏なので参考にしないでください） 2乗損失の期待値の式変形 モデルを作った後、訓練データ、テストデータそれぞれの全データについて、 2乗損失の期待値(MSE)を求め、モデルの当てはまりの良さを調べるらしい。 2乗損失を以下のように式変形する。条件付き期待値(E(t|x))ってなんだ...。 begin{eqnarray} L(y(x),t)^2 &=& (y(x)-t)^2 \\ &=& (y(x)-E(t|x)+E(t|x)-t)^2 \\ &=& left( left( y(x)-E(t|x)right) + left( E(t|x) - y(x)right) right)^2 \\ &=& (y(x)-E(t|x))^2 + 2(y(x)-E(t|x))(E(t|x)-t) + (E(t|x)-t)^2 \\ end{eqnarray} 2乗損失の期待値(MSE)は以下。 第2項は(x)、(t)で積分するとゼロになる! begin{eqnarray} E[L(y(x),t)^2] &=& E[ (y(x)-E(t|x))^2 + 2(y(x)-E(t|x))(E(t|x)-t) + (E(t|x)-t)^2 ] \\ &=& E[ (y(x)-E(t|x))^2 + (E(t|x)-t)^2 ] end{eqnarray} 和の期待値は期待値の和なので、 begin{eqnarray} E[L(y(x),t)^2] = E[ (y(x)-E(t|x))^2 ] + E[(E(t|x)-t)^2 ] end{eqnarray} (x)の出処がテストデータではなく訓練データですよ、と明示するために、 以下みたいな書き方に改める。この式の中で(y(x;D))が学習で得られるモデル。 第2項は学習とは関係なく発生する数値。 begin{eqnarray} E_D[L(y(x;D),t)^2] &=& E_D[ (y(x;D)-E(t|x;D))^2 ] + \\ && E_D[(E(t|x;D)-t)^2 ]] end{eqnarray} 第1項の式変形を続ける。括弧が多すぎて力尽きた..。 余計な項を足して引いて次の式変形の足しにするタイプ。 begin{eqnarray} E_D[ (y(x;D)-E(t|x;D))^2 ] &=& E_D[ ( { y(x;D)-E_D[(y(x;D))] } ] &+& { E_D[ y(x;D)] - E[t|x;D])^2 } \\ &=& E_D [ { (y(x;D))-E_D[ y(x;D)] }^2 ] + \\ &=& E_D [ { E_D[ y(x;D)-E[t|x;D] ] }^2 ] end{eqnarray} バイアス・バリアンスと削減不能誤差 以下はバリアンス項と書かれている。 モデル((y(x;D))による予測が訓練データ集合によって変動する度合いの期待値。 異なる訓練データを使ったときにどの程度モデルが変化するかを表す。 過学習の度合い。 begin{eqnarray} E_D bigl[ bigl{ (y(x;D))-E_D[ y(x;D)] bigr}^2 bigr] end{eqnarray} 以下はバイアス項と書かれている。 複雑な事象を単純なモデルで近似したことによる誤差、と書かれてる。 例えば、3次関数+ノイズから発生するデータを直線で近似すると、モデルが単純すぎて値が大きくなる。 モデルが複雑になればなるほどバイアス項は減っていく様子。 未学習の度合い。 begin{eqnarray} E_D bigl[ bigl{ E_D[ y(x;D)]-E[t|x;D] bigr}^2 bigr] end{eqnarray} で、一番最初に出てきたモデルと関係ない以下。 バイアス、バリアンス共に非負の値だから、2乗損失の期待値は以下より小さくなることはない。 奇跡的にバイアス、バリアンス共にゼロだったとしても、以下は学習とは関係なく発生する。 削減できない誤差。 begin{eqnarray} E_Dbigl[bigl(E(t|x;D)-tbigr)^2 ]bigr] end{eqnarray} 結局よくわからない...。体感の結論.. 訓練データを使ってモデルを複雑にしていけばいくほど、 モデルが訓練データにフィットするようになるが、 その訓練データにフィットしまくったモデルは、未知のテストデータを予測しづらくなる。 モデルの複雑度が\"ある程度\"のところまでは、バリアンスの上昇よりもバイアスの低下が効くから、 訓練データに対する2乗誤差、テストデータに対する2乗誤差ともに減少する。 モデルの複雑度が\"ある程度\"を超えると、バイアスの低下が頭打ちになる一方でバリアンスが上昇し、 訓練データに対する2乗誤差が低下する一方で、テストデータに対する2乗誤差が上昇する。 どう頑張っても、削減不可能な誤差が存在する。 条件付き期待値(E(t|x))の意味を理解できずプロットすることは叶わなかった。

[mathjax] おっさんが入門した軌跡シリーズです。損失関数に関して学んだことをメモしておきます。 入力値(x)、正解(t)からなる訓練データ(T=(x,t))が大量に与えられたときに、 (f(x,w))によって回帰なり分類なりをする。 仮に立てた(f(x,w))と正解(t)の距離(Lleft(f(x,w),tright))を損失関数と呼んだり。 あえて(Lleft(f(x,w),tright))としているのは、 一番わかりやすそうな残差の2乗だけでなく、他があるから。 2乗損失 残差の2乗だと、(f(x,w))と(t)の差が大きい場合に必要以上に大きくなってしまう。 ほとんどのデータで残差が0なのに、特殊なデータで残差が100とかになられたら、 全体の損失は測れそうにないし、異常値に敏感すぎる。 begin{eqnarray} Lleft(f(x,w),tright) = (t-f(x,w))^2 end{eqnarray} Huber損失 単に残差の2乗を使うだけでは不十分で、\"(f(x,w))と(t)の差の大小にあまり影響されないこと\"が必要。 残差の絶対値がある値を超えるまでは残差の2乗、 超えてからは線形にするというのもある(Huber損失)。 begin{eqnarray} Lleft(f(x,w),tright) = begin{cases} (t-f(x,w))^2 & t in [t-delta,t+delta] \\ 2delta (|t-f|-frac{delta}{2}) & それ以外 end{cases} end{eqnarray} 損失関数が微分可能か 損失関数を最小(極小)にすることが目的なので..、 損失関数の1階導関数(勾配ベクトル)を使ってパラメタを更新したりする。 begin{eqnarray} w_{i+1} = w_i - eta left(f\'(x,w),tright) end{eqnarray} 損失関数が微分可能（連続）だと勾配ベクトルがすぐに求まるので、 損失関数は不連続点をなくすように作るらしい。 (Huber損失の2乗から線形に切り替わるところは連続になってる。) 次回以降実際のデータとモデルを使ってやってみる。

[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))の情報量を獲得した、と考えるらしい。 情報ゲインが最大になるような問いを根にもってきて、 再帰的に情報ゲインが大きいものから順に問うことで決定木を作っていく。

同じ出処から取ってきたデータを全て訓練データとして使わずに、 訓練データとテストデータに分割して、訓練データで作ったモデルに対するテストデータの精度を返す、 みたいなことをやるらしい。交差検証(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の時。むぅ. 真ん中の塊の中を左下から右上に向かってギザギザに進んで欲しいのだが、 物凄い引っ張られよう。異常値が原因というよりは、引っ張られ過ぎなんだな。

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

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とか、間取りの区別はしていないから同時に出る。 家賃と占有面積は関係があるということがわかる様子。 広い方が高いよ、は言える様子。 他の因子に引きずられず、築浅であろうが駅近であろうが、広さは正義らしい。 家賃と築年数は、全体的に右肩下がりだけども、 占有面積ほど関係があるとは言えない様子。 現実的には、築浅だとしても駅から遠ければお安い、という話があると思い..。 全体的に右肩下がりなので、そうであっても古ければ概ね安いという。 もっと意外なのは、家賃と最寄駅からの時間。 駅近であることが一番影響してそうに思ってたけどもそんなことはなかった。 確かに、駅近であっても古くて狭いのは安いからな。 占有面積と築年数ほど、家賃に影響を与えないらしい。 思い込みというのは怖いもの..。 多変量解析なら、どれの係数が高いかはここで当たりがつくけれども、 ランダムフォレストで家賃予測するとすると、どうなるか、次回に続く。

[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)が増加していくのがわかる。