カテゴリー: 機械学習

[mathjax] 階層型クラスタリングについてわかりやすい解説を聞いたので頭の中にあるものを書き出してみます。 (せっかく聞いたシャープな解説が台無しになってしまっているかと思いますが...) 本当はもっとアクティブなアウトプット方法があれば良いのですが、今のところはブログしかないのと、 WordPressのエントリという単位が、アウトプットの量として個人的に丁度良いなと思いますので。 初学なので解釈の誤りとか検討不足があるかと思いますが、本日も書いていきます！ 試しに実行してみるためのデータを作る面倒さの洗礼を受けまして、 実際に試すのは別記事になります。 [clink url=\"https://ikuty.com/2019/07/16/hierarchical_clustering_by_method/\"] 階層型クラスタリング 小さいクラスタを大きいクラスタで内包する様が木構造(階層構造)になることから、 階層型クラスタリング、という名前が付いているんでしょうかね。 最初にk個のクラスタに分割する！ みたいな問題であればk-means法、 何個のクラスタに分割すべきかはやってみて考えるのであれば、階層型クラスタリング。 木の深さとクラスタ数を対応させるところに見事な美しさを感じます。 クラスタ間距離として以下のWard距離を使うことで空間拡張性と非単調性に優れた結果を得られやすく、 よく傾向が分かっていないデータに対して一発目に適用すべき方法の様子。 begin{eqnarray} d_{ward} (C,C\') = sum_{iin Ccup C\'} || x_i - bar{x}_{Ccup C\'}||^2 - left( sum_{i in C} ||x_i-bar{x}_C||^2 + sum_{iin C\'}||x_i - bar{x}_{C\'}||^2 right) end{eqnarray} 階層型クラスタリングのアルゴリズム 以下の通り。最初、全てのデータポイントをクラスタとみなすところからスタートする。 クラスタとクラスタを併合してより大きいクラスタを作る処理であって、 繰り返しを止めなければ最終的に一つのクラスタまで成長する。 つまり、どこかで止めることで、その時点のクラスタ分割を得ることができる。 全データポイントそれぞれをクラスタとする 全てのデータポイントが初期クラスタ以外に含まれるまで以下の計算を繰り返す 各クラスタ間の距離を計算する。 最もクラスタ間距離が小さいクラスタを併合する。 クラスタを成長させる際に、うまい具合に併合相手を選ばないとよくない。 全体としてクラスタがバラけて欲しいのだけれども、クラスタが隣のクラスタを併合した結果、 次の併合先としてその隣のクラスタが対象となってしまうかもしれない。 そうすると、単調に端っこから順番に成長することになり、 「大きなクラスタ」と「その他の小さなクラスタ」みたいになる。 何らかの距離に基づいてクラスタとクラスタを併合するのだけれども、 クラスタを併合した結果、併合後のクラスタと別のクラスタの距離が、 より短くなってしまう、というケースも起こってしまう。 クラスタ内のデータポイント間の距離、とか、クラスタの重心間の距離とかでは 上記の問題が起こりやすい。 先に距離があって次に併合を考えるのではなく、併合前後の条件を使って距離を作ることで、 上記の問題が起こりにくくなる。その一つがWard距離。 Ward距離 データポイント間の距離はユークリッド距離を使うとして、クラスタ間の距離をどう決めるのか？ データポイントをクラスタ(C)とクラスタ(C\')に分けるとする。分ける前は(Ccup C\')。 (逆? クラスタ(C)とクラスタ(C\')を併合してクラスタ(Ccup C\')を作る)。 クラスタ(C)とクラスタ(C\')の距離(d_{ward}(C,C\'))を以下のように定義する。 begin{eqnarray} d_{ward} (C,C\') = sum_{iin Ccup C\'} || x_i - bar{x}_{Ccup C\'}||^2 - left( sum_{i in C} ||x_i-bar{x}_C||^2 + sum_{iin C\'}||x_i - bar{x}_{C\'}||^2 right) end{eqnarray} 併合後のクラスタ(Ccup C\')における重心との距離の2乗和( sum_{iin Ccup C\'} || x_i - bar{x}_{Ccup C\'}||^2)から、 併合前のクラスタ(C)における重心との距離の2乗和(sum_{i in C} ||x_i-bar{x}_C||^2)、 (C\')における重心との距離の2乗和(sum_{iin C\'}||x_i - bar{x}_{C\'}||^2)を引いたもの。 当たり前だけれども、大きいクラスタよりも小さいクラスタの方が凝集度が小さい。 つまり、併合前のクラスタ(C)、(C\')を併合してクラスタ(Ccup C\')を作るとなると、 (Ccup C\')の凝集度よりも、(C)、(C\')の凝集度の方が必ず小さい。 (C)の凝集度と(C\')の凝集度の和よりも、(Ccup C\')の凝集度の方が大きくなることがポイント。 要するに(Ccup C\')を(C)、(C\')に分割することは、全体の凝集度が小さくなるような操作を行うということ。 無数にある分割の仕方の中から、凝集度の総和が小さくなるような「マシな分割」を決めたい！ (sum_{i in C} ||x_i-bar{x}_C||^2+sum_{iin C\'}||x_i - bar{x}_{C\'}||^2)よりも( sum_{iin Ccup C\'} || x_i - bar{x}_{Ccup C\'}||^2)の方が大きいから、 (d_{ward})は正の値となる。(d_{ward})が大きければそれはマシな分割。 そもそも、データ達が何をもってクラスタを形成しているのかって考えるとかなり微妙で、 人間が解釈するための合理的な根拠を作り出す手段でしかないことがわかった。 データポイント間の\"距離\"を何らかの方法で定義し、\"距離が近しい塊をクラスタと言う\" という話であるとすれば、確かに合理的なのだろうけども。 モヤモヤしているのは、k個のクラスタ集合とk+1個のクラスタ集合を比較したときに、 前者と後者で異なるクラスタに属するデータポイントが発生するところ。 クラスタリングとは単に境界を引く処理なのであって、 それにどういう意味を付けるかは目的次第なのだろうと改めて認識した次第。 そもそも、これらは単なる道具でしかないので、他の手段も同じく 出力結果にどう意味付けするのかは目的次第なんだろうなと。 最初から何個のクラスタに分割すべきか分かっているときはk-means法を使えば良さそう。 今回は、そもそも何個に分割すべきか分かっていないときに使う手法。 [clink url=\"https://ikuty.com/2019/06/29/k-means_idea/\"]

[mathjax] 教師なし学習の問題。クラスタリング。 クラスタの個数を事前に指定するタイプと、自分でクラスタ数を設定できるタイプがあります。 今回、前者のk-means法をアイデアを聞いたので、まとめなおしてみようと思います。 k-means法 (C_1,cdots,C_K)がクラスタ。各々のクラスタにおける中心は(bar{x_k})(centroid)。 各々のクラスタにおけるcentroidとの距離の2乗を全部足して、 それが最小になるようなクラスタ(C_1,cdots,C_K)を求めるという問題。 begin{eqnarray} newcommand{argmin}[1]{underset{#1}{operatorname{arg},operatorname{min}};} hat{c_1},cdots,hat{c_K} = argmin {c_1,cdots,c_K} sum_{k=1}^K sum_{i in C_k} || x_i - bar{x_k}||^2 end{eqnarray} 距離の2乗なので異常値があると異常値に引っ張られて、クラスタが大きく変わってしまう。 centroid(標本平均)も異常値に引っ張られる。外れ値を除去しないといけない。 距離の2乗を計算するので、ベクトル(v_i)の各次元のスケールが揃っていないと、 計算の過程で単位が大きい要素に引っ張られてしまう。 単位の大きさは重要度とは関係ないのに、単位が大きいことが距離に影響を与えることは 避けないといけない。その対応として単位を揃える（正規化）する。 begin{eqnarray} |x_i-bar{x}|^2 = begin{pmatrix} x_{i1} - bar{x_1} \\ x_{i2} - bar{x_2} \\ vdots \\ x_{in} - bar{x_n} \\ end{pmatrix} = sum_{i=1}^{n} sqrt{x_{i1}-bar{x_1}}^2 end{eqnarray} そのままだとNP困難->近似解法 データポイントの個数、centroidの位置がそれぞれ自由だと、 定義通りに実装しようとしても(mathcal{O}(n^3))となって使えない。 begin{eqnarray} newcommand{argmin}[1]{underset{#1}{operatorname{arg},operatorname{min}};} hat{c_1},cdots,hat{c_K} = argmin {c_1,cdots,c_K} sum_{k=1}^K sum_{i in C_k} || x_i - bar{x_k}||^2 end{eqnarray} 以下のようにランダムなクラスタを初期値として設定し、 より良いクラスタに更新することで計算量を下げる。 データ数(n)、クラスタ数(k)だと(mathcal{O}(nk))。 各データポイントにランダムなクラスタを割り当てる 以下の処理を繰り返す 各データポイントに割り当てられたクラスタのクラスタ中心を計算する 各データポイントについてクラスタ中心との距離を計算する 各データポイントについて距離が最小のクラスタを設定する k-meansで実際にクラスタリングするコードを書いてみようと思いましたが、 そろそろ適当なデータを見つけてくるのではなく、シミュレーション用のデータを作るところも やってみようと思います。次回、pythonでコード書きます。

[mathjax] タイトルの通り。Losso回帰と違って損失関数を偏微分するだけで出来そうなのでやってみる。 Ridge回帰は線形回帰の1種だけれども、損失関数として最小二乗法をそのまま使わず、 (L_2)ノルムの制約を付けたものを使う((L_2)正則化)。 データとモデル 教師データ(boldsymbol{y})、訓練データ(boldsymbol{x})があるとする。 (または目的変数(boldsymbol{y})、説明変数(boldsymbol{x})があるとする。) 例えば(p)次の属性データが(n)個あり、それらと結果の対応が分かっている状況。 begin{eqnarray} boldsymbol{y} &=& begin{pmatrix} y_1 \\ y_2 \\ vdots \\ y_p end{pmatrix} , boldsymbol{x} &=& begin{pmatrix} x_{11} & x_{21} & cdots & x_{n1} \\ x_{12} & x_{22} & cdots & x_{n2} \\ vdots & vdots & ddots & vdots \\ x_{1p} & x_{2p} & cdots & x_{np} end{pmatrix} end{eqnarray} モデルは以下。特徴ベクトル(boldsymbol{w})は訓練データの重み。 特徴空間において損失を最小化する特徴ベクトルを求める問題。 begin{eqnarray} boldsymbol{y} &=& boldsymbol{w} boldsymbol{x} + k \\ boldsymbol{w} &=& begin{pmatrix} w_1 & w_2& cdots &w_p end{pmatrix} end{eqnarray} 損失関数 普通の2乗損失に正則化項((L_2)ノルムを定数倍した値)を付けたものを損失関数として利用する。 正則化項の係数はハイパーパラメータとして調整する値。逆数なのはsklearnに従う。 begin{eqnarray} L(boldsymbol{w}) = |boldsymbol{y} - boldsymbol{w} boldsymbol{x}|^2 +C |boldsymbol{w}|^2 end{eqnarray} 特徴ベクトルは以下。(mathjaxでargminが出せない...) begin{eqnarray} newcommand{argmin}[1]{underset{#1}{operatorname{arg},operatorname{min}};} boldsymbol{w} = argmin w L(boldsymbol{w}) = argmin w |boldsymbol{y} - boldsymbol{w} boldsymbol{x}|^2 + C |boldsymbol{w}|^2 end{eqnarray} 特徴ベクトルを求める 勾配=0と置けば上の式の解を得られる。 損失関数が微分可能だからできる技。 begin{eqnarray} frac{partial L(boldsymbol{w})}{partial boldsymbol{w}} &=& 2 boldsymbol{w}^T (boldsymbol{y} - boldsymbol{w} boldsymbol{x}) + C boldsymbol{w} \\ &=& 0 end{eqnarray} 変形する。 begin{eqnarray} 2 boldsymbol{x}^T (boldsymbol{x}boldsymbol{w}-boldsymbol{y}) + C boldsymbol{w} &=& 0 \\ boldsymbol{x}^T (boldsymbol{x}boldsymbol{w}-boldsymbol{y}) + C boldsymbol{w} &=& 0 \\ boldsymbol{x}^T boldsymbol{x} boldsymbol{w} -boldsymbol{x}^T boldsymbol{y} + Cboldsymbol{w} &=& 0 \\ (boldsymbol{x}^T boldsymbol{x} +C E) boldsymbol{w} &=& boldsymbol{x}^T boldsymbol{y} \\ boldsymbol{w} &=& (boldsymbol{x}^T boldsymbol{x} + C E)^{-1} boldsymbol{x}^T boldsymbol{y} end{eqnarray} テストデータを作る 練習用にsklearnのbostonデータを使ってみる。 ボストンの住宅価格が目的変数、属性データが説明変数として入ってる。 import pandas as pd import numpy as np from pandas import Series,DataFrame import matplotlib.pyplot as plt from sklearn.datasets import load_boston boston = load_boston() boston_df = DataFrame(boston.data) boston_df.columns = boston.feature_names print(boston_df.head()) boston_df[\"PRICE\"] = DataFrame(boston.target) # CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX PTRATIO B LSTAT PRICE # 0 0.00632 18.0 2.31 0.0 0.538 6.575 65.2 4.0900 1.0 296.0 15.3 396.90 4.98 24.0 # 1 0.02731 0.0 7.07 0.0 0.469 6.421 78.9 4.9671 2.0 242.0 17.8 396.90 9.14 21.6 # 2 0.02729 0.0 7.07 0.0 0.469 7.185 61.1 4.9671 2.0 242.0 17.8 392.83 4.03 34.7 # 3 0.03237 0.0 2.18 0.0 0.458 6.998 45.8 6.0622 3.0 222.0 18.7 394.63 2.94 33.4 # 4 0.06905 0.0 2.18 0.0 0.458 7.147 54.2 6.0622 3.0 222.0 18.7 396.90 5.33 36.2 散布図行列を表示してみる。 PRICEと関係がありそうなZN,RM,AGE,DIS,LSTATの5個を使ってみる。 pg = sns.pairplot(boston_df) plt.show() pg.savefig(\'boston_fig.png\') 特徴ベクトルを自力で計算する これを自力で計算してみる。(C=0.01)、(C=0)、(C=100)としてみた。 begin{eqnarray} boldsymbol{w} &=& (boldsymbol{x}^T boldsymbol{x} + C E)^{-1} boldsymbol{x}^T boldsymbol{y} end{eqnarray} X_df = boston_df.drop(columns=[\'CRIM\',\'INDUS\',\'CHAS\',\'NOX\',\'RAD\',\'TAX\',\'PTRATIO\',\'B\',\'PRICE\']) X = X_df.values y = boston.target.T C1 = 0.01 C2 = 0 C3 = 100 e = np.identity(5) w1 = np.dot( np.linalg.inv(np.dot(X.T , X) + C1 * e), np.dot(X.T,y)) w2 = np.dot( np.linalg.inv(np.dot(X.T , X) + C2 * e), np.dot(X.T,y)) w3 = np.dot( np.linalg.inv(np.dot(X.T , X) + C3 * e), np.dot(X.T,y)) print(w1) # [ 0.05338557 5.40396159 -0.01209002 -0.83723303 -0.63725397] print(w2) # [ 0.05338539 5.40403743 -0.01209427 -0.83728837 -0.63725093] print(w3) # [ 0.05612977 4.76664789 0.02374402 -0.38576708 -0.66137596] (C=0)のとき、つまり最小二乗法のとき。 sklearnを使う sklearnのridge回帰モデルを使うと以下みたいになる。 from sklearn.linear_model import Ridge from sklearn.model_selection import train_test_split Xf_train,Xf_test,yf_train,yf_test = train_test_split(X,y,random_state=0) ridge = Ridge().fit(Xf_train,yf_train) print(f\"accuracy for training data:{ridge.score(Xf_train,yf_train):.2}\") print(f\"accuracy for test data:{ridge.score(Xf_test,yf_test):.2f}\") # accuracy for training data:0.68 # accuracy for test data:0.58 print(ridge.coef_) # [ 0.06350701 4.3073956 -0.02283312 -1.06820241 -0.73188192] 出てきた特徴ベクトルを並べてみる 自力で計算したものとsklearnに計算してもらったものを並べてみる。 似てるのか似ていないのかよくわからない .. けど、RMの寄与度が高いというのは似ている。 # 自力で計算 (C=100) # [ 0.05612977 4.76664789 0.02374402 -0.38576708 -0.66137596] # sklearnで計算 # [ 0.06350701 4.3073956 -0.02283312 -1.06820241 -0.73188192] 自力で計算したモデルの正答率を求めてみないとなんとも... そして、正規化項の係数の大小がどう影響するのか、あまり良くわからなかった..。 (L_2)ノルムの制約を付けると、パラメタの大小が滑らかになると言いたかったのだけども。 あと、訓練データに対して68%、テストデータに対して58%という感じで、 大して成績が良くない...。　

[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の時。むぅ. 真ん中の塊の中を左下から右上に向かってギザギザに進んで欲しいのだが、 物凄い引っ張られよう。異常値が原因というよりは、引っ張られ過ぎなんだな。