「はじめての機械学習」一覧

[mathjax] 参考書籍を読んでいく. 今回はロジスティック回帰. 未知のサンプルのクラスの所属確率を見積もることができる, という重要なことが書かれている. また重みの学習は勾配降下法のフレームワークの中で統一されていて一連の理解が進む. 一度は理解しておいた方が良い内容だと思った. パーセプトロンとADALINE 線形和(z=w^T)を考えたとき 「サンプル(x)が決定境界の上にある」 というのは(w^Tx=0). (w^T x > 0), (w^T x < 0) を決定境界の 「こちら側」 か 「向こう側」 か対応させる. パーセプトロンのアイデアは以下の通り. 線形和 (z=w^T x)を決定境界とし, (phi(z)in (-1,+1))を予測ラベルとする 予測が誤った回数, つまり真のラベルと予測ラベルの差の合計がゼロになるように(w)を学習する (phi(z))はステップ関数で連続ではないので解析的な式の変形ができない 決定境界により完全に分離できなければ学習が収束することはない 決定境界から離れている度合いを連続で凸な関数として定義し, その関数が極小となる(w)を求めていこう, というアイデアに拡張したのがADALINE. (phi(z) = z). コスト関数(J(w) = frac{1}{2} sum_i left( hat{y^{(i)} - phi ( z^{(i)})} right)^2 ) トレーニングセット全体から(J(w))の極小化を考えることが可能. (勾配降下法) トレーニングセットの一部を使って(J(w))の極小化を考え, 他のトレーニングセットで反復することも可能. (確率的勾配降下法) ロジスティック回帰 決定境界を表す線形和(z=w^T x)について(phi(z) in (mathbb{R}| 0 leq z leq +1))を考える. 線形和(z)がクラス1に分類される確率を表す(phi(z))を考える. (phi(z))をどう作るのか 線形和を(0)から(1)の実数全体に変換する関数を考える. オッズ比の対数をとったロジット関数(logit(p))を使う. ここで (p)は(x)が与えられたときにサンプルデータがクラス1に属する確率. begin{eqnarray} logit(p(y=1 | x)) &=& log frac{p}{1-p} \\ &=& w_0 x_0 + w_w x_w + cdots + w_m x_m \\ &=& sum_{i=0}^m w_i x_i \\ &=& w^T x end{eqnarray} つまり, 線形和(z = w^T x in (mathbb{R} | 0 leq z leq 1) ) (p)は未知であるということが重要. この式から(p)を予測していく. (f(p) = log frac{p}{1-p}) の逆関数が分かれば良い. begin{eqnarray} f(p) &=& log frac{p}{1-p} \\ e^{f(p)} &=& frac{p}{1-p} \\ (1-p)e^{f(p)} &=& p \\ e^{f(p)} &=& p(1+e^{f(p)}) \\ p &=& frac{e^{f(p)}}{1+e^{f(p)}} \\ p &=& frac{1}{1+e^{-f(p)}} end{eqnarray} (f(p)=z)としていたので以下のようになる. (zとpが混在していてたぶん数式としての表現は間違ってる.. ) begin{eqnarray} p &=& frac{1}{1+e^{-z}} end{eqnarray} (p)は確率なので(p in (mathbb{R}|0 leq p leq 1)). (phi(z) in (mathbb{R}|0 leq p leq 1) )とすると, begin{eqnarray} phi(z) = frac{1}{1+e^{-z}} end{eqnarray} (z-phi(z))をGeoGebraでプロットすると以下のような感じ. (z)を大きくしていくと(phi(z))は限りなく(1)に近づく. (z)を小さくしていくと(phi(z))は限りなく(0)に近づく. (phi(z)=)の出力は正事象が起こる確率であると解釈される. (phi(z)=0.8)である場合, 80%の確率で正事象が起こり,20%の確率で正事象が起こらないことを表す. トレーニングデータからパラメタを学習し, ある(phi(z))が得られたとして, 未知のサンプルのクラスの所属関係を確率で見積もることができる. 天気予報で晴れの確率70%、晴れでない確率30%のように. 予測の際に,クラスの所属関係の確率を見積もることができるのが神がかるレベルで重要!! もし(phi(z))を2値分類に当てはめるのであれば, (phi(z)=0.5)を境に上半分を(1)、下半分を(0)と予測するものとする. begin{eqnarray} hat{y} = begin{cases} 1 & phi(z) ge 0.5 \\ 0 & phi(z) lt 0.5 end{cases} end{eqnarray} それって, (z=0)を境に左半分, 右半分に割り付けるのと同じなので, begin{eqnarray} hat{y} = begin{cases} 1 & z ge 0 \\ 0 & z lt 0 end{cases} end{eqnarray} ロジスティック回帰の重みの学習 では(w)をどうやって学習するか. 式をイジるだけで面倒なだけで,これ以上は洞察得られないなと思うので, 流れだけ書いてみる. 学生時代に形態素解析っぽい何かのアルゴリズムを考えたことがあって, 形態素どうしのベターな接続を求め際に, 接続コストをコーパスから作った確率で表して, 尤度を最大化する問題として解いた記憶があって, 今更合点が行くという謎. 勉強不足で尤度最大化のアイデアがあってそうした訳ではないので, もっと勉強していたらよかったな..と. ADALINEの説明の中で, 真の値と予測値の差の2乗誤差をコスト関数として, 連続で凸な関数の極小化の問題として(w)を求めていった. ロジスティック回帰の場合, 予測は(phi(z))も連続で凸な関数なので, 同じロジックで(w)を求めることもできる. コスト関数(J(w))は同様に以下. begin{eqnarray} J(w) &=& sum_i frac{1}{2} left( phi(z^{(i)}) - y^{(i)} right)^2 end{eqnarray} (phi(z)=frac{1}{1+e^{-z}})だと, これを解析的に微分するのは不可能. トレーニングデータが(m)個あった場合に(i)番目のトレーニングデータにおける確率は以下. (x)も(w)もパラメタですよ,って書き方. begin{eqnarray} P(y|x;w)=P(y^{(i)}|x^{(i)};w) end{eqnarray} 求めたいのは, 全てのトレーニングデータについて(P(y))の積が最大になる(w). 積,つまり尤度を書き下すと以下のようになる. begin{eqnarray} L(w) &=& prod_{i=1}^m left( P(y^{(i)}) | x^{(i)};w right) \\ &=& prod_{i=1}^m left( phi(z^{(i)}) right)^{y^{(i)}} + left( 1-y^{(i)}right) log left( 1-phi(z^{(i)}) right) end{eqnarray} 尤度そのものの最大化するよりも対数をとったものを最大化する方が楽になる. (対数をとると指数の肩が定数倍になり積が和になる. 相当簡単になる. ) begin{eqnarray} l(w) &=& log L(w) \\ &=& sum_{i=1}^{n} left( y^{(i)} log( phi(z^{(i)}) ) + (1-y^{(i)}) log (1-phi(z^{(i)})) right) end{eqnarray} これを最大化するということは, これに(-1)を掛けたものを最小化するということ. そもそもコスト関数(J(w))を立てて最小化しようとしていたので, これを(J(w))と置いてしまう. begin{eqnarray} J(w) = sum_{i=1}^{n} left( -y^{(i)} log( phi(z^{(i)}) ) - (1-y^{(i)}) log (1-phi(z^{(i)})) right) end{eqnarray} ADALINEの勾配降下法で, コスト関数を最小にするように(w)を更新していった. つまり, (w := w - Delta w)という更新をしていった. ロジスティック回帰において, コスト関数を対数尤度で書き直していて, 対数尤度を最大化するように(w)を更新する. つまり, (w := w + Delta w)という更新をおこなう. 対数尤度を最大化する(w)を更新する話は, コスト関数を最小化する(w)を更新する話と同じ, なので, 結局のところ勾配降下法と同じ式となる. begin{eqnarray} Delta w_j &:=& - eta frac{partial J}{partial w_j} end{eqnarray} 尤度関数の偏導関数はシグモイド関数の偏導関数を使って以下のようになるらしい. (省略) シグモイド関数の偏導関数は以下 begin{eqnarray} frac{partial phi(z)}{partial z} &=& frac{1}{partial z} frac{1}{1-e^{-z}} \\ &=& frac{1}{1+e^{-1}} left( 1- frac{1}{1+e^{-z}} right) \\ &=& phi(z) (1-phi(z)) end{eqnarray} 対数尤度の偏導関数はシグモイド関数の偏導関数を使って以下. begin{eqnarray} frac{partial l(w)}{partial w_j} &=& left( y frac{1}{phi(z)} - (1-y) frac{1}{1-phi(z)} right) frac{partial}{partial w_j} phi(z) \\ &=& left( yfrac{1}{phi(z)} - (1-y) frac{1}{1-phi(z)} right) phi(z)(1-phi(z)) frac{partial}{partial w_j} z \\ &=& left( y(1-phi(z))-(1-y)phi(z) right) x_j \\ &=& left( y-phi(z) right) x_i end{eqnarray} コスト関数の最小化と対数尤度関数の最大化が同じ,というロジックが効いて, ADALINEの更新と全く同じになる. begin{eqnarray} w_j &:=& - eta frac{partial J}{partial w_j} \\ &=& eta sum_{i=1}^n left( y^{(i)} -phi(z^{(i)}) right) x_j^{(i)} end{eqnarray} 奇跡的..

[mathjax] 参考にした書籍でこの順序で誘導されていて理解しやすかったです. パーセプトロンによる一番簡単な教師あり学習を理解する ADALINEにより学習を凸で連続なコスト関数の最小化問題として捉える パーセプトロンの学習 2値のラベル((+1),(-1))が付与されているサンプルデータが与えられたとする. ラベル値が(-1)であるデータと, ラベル値が(+1)であるデータに分類できる. それらの境界が線形和で表される問題である場合, その境界を求めることができれば, 未知のデータに対してその境界の(+1)側か(-1)側かを判定できる. (m)次のサンプルデータと重みの線形和を考える. [ boldsymbol{z} = boldsymbol{w}^T boldsymbol{x} = w_1 x_1 + w_2 x_2 + cdots x_m x_m ] (z)に関するステップ関数を用意する. こうすることで (z)-(phi(z)) 空間において(z=theta)を境に線形和(z)が(-1),(+1)に分かれることを表現できる. [ varphi(z) = begin{cases} 1 & (z gt theta) \\ -1 & (z leq theta) end{cases} ] (theta)を左辺に移動し, (x_0 = 1), (w_0 = -theta) とすると, 線形和に(theta)の項を組み入れることができる. 新しい線形和(z\')が(0)になる箇所を境にステップ関数の応答値が切り替わる. [ varphi(z) = begin{cases} 1 & (z-theta gt 0) \\ -1 & (z-theta leq 0) end{cases} ] [ z\' = z - theta = w_0 x_0 + w_1 x_1 + w_2 x_2 + cdots x_m x_m ] あるパラメータ(boldsymbol{w})が存在するとして, データ(boldsymbol{x})について, (boldsymbol{w}^T boldsymbol{x}=0)を境界にステップ関数の応答値が切り替わる. 未知のデータ(x)に対して, (boldsymbol{w}^T boldsymbol{x})はラベルの予測値(+1), (-1)を応答する. 予測値(varphi(w^Tx))が正解であるサンプルデータに含まれる(y^{i})と同じとなる(w)を探したい. それを探していく. 最初,(w)に初期値を設定し, 以下の手続きによって(w)を更新していく. 上付きの添字はサンプルデータ内の順序数を表している. (y^{(i)})は(i)番目のサンプルデータ. 右下の添字は該当サンプルデータの各次元を表している. (y^{(i)_m})は(y^{(i)})の(m)番目の要素. 現在の(w)を使って予測値を計算する. (hat{y}=w^T x) (w)を更新する. (Delta w_j := eta(y^{i}-hat{y})x^{i}_j) 1個のサンプルデータで1回(w)が更新される. そもそもデータセットを綺麗に線形和が表す決定境界で分離できないような状態だと、 何度やっても予測したクラスラベルと真のクラスラベルの差が埋まらない。 あっちを立てればこっちが立たず、みたいになる。 予測したクラスラベルと真のクラスラベルを比較するこの方法だと、 (Delta w_j := eta(y^{i}-hat{y})x^{i}_j) が良い値なのか悪い値なのか、繰り返さないとわからない。 予測と真の値の比較を凸で連続な関数で表すことができれば、その関数の凸の底に向かうやり方で、 「今より良い次の値」に更新する方法を解析的に求めることができて都合が良い。 真の値と予測した値の差が「損失」とか「コスト」とか呼んで、 「コスト関数」を最小化するパラメータが良いパラメータなんだろう、という話が進んでいく。 ADALINEの学習 線形和をステップ関数に変換する部分があった。 こうして(w^T x varphi(w^T x))空間で (w^T x)が(-1),(+1)に分かれる(w)を求めようとした. [ varphi(z) = begin{cases} 1 & (z-theta gt 0) \\ -1 & (z-theta leq 0) end{cases} ] そうではなく以下の恒等式(左辺と右辺が同じ)を使い,真の値と予測した値の差の求め方を変えると, 差が凸で連続な関数で表せるようになり (解析的に微分できるようになり), 一番小さいところは? という話がしやすくなる. [ varphi(z)=z ] 最小二乗法のときやったやつで、差の2乗を足し合わると符号がキャンセルされてよくて、差を式で表せる。 [ J(w) = sum_i left( y^{(i)} - varphi(w^T x^{(i)}) right)^2 ] こういうのを誤差平方和と言うようで(frac{1}{2})倍するらしい. [ J(w) =frac{1}{2} sum_i left( y^{(i)} - varphi(w^T x^{(i)}) right)^2 ] (J(w))がなんで連続で凸なのかは知ったことではないが、 連続で凸であれば(J\'(w)=0)を解くと極小となる(w)が求まるのは高校生のときにやった話. (J(w))の接線(多次元だと接っする面)の傾きが(J\'(w))。 ただ(w)には次数があって傾きは以下(それぞれの次数毎の極限). [ nabla J(w) = frac{partial J(w)}{partial w_j} ] (w J(w))空間で(nabla J(w))は(w)における傾きなので, (w)に(-nabla J(w) cdot 1)を足すと, (J(w))が極小になる箇所に近づく. どれだけ行けば良いかわからないので(-nabla J(w) cdot eta)を足すことにする. (w)を以下の通り更新する. begin{eqnarray} w &:=& w + nabla J(w) \\ &:=& w -nabla J(w) cdot eta end{eqnarray} ちなみに(nabla J(w))は式をこねくり回すと計算できる. 最終的に更新は以下の通りとなる. [ w := w - eta sum_i left( y^{(i)} -varphi left( z^{(i)} right) right) x_j^{(i)} ] パーセプトロンの更新とそっくりな形が出てきて感動する... あっちは(varphi(z))が不連続なステップ関数でこっちは連続な関数. 式をこねくり回す際に, (varphi(z) = z)の恒等式の関係を使っていない. 何か連続な関数であればこれが成り立つ. 形式上,例えば(varphi(z))を以下(シグモイド)とすることでロジスティック回帰 (本当か? 次回確認.). [ varphi (z) = frac{1}{1+e^{-z}} ]

multivariate_normalを使って6個のデータポイント近隣にデータ達を発生させる。 薄くクラスタを見つけられそうだけれど境界は曖昧で大分被っているという状況。 そんなデータ達に6-meansをかけてみたとき、どうクラスタが出来るのかという実験。 被っている部分がどう別のクラスタに入るのか確認するため。 import numpy as np import matplotlib.pyplot as plt import seaborn as sns mu = [[0,0], [20,20], [50,50], [40,30], [40,10], [20,40]] sigma = [ [[30,20],[20,50]], [[20,30],[10,20]], [[60,40],[20,20]], [[60,20],[20,60]] ,[[30,10],[10,30]],[[50,20],[20,50]] ] points = 100 clusteres = [] for index in range(len(mu)): cluster = np.random.multivariate_normal(mu[index], sigma[index], points) dig = np.full((points,1),index+1, dtype=int) cluster = np.hstack((cluster,dig)) clusteres = np.r_[clusteres,cluster] if len(clusteres) > 0 else cluster plt.scatter(x=clusteres[:,0], y=clusteres[:,1],c=clusteres[:,2]) おもむろにsklearnのKMeansを使ってみる。 n_clustersは6、max_iterを10に設定してみた。 クラスタ中心がわかるように重ねてみた。 from sklearn.cluster import KMeans kmeans_model = KMeans(n_clusters=6, init=\'random\',max_iter=10).fit(clusteres[:,:2]) labels = kmeans_model.labels_ centers = kmeans_model.cluster_centers_ centers fig = plt.figure() ax1 = fig.add_subplot(1,1,1) ax1.scatter(x=clusteres[:,0], y=clusteres[:,1],c=labels) ax2 = fig.add_subplot(1,1,1) ax2.scatter(x=centers[:,0], y=centers[:,1], alpha=0.5,s=600,c=\"pink\",linewidth=2,edgecolors=\"red\") 当然のごとく、微妙に被っていた部分は異なるクラスタに分類された。 今回は当たり前なことを確認して終了。 2次元データじゃつまらない...。

[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%という感じで、 大して成績が良くない...。　

線形サポートベクトル分類器で画像認識する流れを理解したので、 定着させるために記事にしてみます。 当然、モデルの数学的な理解がないとモデルを解釈することは不可能だし、 正しいハイパーパラメータを設定することも不可能なので、数学的な理解は不可欠。 NumPy、pandas、matplotlibに慣れないと、そこまで行くのに時間がかかります。 こちらはPythonプログラミングの領域なので、数こなして慣れる他ないです。 機械学習用のサンプル画像で有名なMNISTを使ってNumPy、pandasの練習。 手書き文字認識用の画像データを読み込んでみる。サイズは28x28。各々1byte。 MNISTの手書き文字認識画像の読み込み まず読み込んでみて、データの形を出力してみる。 X_trainは、要素が3個のTupleが返る。3次。 1番外が60000。28x28の2次のndarrayが60000個入っていると読む。 1枚目の画像データはX_train[0]によりアクセスできる。 import tensorflow as tf minst = tf.keras.datasets.mnist (X_train,y_train),(X_test,y_test) = mnist.load_data() print(X_train.shape) # (60000, 28, 28) y_trainは要素が1個のTupleが返る。1次。 1枚目から60000枚目までの画像が0から9のいずれに分類されたかが入っている。 y_train[0]が4なら、1枚目の画像が4に分類された、という意味。 print(y_train.shape) # (60000,) データセットの選択 X_train,y_train、X_test,y_testから、値が5または8のものだけのViewを取得する。 そのために、まず値が5または8のものだけのインデックスを取得する。 NumPyのwhereはndarrayのうち条件を満たす要素のインデックスを返す。 X_trainに入っている60000件の2d arrayのうち、 値が5または8のインデックス(0-59999)を取得するのは以下。 index_train = np.where((X_train==5)|(X_train==8)) print(index_train) # (array([ 0, 11, 17, ..., 59995, 59997, 59999]),) index_test = np.where((X_test==5)|(X_test==8)) print(index_test) # (array([ 8, 15, 23, ..., 9988, 9991, 9998]),) インデックスを使って絞り込む。 X_train,y_train = X_train[index_train],y_train[index_train] X_test,y_test = X_test[index_test],y_test[index_test] print(X_train.shape) # (11272, 28, 28) print(X_test.shape) # (1866, 28, 28) 前処理 0-255の間の値を0-1の間の値に変換する（正規化）。 28x28の画像(2darray)を1x784(1darray)に整形する（平坦化)。 X_train,X_test = X_train / 255.0, X_test / 255.0 X_train = X_train.reshape(X_train.shape[0], X_train.shape[1] * X_train.shape[2]) X_test = X_test.reshape(X_test.shape[0], X_test.shape[1] * X_test.shape[2]) ベストなハイパーパラメータの選択 線形サポートベクトル分類器を作成する。 from sklearn.svm import LinearSVC linsvc = LinearSVC(loss=\"squared_hinge\",penalty=\"l1\",dual=False) 線形サポートベクトル分類器のハイパーパラメータCの選択 逆正則化パラメータCをGridSearchCVで探す。MBP2013Laterで学習(fit)に5分くらいかかった。 GridSearchCVからはC=0.2がbestと返ってくる。 from sklearn.model_selection import GridSearchCV param_grid = {\"C\":[0.025,0.05,0.1,0.2,0.4]} model = GridSearchCV(estimator=linsvc, param_grid=param_grid,cv=5,scoring=\"accuracy\",return_train_score=True) model.fit(X_train,y_train) print(model.cv_results_[\"mean_train_score\"]) # array([0.96291693, 0.96775192, 0.97059085, 0.97340754, 0.97626859]) print(model.cv.results_[\"mean_test_score\"]) # array([0.95626331, 0.95990064, 0.96158623, 0.9625621 , 0.96105394]) print(model.best_params_) # {\'C\': 0.2} 学習、精度評価 C=0.2を使って新しく学習させる。 linsvc = LinearSVC(loss=\"squared_hinge\",penalty=\"l1\",dual=False,C=0.2) linsvc.fit(X_train,y_train) 訓練データ、テストデータに対して正答率を求める。 訓練データについて97.2%、テストデータについて96.2%。 過学習すると訓練データが高くテストデータが低くなる。 from sklearn.metrics import accuracy_score pred_train = linsvc_best.predict(X_train) acc = accuracy_score(y_true = y_train,y_pred = pred_train) print(acc) # 0.9723207948899929 pred_test = linsvc_best.predict(X_test) acc = accuracy_score(y_true = y_test,y_pred = pred_test) print(acc) # 0.9619506966773848 モデルの解釈可能性 [mathjax] 線形SVMの決定境界(f(x))の係数をヒートマップっぽく表示して、どの係数を重要視しているかを確認する。 基本的に真ん中に画像が集まっているので、28x28の隅は使わないのが正しそう。 正則化パラメータによって係数の大きさを制御しているため、正則化パラメータを変えると係数が変わる。 今回のは(L_1)正則化なので、係数が0のものが増える..らしい（..別途調べる..)。 (f(x) = w_0 + w_1 x_1 + w_2 x_2 + cdots w_{784} x_{784}) import matplotlib.pyplot as plt weights = linsvc_best.coef_ plt.imshow(weights.reshape(28,28)) plt.colorbar() plt.show()

[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乗から線形に切り替わるところは連続になってる。) 次回以降実際のデータとモデルを使ってやってみる。