ラビットチャレンジレポート：ステージ2 機械学習
1.線形回帰モデル
1-1 要点
ある入力（説明変数）から出力（目的変数）の振る舞いを直線で予測するモデル。入力を$${\bm{x}}$$出力を$${y}$$とし、$${\bm{x}}$$は以下のように表される。
$$
\bm{x}=(x_1, x_2,\dots, x_m)^\mathsf{T}
$$
重みパラメータを$${\bm{w}}$$、モデルの予測値を$${\^{y}}$$とすると
予測値$${\^{y}}$$は入力$${\bm{x}}$$とパラメータ$${\bm{w}}$$の
線形結合により、以下のように表される。
$$
{\^{y}}={\bm{w}}^\mathsf{T}{\bm{x}}+w_0
$$
ここで、出力$${y}$$は回帰直線に誤差が加わった$${{\bm{w}}^\mathsf{T}{\bm{x}}+w_0+\epsilon}$$で表されると仮定し、予測値$${\^{y}}$$との誤差が最小になるようにパラメータを決定することでモデル式を作る。
最小にする誤差関数としてMSE（平均二乗誤差）を用いる。（最小二乗法）
$$
MSE=\frac{1}{n}\displaystyle\sum_{i=1}^n(\^{y_i}-y_i)^2
$$
1-2 実装演習
ボストン住宅データセットを用いて線形単回帰分析を行った。
#ライブラリの導入
from sklearn.datasets import load_boston
from sklearn.linear_model import LinearRegression
from pandas import DataFrame
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt#ボストン住宅データの読み込みと、データ確認。
boston = load_boston()
print(boston)
print(boston.keys())ボストン住宅データ出力結果（一部省略）
{'data': array([[6.3200e-03, 1.8000e+01, 2.3100e+00, ..., 1.5300e+01, 3.9690e+02,
4.9800e+00],
[2.7310e-02, 0.0000e+00, 7.0700e+00, ..., 1.7800e+01, 3.9690e+02,
9.1400e+00],
[2.7290e-02, 0.0000e+00, 7.0700e+00, ..., 1.7800e+01, 3.9283e+02,
4.0300e+00],
...,
[6.0760e-02, 0.0000e+00, 1.1930e+01, ..., 2.1000e+01, 3.9690e+02,
5.6400e+00],
[1.0959e-01, 0.0000e+00, 1.1930e+01, ..., 2.1000e+01, 3.9345e+02,
6.4800e+00],
[4.7410e-02, 0.0000e+00, 1.1930e+01, ..., 2.1000e+01, 3.9690e+02,
7.8800e+00]]), 'target': array([24. , 21.6, 34.7, 33.4, 36.2, 28.7, 22.9, 27.1, 16.5, 18.9, 15. ,
18.9, 21.7, 20.4, 18.2, 19.9, 23.1, 17.5, 20.2, 18.2, 13.6, 19.6,
15.2, 14.5, 15.6, 13.9, 16.6, 14.8, 18.4, 21. , 12.7, 14.5, 13.2,
13.1, 13.5, 18.9, 20. , 21. , 24.7, 30.8, 34.9, 26.6, 25.3, 24.7,
21.2, 19.3, 20. , 16.6, 14.4, 19.4, 19.7, 20.5, 25. , 23.4, 18.9,ボストン住宅データの項目出力結果
dict_keys(['data', 'target', 'feature_names', 'DESCR', 'filename', 'data_module'])ボストン住宅データをデータフレームに変換した後、ターゲット（住宅価格）列を追加。最初の5行を確認した。
1行目はデータ出力が省略されて表示されないようにしている。
pd.set_option('display.max_columns', 14)
df = DataFrame(data=boston.data, columns=boston.feature_names)
df['PRICE'] = np.array(boston.target)
print(df.head())データフレーム出力結果
CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX \
0 0.00632 18.0 2.31 0.0 0.538 6.575 65.2 4.0900 1.0 296.0
1 0.02731 0.0 7.07 0.0 0.469 6.421 78.9 4.9671 2.0 242.0
2 0.02729 0.0 7.07 0.0 0.469 7.185 61.1 4.9671 2.0 242.0
3 0.03237 0.0 2.18 0.0 0.458 6.998 45.8 6.0622 3.0 222.0
4 0.06905 0.0 2.18 0.0 0.458 7.147 54.2 6.0622 3.0 222.0 'RM'（部屋数）を説明変数、'PRICE'（住宅価格）を目的変数に取り、単回帰分析を行った。
#線形単回帰分析
data = df.loc[:,['RM']].values
target = df.loc[:,'PRICE'].values
model = LinearRegression()
model.fit(data,target)出来上がったモデルを使用し、部屋数1～9までの予測値を出力した。
for room in range(1,10):
print(f'予測（部屋数{room}）：{model.predict([[room]])}')出力結果
予測（部屋数1）：[-25.5685118]
予測（部屋数2）：[-16.46640281]
予測（部屋数3）：[-7.36429383]
予測（部屋数4）：[1.73781515]
予測（部屋数5）：[10.83992413]
予測（部屋数6）：[19.94203311]
予測（部屋数7）：[29.04414209]
予測（部屋数8）：[38.14625107]
予測（部屋数9）：[47.24836005]部屋数3以下の時に住宅価格がマイナスになってしまった。
このようになってしまう原因として
①線形モデルを仮定することに無理がある
②用いたデータに何かしらの過不足があった
といったことが考えられる。
①については、部屋数が少ない範囲では、部屋数が少なくなるほど、0でないいくらかの価格に向かって収束していくのが自然であると思われる。また、部屋数が大きくなった場合も、直線的に価格が上がり続けるのは不自然であるように思われる。
②については、部屋数だけでは価格を説明するのに不十分であることや、データフレーム変換後のデータを見る限り、部屋数のデータそのものが6～7の辺りに偏っている可能性がある。
②の検証のため使用したdataとtargetの散布図を確認した。
plt.scatter(data, target, facecolor="none", edgecolor="b", s=50)
plt.show()
これを見ると、やはりデータは部屋数6～7の辺りに偏っているようである。また、住宅価格が50＄以上の場合、全て50.0＄としてデータ入力されていると思われ、部屋数4、部屋数9の辺りでは外れ値として処理した方が適切であると思われるデータも見られる。
ただ、こういったイレギュラーなデータを除外し、適用範囲を限定するのであれば線形モデルでの説明は妥当かもしれない。
このことから、線形単回帰の精度を改善するのであれば、住宅価格50＄以上のデータと外れ値は除外した上で、部屋数も5～8の間までで限定して回帰式を用いるのが良いのではないかと思われる。
2.非線形回帰モデル
2-1 要点
線形では捉えることができない、複雑な構造を持った現象に対して回帰モデルを当てはめるための手法。入力$${\bm{x}}$$をそのまま線形モデルで扱うのではなく、一旦、基底関数と呼ばれる非線形の関数$${\phi}$$で変換し、$${\phi}$$と重みパラメータ$${\bm{w}}$$の線形結合によりモデルを作る。
$$
y_i=w_0+\displaystyle\sum_{j=1}^mw_j\phi_j(\bm{x}_i)+\epsilon_i
$$
ここで$${m}$$は基底関数の数である。
基底関数には、多項式、ガウス型基底等がよく用いられる。
パラメータは線形回帰と同じく、$${\^{y}}$$と$${y}$$の二乗誤差の和が最小になるように決定される。
2-2 実装演習
sin関数を用いてノイズの入ったサンプルデータを作成し、これに対して非線形回帰を行った。ここでは、リッジ回帰（α=0.001, α=0.1）、ラッソ回帰（α=0.001）について、過学習抑制の度合いを比較した。
import numpy as np
import matplotlib.pyplot as plt
def true_func(x):
z = np.sin(x)
return z
n = 100
data = 2 * 3.14 * np.random.rand(n).astype(np.float32) #0～2πの範囲でx座標をランダムに100点取得
data = np.sort(data) #昇順にx座標をソート
target = true_func(data) #正解値としてsin(x)の値を取得
noise = 0.3 * np.random.randn(n) #標準正規分布に従うノイズを生成
target += noise #データ点としてtagetにノイズを乗せる
#リッジ回帰（alpha=0.001)
from sklearn.kernel_ridge import KernelRidge
clf = KernelRidge(alpha=0.001, kernel='rbf')
data = data.reshape(-1,1)
clf.fit(data, target)
p_kridge = clf.predict(data)
plt.scatter(data, target, color='blue', label='data')
plt.plot(data, p_kridge, color='orange', linestyle='-', linewidth=3, markersize=6, label='kernel ridge(alpha=0.001)')
#alpha=0.1
clf2 = KernelRidge(alpha=0.1, kernel='rbf')
clf2.fit(data, target)
p_kridge2 = clf2.predict(data)
plt.plot(data, p_kridge2, color='red', linestyle='-', linewidth=3, markersize=6, label='kernel ridge(alpha=0.1)')
#ラッソ回帰
from sklearn.metrics.pairwise import rbf_kernel
from sklearn.linear_model import Lasso
data2 = rbf_kernel(data, data)
lasso = Lasso(alpha = 0.001).fit(data2, target)
p_lasso = lasso.predict(data2)
plt.plot(data, p_lasso, color='green', linestyle='-', linewidth=3, markersize=6, label='lasso(alpha=0.001)')
plt.legend()
plt.show()出力結果
リッジ回帰（α=0.001）では、やや過学習気味の結果が得られた。正則化率を0.1に変更したことで、より滑らかな曲線が得られ、過学習の抑制がみられた。ラッソ回帰では、α=0.001の条件でも、リッジ回帰（α=0.1）と同等に過学習を抑えたモデルが得られた。
次に、データ点の数を100から5000に増やし、α=0.001の条件でリッジ回帰を行った。
#データ点5000
n2 = 5000
data2 = 2 * 3.14 * np.random.rand(n2).astype(np.float32)
data2 = np.sort(data2)
target2 = true_func(data2)
noise = 0.3 * np.random.randn(n2)
target2 += noise
clf3 = KernelRidge(alpha=0.001, kernel='rbf')
data2 = data2.reshape(-1,1)
clf3.fit(data2, target2)
p_kridge3 = clf3.predict(data2)
plt.scatter(data2, target2, color='blue', label='data')
plt.plot(data2, p_kridge3, color='orange', linestyle='-', linewidth=3, markersize=6, label='kernel ridge(alpha=0.001)')
plt.legend()
plt.show()出力結果
データ点を増やすことで過学習が抑えられる様子が確認できた。
3.ロジスティック回帰モデル
3-1 要点
教師あり学習の分類問題を扱う手法である。
入力（説明変数または特徴量と呼ぶ）
$${\bm{x}=(x_1, x_2,\dots, x_m)^\mathsf{T}}$$を基に
出力（目的変数）
$${y\in\lbrace0,1\rbrace}$$
を正しく分類するモデルを作ることが目的となる。
分類の結果は、パラメータ$${\bm{w}}$$を用いた$${\bm{x}}$$の線形結合をシグモイド関数$${\sigma}$$に入力することで、「1と分類される確率」という形で出力される。シグモイド関数は以下の形で表され、その値は0～1を取る。
$$
\sigma(x)=\frac 1 {1+\exp(-ax)}
$$
パラメータの決定には最尤推定を用いる。
最尤推定では、尤度関数を$${L(\bm{w})=\displaystyle\prod_{i=1}^n\sigma(\bm{w}^\mathsf{T}\bm{x}_i)^{y_i}(1-\sigma(\bm{w}^\mathsf{T}\bm{x}_i))^{1-y_i}}$$
とし、$${E(\bm{w})=-\log(L(\bm{w}))}$$が最小となる$${\bm{w}}$$を求める。
3-2 実装演習
アイリスデータセットに対してロジスティック回帰での分類を行った。
ライブラリのインポート
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn import datasets, linear_model
from sklearn.model_selection import train_test_splitアイリスデータセットをダウンロードし、データフレームに変換。
最初の10行と、seabornでのpairplotを確認。
iris = datasets.load_iris()
df = pd.DataFrame(iris.data, columns = iris.feature_names)
df['target'] = np.array(iris.target)
pd.set_option('display.max_columns', 5)
print(df.head(10))
sns.pairplot(df, hue = 'target')
plt.show()出力結果
sepal length (cm) sepal width (cm) petal length (cm) petal width (cm) \
0 5.1 3.5 1.4 0.2
1 4.9 3.0 1.4 0.2
2 4.7 3.2 1.3 0.2
3 4.6 3.1 1.5 0.2
4 5.0 3.6 1.4 0.2
5 5.4 3.9 1.7 0.4
6 4.6 3.4 1.4 0.3
7 5.0 3.4 1.5 0.2
8 4.4 2.9 1.4 0.2
9 4.9 3.1 1.5 0.1
target
0 0
1 0
2 0
3 0
4 0
5 0
6 0
7 0
8 0
9 0 
データセットを訓練データとテストデータに分割。
データフレームのtarget列はデータが昇順に並んでいると思われるため、データをシャッフルしている。
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, shuffle=True)訓練データの中身とデータ数を確認。
print(X_train)
print(y_train)
print(f'訓練データ数：{len(X_train)}')
print(f'テストデータ数：{len(X_test)}')出力結果（一部省略）
[[7.7 2.8 6.7 2. ]
[5. 3.3 1.4 0.2]
[5.9 3. 4.2 1.5]
[5.1 3.3 1.7 0.5]
・・・
[5.7 2.8 4.1 1.3]
[4.6 3.4 1.4 0.3]
[6.4 3.1 5.5 1.8]]
[2 0 1 0 1 0 0 0 1 0 0 1 1 1 2 2 0 0 1 0 1 2 1 2 0 2 1 1 0 0 2 0 0 2 2 2 2
2 2 1 1 0 1 1 1 2 0 2 1 1 0 1 1 0 2 1 2 0 2 0 0 2 2 2 0 2 2 0 1 0 1 2 1 1
1 2 0 0 2 0 1 1 2 2 1 0 1 1 0 0 1 2 2 1 1 2 0 2 2 2 2 0 0 1 0 0 1 2 0 1 0
2]
訓練データ数：112
テストデータ数：38y_trainのデータを見ると、きちんとシャッフルされ、データ数も112個と38個に分割されている。
次に、モデルの学習と、テストデータでの精度検証を行った。
clf = linear_model.LogisticRegression()
clf.fit(X_train,y_train)
acc = clf.score(X_test, y_test)
print(f'精度：{acc:.4f}')出力結果
精度：0.9737精度97.4%であった。
4.主成分分析
4-1 要点
多変量データの持つ情報を、より少ない個数の変数に圧縮して表現することを目的とする（教師無し学習の次元削減）。これにより、分析者にとって分かりやすい形に変換したり、グラフ上にプロットしての可視化が可能になる。
ただし変数を減らすため、幾分かの情報は減損する。損失を抑えるため、変換後の軸はデータの分散が大きいものから順に選んでいく。
4-2 実装演習
乳がんデータセットのデータを、主成分分析を用いて2次元のデータに削減し、全特徴量を用いたロジスティック回帰と2次元データを用いたロジスティック回帰で比較した。
ライブラリの導入
import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.datasets import load_breast_cancer
from sklearn import linear_model
from sklearn.model_selection import train_test_split乳がんデータセットの導入。
データフレームに変換し、データの頭5行、悪性・良性のデータ数を確認。
breast_cancer = load_breast_cancer()
pd.set_option('display.max_columns', 31)
df = pd.DataFrame(data=breast_cancer.data, columns=breast_cancer.feature_names)
df['diagnosis'] = np.array(breast_cancer.target)
print(df.head())
B = df['diagnosis'].sum()
M = len(df['diagnosis'])-B
print(f'良性：{B}')
print(f'悪性：{M}')出力結果
mean radius mean texture mean perimeter mean area mean smoothness
0 17.99 10.38 122.80 1001.0 0.11840
1 20.57 17.77 132.90 1326.0 0.08474
2 19.69 21.25 130.00 1203.0 0.10960
3 11.42 20.38 77.58 386.1 0.14250
4 20.29 14.34 135.10 1297.0 0.10030
mean compactness mean concavity mean concave points mean symmetry
0 0.27760 0.3001 0.14710 0.2419
1 0.07864 0.0869 0.07017 0.1812
2 0.15990 0.1974 0.12790 0.2069
3 0.28390 0.2414 0.10520 0.2597
4 0.13280 0.1980 0.10430 0.1809
mean fractal dimension radius error texture error perimeter error
0 0.07871 1.0950 0.9053 8.589
1 0.05667 0.5435 0.7339 3.398
2 0.05999 0.7456 0.7869 4.585
3 0.09744 0.4956 1.1560 3.445
4 0.05883 0.7572 0.7813 5.438
area error smoothness error compactness error concavity error
0 153.40 0.006399 0.04904 0.05373
1 74.08 0.005225 0.01308 0.01860
2 94.03 0.006150 0.04006 0.03832
3 27.23 0.009110 0.07458 0.05661
4 94.44 0.011490 0.02461 0.05688
concave points error symmetry error fractal dimension error
0 0.01587 0.03003 0.006193
1 0.01340 0.01389 0.003532
2 0.02058 0.02250 0.004571
3 0.01867 0.05963 0.009208
4 0.01885 0.01756 0.005115
worst radius worst texture worst perimeter worst area worst smoothness
0 25.38 17.33 184.60 2019.0 0.1622
1 24.99 23.41 158.80 1956.0 0.1238
2 23.57 25.53 152.50 1709.0 0.1444
3 14.91 26.50 98.87 567.7 0.2098
4 22.54 16.67 152.20 1575.0 0.1374
worst compactness worst concavity worst concave points worst symmetry
0 0.6656 0.7119 0.2654 0.4601
1 0.1866 0.2416 0.1860 0.2750
2 0.4245 0.4504 0.2430 0.3613
3 0.8663 0.6869 0.2575 0.6638
4 0.2050 0.4000 0.1625 0.2364
worst fractal dimension diagnosis
0 0.11890 0
1 0.08902 0
2 0.08758 0
3 0.17300 0
4 0.07678 0
良性：357
悪性：212良性・悪性とも分類に使うためには十分なデータ数。
まず、全特徴量を使ってロジスティック回帰。
#全データでロジスティック回帰
X = df.drop(columns='diagnosis')
y = df['diagnosis']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, shuffle=True)
clf = linear_model.LogisticRegression()
clf.fit(X_train,y_train)
acc = clf.score(X_test, y_test)
print(f'精度：{acc:.4f}')精度：0.9371次に2次元にデータ削減してロジスティック回帰。
pca = PCA(n_components = 2)
X_train2 = pca.fit_transform(X_train)
X_test2 = pca.fit_transform(X_test)
clf2 = linear_model.LogisticRegression()
clf2.fit(X_train2, y_train)
acc2 = clf2.score(X_test2, y_test)
print(f'精度：{acc2:.4f}')精度：0.8601精度は93.7％から86.0％に落ちた。
乳がんの検診ということから考えて好ましくない結果であり
この例においては2次元への削減はするべきではないと思われる。
4-3 関連する内容
学習データを$${\bm{x}=(x_1, x_2,\dots, x_m)^\mathsf{T}}$$、その平均を$${\bm{\=x}}$$、データ行列を
$${\bm{\=X}=(\bm{x}_1-{\bm{\=x}},\dots,\bm{x}_n-{\bm{\=x}})^\mathsf{T}}$$とすると、分散共分散行列$${\frac 1 n\bm{\=X}^\mathsf{T}\bm{\=X}}$$の固有値が射影軸の分散になっており、大きい順にk番目の固有値に対応する固有ベクトルを第k主成分と呼ぶ。
第k主成分の分散の、全分散に対する割合を寄与率と呼び、第k主成分が持つ情報の割合（その軸で対象となる現象を説明できる割合）を意味する。
分散の大きいものから順に選んでいるので、これは情報をより説明できる成分から優先して選んでいると言える。
また、第1～kまでの主成分の分散を足し合わせたものの全分散に対する割合を累積寄与率と呼び、第k主成分までの情報で、元々の情報の何％を説明できるかを示す数値となる。
5.サポートベクターマシン（SVM）
5-1 要点
サポートベクターマシンは回帰・分類を行う手法の一つであり、分類問題に用いられることが多い。分類のために、データ点間に直線（3次元空間では平面、4次元以上では超平面）を引く。この時、境界に最も近いデータ点をサポートベクトル、境界とサポートベクトルとの距離をマージンと呼び、境界はマージンを最大化するように決定される。
分類誤りを許容しないものをハードマージンと言う。しかし、これは完全に直線（または超平面）で分離できるデータでのみ実現できる。現実にはこのようなデータが得られることは少なく、ある程度の分類誤りを許容する必要が出てくる。これをソフトマージンと呼ぶ。
また、そもそも線形の境界では上手く分類できそうもない場合も考えられる。（例：下図）
このような場合、一旦データを非線形な関数に移すことで線形分離が可能になる場合がある。（例：下図）
5-2 実装演習
乳がんデータセットを、SVMを用いて分類した。
その際、正則化係数（c）の値によって検証データの正解率がどのように変化するかを見た。
ライブラリの導入。
import numpy as np
import pandas as pd
from sklearn.datasets import load_breast_cancer
from sklearn.svm import LinearSVC
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split乳がんデータセットの導入とデータフレーム変換。
breast_cancer = load_breast_cancer()
pd.set_option('display.max_columns', 31)
df = pd.DataFrame(data=breast_cancer.data, columns=breast_cancer.feature_names)
df['diagnosis'] = np.array(breast_cancer.target)データを訓練データとテストデータに分割。
X = df.drop(columns='diagnosis')
y = df['diagnosis']
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, shuffle=True)データの値を標準化。
scaler = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)正則化係数cの値を0.00001から1000まで変化させながら訓練データの精度と検証データの精度の変化を見た。
cが大きくなるほど、ハードマージンに近い条件になっている。
c = 0.00001
while c < 1001:
svc = LinearSVC(C = c)
svc.fit(X_train, y_train)
acc_train = svc.score(X_train, y_train)
acc_test = svc.score(X_test, y_test)
print(f'精度（訓練データ）：{acc_train:.4f}, 精度（検証データ）：{acc_test:.4f}, c={c}')
c *= 10出力結果
精度（訓練データ）：0.9296, 精度（検証データ）：0.9580, c=1e-05
精度（訓練データ）：0.9437, 精度（検証データ）：0.9650, c=0.0001
精度（訓練データ）：0.9671, 精度（検証データ）：0.9860, c=0.001
精度（訓練データ）：0.9836, 精度（検証データ）：0.9860, c=0.01
精度（訓練データ）：0.9883, 精度（検証データ）：0.9790, c=0.1
精度（訓練データ）：0.9883, 精度（検証データ）：0.9720, c=1.0
精度（訓練データ）：0.9906, 精度（検証データ）：0.9720, c=10.0
精度（訓練データ）：0.9836, 精度（検証データ）：0.9441, c=100.0
精度（訓練データ）：0.9883, 精度（検証データ）：0.9720, c=1000.0訓練データの精度はcが大きくなるほど同じように大きくなるのかと思っていたがそうではなかった。検証データの精度については、c=0.001と0.01の時に最大になった。何度か、元のデータを分割し直して同じように精度の変化を見てみたが、今回は、ほとんどの場合c=0.001またはc=0.01で精度最大となった。
この記事が気に入ったらサポートをしてみませんか？