線形回帰モデル

以下のような m 個のデータセットが与えられたとする。(i=1,2,...,m)

(y_i,x_{1i},x_{2i},\ldots,x_{ni})

yを目的変数、xを説明変数と呼ぶ。
y を x で回帰するとは、以下の式に対し、w を求めることである。

y = w_0 + w_1 * x_1 + \cdots w_n * x_n

wの求め方は以下で示す平均二乗誤差 Q の値が最小になるように定める。

Q = \sum_{i=1}^{m} \{y_i - w_0 - w_1 * x_{1i} - \ldots - w_n * x_{ni}\} ^ 2

w0,w1,...,wn で Qを偏微分したものを0とすることで、以下の決定方程式を得る。ただし、上部のバーはデータセットの平均をあらわす。

\begin{pmatrix}
\overline{y} \\
\overline{y x_1} \\
\vdots \\
\overline{y x_n} \\
\end{pmatrix}

=

\begin{pmatrix}
1 & \overline{x_1} & \ldots & \overline{x_n} \\
\overline{x_1} & \overline{x_1^2} & \ldots & \overline{x_1 x_n}\\
\vdots & & \ddots & \vdots\\
\overline{x_n} & \overline{x_n x_1} & \ldots & \overline{x_n^2}
\end{pmatrix}

\begin{pmatrix}
w_0 \\
w_1 \\
\vdots \\
w_n \\
\end{pmatrix}

線形回帰はその名の通り、目的変数・説明変数間の線形関係を見るものであり、線形的ではない関係を見るのには不適切。

決定係数

データがどれだけフィットしているかを知るために、決定係数という概念を導入する。
まず、決定方程式により求めた w0,w1,...,wn とデータセットを用い、内挿値を以下で定義する。

\widehat{y_i} = w_0 + w_1 * x_{1i} + \ldots + w_n * x_{ni}

直感的にはこれが yi と近いほど、データがモデルにフィットしているといえる。(ある i　においてだけではだめで、全体的に近いほうが良い）
これを表す指標として、決定係数を定義する。

R^2 = \frac{\sum (\bar{y} - \widehat{y_i})^2}{\sum (y_i - \bar{y})^2}

分母を回帰変動、分子を全変動と呼ぶ。
また全変動 - 回帰変動を残差変動と呼び、これが上記でいう yi と内挿値の近さにまつわる量になる。
つまり、決定係数が1に近いほど回帰直線がデータによく当てはまっているといえ、逆に決定係数が0に近いほど回帰直線はデータに当てはまっていないといえる。

決定係数には、説明変数が増えると大きくなるという特徴がある。
実際にはデータがフィットしていなくても、適当な説明変数を増やすだけでフィットしているように見えてしまう。
これを解決する方法として、自由度修正決定係数を用いる方法があり、今回はその効果を実装演習で見てみることとする。

実装演習

プログラムは以下。
説明変数の数を増減させ、決定係数の値を見てみる。data の部分を増やしていく。

from sklearn.datasets import load_boston
from pandas import DataFrame
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
import numpy as np

boston = load_boston()
# 説明変数追加
df = DataFrame(data=boston.data,columns=boston.feature_names)
# 目的変数追加
df['PRICE'] = np.array(boston.target)

# 目的変数部分抜出
target = df.loc[:,'PRICE']

# モデル定義
model=LinearRegression()
# 線形回帰
model.fit(data,target)

# データ総数
m=506
# 説明変数
k = 2
# y の平均
ym = DataFrame.mean(target)
yma =  [ym for i in range(n)]
# 内挿値をだす
y_in = model.predict(data)
# 残差変動を出す
e = m * mean_squared_error(y_in, target)
# 全変動を出す
eall = m * mean_squared_error(target, yma)
# 決定係数をだす
R = 1 - (e / eall)
print(R)
# 自由度修正決定係数を出す
Rm = 1 - ( e / (m - k - 1)) / (eall / (m - 1))
print(Rm)

結果

説明変数の数	決定係数	自由度修正決定係数
1	0.15078046904975706	0.14909550966295104
2	0.23398843834155303	0.23094266672462094
3	0.29371357158118694	0.28949273635159245
5	0.33131270122233847	0.3246258282345619
10	0.6396628026756687	0.6323832633357832
13	0.7406426641094094	0.7337897263724629

考察

決定係数から、説明変数の数を増やすにつれてよりよいモデルになっていることがわかる。
13 個の説明変数すべてを使った場合の決定係数が 0.74 もあるのでかなりフィットしているといえる。
自由度修正決定係数については、思ったより決定係数と差がなかった。
説明変数の数をデータの数に比べて著しく大きくしない限りは決定係数の値を見れば十分そうであった。

非線形回帰モデル

基本的には線形回帰モデルと同じ。線形回帰モデルの拡張といってもいい。
以下の式で回帰する。

y_i = w_0 + w_1 * \phi_1 (x_{1i}) + w_2 * \phi(x_{2i}) + \cdots + w_n * \phi(x_{ni})

回帰関数は例えば以下のようなものがある。

\phi(x) = x \\
\phi(x) = x ^ j \\
\phi(x) = exp( \frac{(x-m)^2}{2h} )

上から線形回帰(多項式関数の一種)、多項式関数、ガウス型基底関数である。
これらを組み合わせて回帰する。
ガウス型基底関数はどうやら以下のように設定することが多いらしい。

\phi_i(x) = exp( \frac{(x-x_i)^2}{2h} )

なぜこのように設定するのかはよくわからないので、文献をあさる予定。
実装演習で一部この辺りも考察する。

決定方程式の形が個人的に好きなため、行列表現の形を示していなかったが、せっかくなのでこちらで示しておく。

X^t X w = X^t y

決定方程式の両辺m倍したものになっている。
一応それぞれの変数の中身を以下に示す。

X_i = (1,  \phi_1 (x_{1i}) ,  \phi_1 (x_{2i}), \cdots ,  \phi_1 (x_{ni})  )　\\
w^t = (w_0, w_1, \cdots , w_n) \\
y^t = (y_1, \cdots , y_m) \\

正則化法

過学習の対策。表現力を抑止することで、過学習を抑える。
Q に正則化項R 追加して回帰する。

Q + \alpha R(w)

R は以下の二つが例として挙げられることが多い。

R(w) = \sum_{i=1}^{n} | w_i | \\
R(w) = \sum_{i=1}^{n} w_i^2

上から、L1ノルム(Q+RをLasso推定量、L2ノルム(Ridge推定量という。
講義で聞いたときは

R(w) = \sum_{i=8}^{n} |w_i|

などとして、例えば高次項の影響を意図的に低減させるのかと思ったが、実際には全体的に表現力を落とすらしい。
高次項の影響が大きいモデルだった場合、自然にその影響が顕著に低減させられるといったほうが適切かもしれない。

実装演習

今回はモデルをわざと過学習にさせる -> 正則化法でそれを取り除けるかを確認する。ことを目的とする。

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.linear_model import LinearRegression

%matplotlib inline

#seaborn設定
sns.set()
#背景変更
sns.set_style("darkgrid", {'grid.linestyle': '--'})
#大きさ(スケール変更)
sns.set_context("paper")

n=100

def true_func(x):
    z = 1-48*x+218*x**2-315*x**3+145*x**4
    return z

def linear_func(x):
    z = x
    return z

# 真の関数からノイズを伴うデータを生成

# 真の関数からデータ生成
data = np.random.rand(n).astype(np.float32)
data = np.sort(data)
target = true_func(data)

# 　ノイズを加える
noise = 0.5 * np.random.randn(n)
target = target  + noise

# ノイズ付きデータを描画

plt.scatter(data, target)

plt.title('NonLinear Regression')
plt.legend(loc=2)

# ガウス型基底を用いるケース
from sklearn.metrics.pairwise import rbf_kernel

data = data.reshape(-1,1)
target = target.reshape(-1,1)
# ガウス型基底に変換 exp{(x_i - x_j)^2}
kx = rbf_kernel(X=data,Y=data,gamma=1)

from sklearn.linear_model import LinearRegression
model1 = LinearRegression()
model1.fit(kx, target)

y_in = model.predict(kx)

plt.scatter(data, target)
plt.plot(data, y_in)
print(model1.score(kx, target))

from sklearn.kernel_ridge import KernelRidge

model2 = KernelRidge(alpha=0, kernel="rbf")
model2.fit(data, target)

y_in = model2.predict(data)

plt.scatter(data, target)
plt.plot(data, y_in)
print(model2.score(data, target))


model3 = KernelRidge(alpha=0.0002, kernel='rbf')
model3.fit(data, target)

y_in = model3.predict(data)

plt.scatter(data, target)
plt.plot(data, y_in)
print(model3.score(data, target))

考察

model1

ガウス型基底のmを自分自身であてはめたモデル。
正則化項をつけていないので、もっとあからさまに過学習になるかと思ったが、
見た目はかなりいい感じにfitした。
mを自分自身に設定するのがかなりいいパラメータの決め方ということなのだろうか？

決定係数 = 0.8334929586559503

model2

KarnelRidge(rbf) に、正則化項なしであてはめたモデル。
このときのガウス型基底のmがいくつに設定されているのかが、ドキュメントからはわからなかった。
実装読むしかないのか・・・？
おそらく m の設定がめちゃめちゃfitするように設定されてしまい、以下のように過学習が起こると考えられる。

決定係数 = 0.9999999999999978

model3

正則化項として、L2ノルムを採用したケース。
model2 と比べると、過学習が抑えられていることがわかる。
ハイパーパラメータは0.0002 程度で十分であることも分かった。
ただ、 model1 よりは決定係数が小さくなるのかと思わなくもない。
学習データ以外の未知データでモデルの正確さは判断するべきだということは実感できた。

ロジスティック回帰

分類問題を解く教師あり機械学習モデル。
目的変数が2値(クラスに入る or 入らない) なのが特徴。

目的変数 Y=1となる確率として、実測値ペアから以下のように定義する。(Yはベルヌーイ分布で確率が以下のようになると仮定する)

p_i = \sigma ( w_0 + w_1 * x_{i1} + \cdots + w_m * x_{im})

シグマはシグモイド関数である。

\sigma(x) = \frac{1}{1+exp(-ax)}

これらの p から、尤度関数Lを作り、最尤推定法で w を推定する。
解析的には解けないため、勾配降下法を用いる。
-L としているのは極大と極小を入れ替えてみただけで大きな意味はない。

w(k+1) = w(k) - \alpha \frac{\partial (-L(w))}{\partial w}

この式を繰り返し計算し、

w(k+1) = w(k)

となる、もしくはそれに近しいものになったら終了する。
そのときの w が極小値を与える。つまり、尤度関数の極大値を与える。

実装演習

今回はデータ欠損をテーマにしようと思う。使うデータはタイタニックのやつ。

from google.colab import drive
drive.mount('/content/drive/')

# from モジュール名 import クラス名（もしくは関数名や変数名）
import pandas as pd
from pandas import DataFrame
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression

%matplotlib inline

# titanic data csvファイルの読み込み
titanic = pd.read_csv('/content/drive/My Drive/Colab Notebooks/titanic_train.csv')

# age 欠損値平均値うめ
titanic["AgeFill"] = titanic["Age"].fillna(titanic["Age"].mean())
# gender 2値変換
titanic["Gender"] = titanic["Sex"].map({"female" : 0, "male" : 1}).astype(int)

x_fill = titanic.loc[:, ["AgeFill", "Gender"]]
y = titanic.loc[:, "Survived"]

logistic = LogisticRegression()
logistic.fit(x_fill, y)
print(logistic.predict_proba([[30, 1]]))

# 欠損値削除
t = titanic.loc[:, ["Age", "Gender", "Survived"]].dropna()
x_na = t.loc[:, ["Age", "Gender"]]
y_na = t.loc[:, "Survived"]

logistic_na = LogisticRegression()
logistic_na.fit(x_na, y_na)
print(logistic_na.predict_proba([[30, 1]]))

考察

欠損値を平均で埋めたケースの 30代男性生存率予測値 => 19.331898 %
欠損データを無視したケースの 30代男性生存率予測値 => 21.091533 %

個人的な感覚として、年齢はかなり影響を及ぼしそうであるのにもかかわらず、
欠損値を平均で一律に突っ込むのはまずそうと思っていた。そこで欠損値を平均埋めしたケースと
欠損値があるデータを無視したケースを用意した。
やってみてから気が付いたが、ロジスティック回帰の評価がかなり難しく、どちらが良いとも言えなさそうである。

ここに労力を割くよりは、例えば年齢から大人、子供、老人といった新しい特徴量を作って分析したほうが有意義な気がする。

主成分分析

データの次元を削減するための手法。

例えば、国語、数学、理科、社会、英語、体育の6つのテストのデータがあるとする。
さらに、それとは別に直近の模試の合計点がわかっているとする。
普通に解析すると6個の説明変数で模試の合計点を目的変数として、回帰すればよさそうだが、
次元を削減してみる。
例えば、
文系科目 = (国語　+ 英語　+ 社会) / 3
理系科目 = (数学 + 理科) /2
とした二つの説明変数だけで目的変数をうまく説明できるかもしれない。

こういった、次元の削減を理論的に行うのが主成分分析である。（上記例は特徴量の抽出といった感じになってしまってはいるが・・・
実際には多すぎる説明変数を削減するために用いられる。

理論は簡単に言うと、データ値の平均値からのずれの行列に線形変換させたもののうち、
分散が最大になるものをとってくる。

a_j^tVar(X)a_j

を最大化する。

個人的には、理論にこちらの工夫の余地があまりなさそうで、どう使うかのほうが重要だと
思っているため、理論はこの辺で終わる。
せいぜい分散が最大になるという部分をいじれるかどうかだが、あまりいい案が思いつかない。

実装演習

次元圧縮したデータと、しないデータを比較して、どのくらい変化があるのかを見てみるのが今回の目的。
使用データは乳がん検査データ。

from pandas import DataFrame
import numpy as np
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression

%matplotlib inline

cancer = load_breast_cancer()

# 説明変数追加
df = DataFrame(data=cancer.data,columns=cancer.feature_names)
# 目的変数追加
df['Y'] = np.array(cancer.target)

target = df.loc[:,"Y"]
data = df.loc[:, "mean radius":]

# 学習用とテスト用でデータを分離 25%test
X_train, X_test, y_train, y_test = train_test_split(data, target, test_size=0.25)

# 標準化 (x-m)/h
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)

logistic = LogisticRegression()
logistic.fit(X_train, y_train)
y_pred_before = logistic.score(X_test, y_test)

pca = PCA(n_components = 2)
pca_train = pca.fit_transform(X_train)
pca_test = pca.fit_transform(X_test)

logistic.fit(pca_train, y_train)
y_pred_after = logistic.score(pca_test, y_test)

print(y_pred_before, y_pred_after)

考察

なにもしない = 0.993006993006993
次元圧縮2次元 = 0.916083916083916

次元圧縮して2次元にしたにも関わらず、これほどの精度がでていておどろいた。
今回から学習データとテストデータを分離してモデルを評価するホールドアウト法(上位互換としてクロスバリデーション法がある)を用いてみたが、これまた予想外に精度が良すぎる。
今回の主題とはそれるため追調査はやめておくが、そんなに回帰に当てはまりやすいデータだったのだろうか・・・？

サポートベクターマシーン