カテゴリー: 回帰

[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] 前の記事で線形単回帰において訓練データから回帰係数を求める方法を書いてみた。 標本平均を使って母平均を推測する話とリンクさせることで、 回帰係数の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}