ホーム » 「線形回帰」タグがついた投稿

タグアーカイブ: 線形回帰

GLM（一般化線形回帰）をPythonでやる（Poisson）

2019年5月12日 / コメントする

こんにちは。

今回はたまに聞くであろうGLM、すなわち、一般化線形回帰についてです。回帰といえば今まで線形回帰とかちょろっとやりました。せっかくなので回帰についてちょっとだけ復習してから本題に入りましょう。

上のデータを回帰することを考えます。今まで線形回帰ばかりしてきました。しかし、今回のデータは線形っぽくない。このような時はどうすればいいのか？一般化線形モデルの登場です。

データが線形ではないとき
誤差が $\epsilon \sim \mathcal{N}(0,\sigma)$ でない
$y$ がカテゴリカルとかだったら

上記のような場合、線形回帰よりも一般化線形回帰が有効です。一般化線形モデルに入る前に単回帰、重回帰についてもついでに復習しておきましょう。

まあ、簡単ですね。では一般化線形モデルへ。

GLMでは上記の３つのファクターが重要になります。これらによりモデルが決まります。

線形予測子
リンク関数
確率分布

それぞれ、

説明変数の線形モデル
線形モデルに移す関数
目的変数の従う分布

になります。具体的に例を見ていきましょう。

普通の線形モデルは上のようになります。これだけだとわからないのでポアソン分布を例にとってみましょう。

目的変数が非負ということが１つの特徴でしょう。

$\beta_0 + \beta_1 x$ を $exp$ に食わせることで指数関数を用いてモデルを構成しています。線形予測子を得るためには $exp$ を外せばいいので $\log$ つまり対数がリンク関数になることがわかります。言い換えれば、 $y_i = \beta_0 + \beta_1 x$ と単純にモデリングしてしまうと目的変数が正規分布に従わないせいで誤差が大きくなってしまいます。それを避けるためにリンク関数により $\log \lambda_i = \beta_0 + \beta_1 x$ と変形していわば補正をかけます。（ $y_i$ が平均 $\lambda_i$ に従うポアソン分布から生成されたと仮定します）

これを用いようとするとPythonよりもRがメジャーな選択肢になりますが、僕はPython派なのでこっちを使います。以下がコードになります。

北大のレジュメの問題はRで解かれておりこれをPythonで解いた結果近い値が得られたのでうまく回帰できたと思います。しかしもう世の中はニューラルネットワークなのか、、、、

でわ

Ridgeを数式から実装まで(理論)

2019年2月10日 / コメントする

こんにちは。
料理にはまっているkzです。

本日はリッジ回帰行きましょう。次回実装します。🐷

Ridge Regression

過学習を防ごうということで導入された正則化回帰の一例。
– $L_2$ 正則化項を使用した回帰model
– パラメータに上限あり
たまにリッジ回帰はペナルティーありの線形回帰という記事を見かけますがそうではなく多項式回帰でもあります。つまりモデルが

(1) $\begin{align*} y = \beta_0 + \beta_1x_1 + \beta_2x_2 \cdots + \beta_ix_i + \epsilon \end{align*}$

だけでなく

(2) $\begin{align*} y = \beta_0 + \beta_1x + \beta_2x^2 \cdots + \beta_ix^i + \epsilon \end{align*}$

もOKということ。

Notation

Number of data: $N$
Dimension: $n$
Parameters: $b$

つまり各データは

(3) $\begin{align*} x_i := x^{(i)} = (x_1, x_2, \cdots, x_n)^T \end{align*}$

と表されます。 $b$ も同様。

Constrained minimizationは

(4) $\begin{align*} \texttt{minimize } \sum_{i=1}^N (y_i - x_i^Tb)^2 \end{align*}$

(5) $\begin{align*} \texttt{subject to } \sum_{i=1}^n b_i^2 \leq t ~~ , ~~ 0\leq t \end{align*}$

となりますが扱いずらいのでベクトル表示に変更し同時にKKT条件より最適化問題を同値なものに書き換えます。

(6) $\begin{align*} \texttt{minimize } |(y-Xb)|_2^2 + \lambda |b|_2^2 ~~ , ~~ 0\leq \lambda \end{align*}$

ではまず解を求めましょう。いつも通り微分をゼロとおいて解くだけ。今回のコスト関数は次の通り

(7) $\begin{align*} g: b \rightarrow |(y-Xb)|_2^2 + \lambda |b|_2^2 \end{align*}$

として $b$ で微分する。
– 「ベクトルで微分・行列で微分」公式まとめ

(8) $\begin{align*} \frac{\partial g}{\partial b} = -2X^Ty + 2bX^TX+2\lambda \end{align*}$

したがって

$b_{ridge} = (X^TX + \lambda I)^{-1}X^Ty$

$\lambda$ は正則をキープするためのものでしたね。

ここまでは前回の内容をちょっとだけ丁寧にやっただけのおさらいのようなものです。ここからが数学チック。

さて、この解が最適解かどうかを関数の凸性から確認してみましょう。凸性は前回やりましたがお椀型なら極小値が最小値になることでしたね。
– 第2章凸関数
– Convex functions

(9) $\begin{align*} f \in \mathcal{C}^\infty: \mathbb{R}^n \rightarrow \mathbb{R} \end{align*}$

に対して次は同値である。
– $f$ が凸である
– $\forall x\in \mathbb{R}^n$ に対してヘッセ行列 $\nabla^2 f(x)$ がpositive-semi-define

では

(10) $\begin{align*} \frac{\partial^2 g}{\partial b^2} = 2X^TX + 2I_n \end{align*}$

となりますね。ここでこれがpositive-semi-defineであることは次からわかります。 $\forall y\in \mathbb{R}^n$ に対して

(11) $\begin{align*} y^T X^TX y = |Xy|_2^2 \geq 0 \end{align*}$

とノルムの定義よりわかりました。

$g$ はほんまに凸か？

今見たのは

(12) $\begin{align*} H: b \rightarrow |(y-Xb)|_2^2 \end{align*}$

の凸性ですよね。そして

(13) $\begin{align*} K: b \rightarrow \lambda |b|_2^2 \end{align*}$

こちらは明らかに凸。では $g = H + K$ が本当に凸かどうか。これは直感的には明らかですが証明が簡単なのでやってみましょう。

集合 $X$ を定義域とする関数 $H,K$ がともに凸関数であると仮定する。
このとき $\forall x,y \in X ~~ , ~~ \forall \lambda \in [0,1]$ に対して $g$ が凸であることを示す。

(14) $\begin{align*} (H+K)(\lambda x + (1-\lambda)y) = H(\lambda x + (1-\lambda)y) + K(\lambda x + (1-\lambda)y) \end{align*}$

(15) $\begin{align*} \leq \lambda H(x) + (1-\lambda)H(y) + \lambda K(x) + (1-\lambda)K(y) \end{align*}$

(16) $\begin{align*} = \lambda (H+K)(x) + (1-\lambda)(H+K)(y) \end{align*}$

Q.E.D.

Lassoは結構めんどくさいんですよね。。
でわ

線形回帰って何？

2018年11月12日 / コメントする

こんにちは、本日はとりあえず実装をする前に、パパッとデータの確認です。

$y = 3x - 4 + noise$

で二つ作りました。見ての通り上の方がノイズが大きいです。今回の予測モデルは直線

$y=mx+b$

とします。サクッとPYTHONで書くと下のような結果がとれます。まずノイズが小さい方のデータに対して

コストが減っており、それに合わせてパラメータ $(m,b)$ も本来の $(3,-4)$ に収束しようとしています。上のプロットが見にくいのですが、学習の際のパラメータ履歴で直線をプロットしたものです。次にノイズが大きい方の結果を見てみましょう。

同様にコストも減り、パラメータも本来の値に近ずいてはいますがどおも届きそうにないです。iteration=150あたりで学習が止まっているのでこれ以上は赤点に近づかないでしょう。これはノイズが大きいことが原因です。

とりあえず線形回帰を二変数で行いましたが多変数になってもやり方は全く変わりません。コードを書くのが面倒になりますが　sklearnにはすでにlinear regression含め他にも多くの手法があるので簡単に回帰ができます。詳しくは「python 線形回帰 sklearn」などでググればいっぱい出てくるので見てください。

ところで

線形回帰をいい例題にして行列の勉強を簡単にしましょう。というのは何かというと、実は、勾配法使わなくても一発で解が出るんです。

$y_i = 1 * \beta_0 + x_{i1} * \beta_1 + x_{i2} * \beta_2 + \cdots + x_{ik} * \beta_k + \epsilon_i$

というモデルを考えます。これを行列表記すると

$\begin{bmatrix} y_1 \ y_2 \ \vdots \ y_n \ \end{bmatrix} = \begin{bmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1k} \ 1 & x_{21} & x_{22} & \cdots & x_{2k} \ \vdots & \vdots & \vdots & \vdots & \vdots \ 1 & x_{n1} & x_{n2} & \cdots & x_{nk} \ \end{bmatrix} \begin{bmatrix} \beta_0 \ \beta_2 \ \vdots \ \beta_k \ \end{bmatrix} + \begin{bmatrix} \epsilon_1 \ \epsilon_2 \ \vdots \ \epsilon_n \ \end{bmatrix}$

$y = X \beta + \epsilon$

ただし、 $n$ はデータ数、 $k$ は次元（変数の数）とした。ここで前回の理論編を思い出すと今回の最適化は誤差二乗和の最小化であった。よってノルムの定義を思い出せば

$\frac{1}{2} | \epsilon |^2 = \frac{1}{2} | y - X \beta |^2$

を最小化すればいい(1/2は計算上の都合、別になくても良い)、これを計算する。ベクトルの微分については

ここあたりがいい教材だと思う。もちろん「ベクトル　微分」で検索しても多く出てくる。

今回使うのは次の3つ。

$\frac{\partial }{\partial x} a^Tx = a$

$\frac{\partial }{\partial x} x^Ta = a$

$\frac{\partial }{\partial x} x^TAx = (A + A^T)x$

では計算すると

$\frac{1}{2} | y - X \beta |^2 = \frac{1}{2} (y - X \beta)^T (y - X \beta)$

$\frac{1}{2} | y - X \beta |^2 = \frac{1}{2} ( \beta^TX^TX \beta - 2y^TX\beta +y^Ty )$

$\nabla_\beta \frac{1}{2} | y - X \beta |^2 = (X^T X + (X^T X)^T)\beta - 2X^Ty$

これを0と置くと

$X^T X \beta = X^Ty$

ここで

$X^T X$ が正則行列を仮定すると

(機械学習の観点から「正則行列」を説明すると「どの列もほとんど似通ってない」になる)

$\beta = (X^T X)^{-1} X^Ty$

と求まる。これが最も誤差が最小であるときのパラメータである。

なんの射影？どこへの射影？

そう、この疑問は大切だ。今私たちが扱っているのは $\frac{1}{2} | y - X \beta |^2$ である。 $X \beta$ についてじーっと見るとこれは $X$ の各列ベクトルの線型結合(四則演算)となることがわかる。たとえば $X = [X_1 X_2 \cdots X_n]$ と表し、 $\beta=[ \beta_1 \beta_2 \cdots \beta_n ]^T$ とすると

$X \beta = X_1 \beta_1 + X_2 \beta_2 + \cdots X_n \beta_n +$

となることがわかる。ここで各 $X_i$ はベクトルであり $\beta_i$ はスカラーであることに注意する。したがって $X$ の各列ベクトルの線型結合で出来上がったベクトルの中でもっとも $y$ に近いやつ(似ているやつ)が

$X\beta = X(X^T X)^{-1} X^Ty$

となるのだ。ここで

$X(X^T X)^{-1} X^T$

を

$X$ の列ベクトルが張る空間への射影行列

と呼ぶ。「射影」という考え方は機械学習において非常に重要なので、ぜひ覚えておきましょう。

ねぇPython、線形回帰って何？（理論編）

2018年11月8日 / コメントする

こんにちは。

勾配法を終えました。差分法も終えました。目的関数も終えました。いい流れですねえ。余談ですが最近、友人に機械学習の入門書を選んでくれと言われ書店へ行ったのですがブームってすごいですね。今やいい本がたくさんあります。コードも丁寧に説明してあって数式も導出してあって、、、

今回のデータセット

では勾配法を使って回帰します。まずはモデルを定義します。今回はこのデータにいい感じに刺さる直線が欲しいですよね。よって今回はモデルを次のように定義します。

$\hat{y} = m x + b$

はい、次は目的関数を定義します。再度状況確認しましょう。

みての通りどんなにいい線を引いても「誤差(赤い線)」がありますよね。仕方ないです。現実のデータには誤差があります。これをよくNOISEと言います。とはいえ、各データの誤差の和（赤い線の和）が最小の時が一番いい線です。よって

$error = y - \hat{y} = y - ( mx + b )$

を最小化したいですよね。けどこのままだと扱いずらいので損失関数を $error = |y - ( mx + b )|$ とするとどうでしょうか。微分できないのであんまり良くないです。なのでこうしましょう。

$error = (y - ( mx + b ))^2$

微分可能であり、さらに下に凸の凸関数ですね。凸関数って何がいいんでしたっけ？そうです、お椀型なので勾配法をかけると大域最適解がゲットできます。よって今回のコスト関数はこれの和になるので

$E = \frac{1}{N} \sum_{I=1}^{N} (y_i - ( mx_i + b ))^2$

ではそれぞれ勾配を求めましょう。

$\frac{\partial E}{\partial m} = \frac{2}{N} \sum_{I=1}^{N} (-2x_i) (y_i - ( mx_i + b ))$

$\frac{\partial E}{\partial b} = \frac{2}{N} \sum_{I=1}^{N} (-2) (y_i - ( mx_i + b ))$

今回は２次元のデータだったので簡単でした、変数が増えてもやることは同じです。ちなみにこれは線形ですが曲げたければモデルを２次式・３次式、、、にすればいいですし、周期性があれば三角関数にしてもいいですし、そこは自由です。ただ現実のデータは可視化できないのが多いと思うので今回のように予想は立てられないと思います。

データが可視化できない？

いえいえ、機械学習の力を使えばできるんです。その名も

PCA

主成分分析とも言います。高次元空間のデータを低次元空間へピョンと飛ばすテクニックです。例えば500次元のデータを２次元に変換するとx-y平面で可視化できますね。便利便利。

おっと話が逸れてしまいました。ところで今回の実装テーマの「線形回帰」についてですがそもそも「回帰」ってなんでしょう？

いろんな捉え方がありますが、今回の線形回帰で得られるのはいい直線でしたよね。言い換えると「xとy」をうまく結びつけてるヤツですよね？なので新しく $(x')$ というデータを得た時このモデルに適応させると $(y')$ が得られますよね。逆も同じです。
他には　
$y = 12x_1^2 + 8x_2 - 62$

というモデルだとわかった時に、目的変数 $(y)$ の値が説明変数 $(x_2)$ に比べて $(x_1)$ に強く依存することがわかります。つまり目的変数に対する寄与度が相対的にわかるんです。ちなみにこれも手法があってLasso Regression, Ridge Regression, ElasticNet Regressionなどあります。
あとはこれは僕の中では回帰っぽくないのですがロジスティック回帰です。これは分類に使われます。一見回帰はデータをうまく表現するモデルを得ることと思っちゃうんですが分類にも使えるんですね。