A. 講義内容
1. 線形回帰モデル
1. 概要
回帰問題を解くための機械学習（教師あり学習）のモデルの一つで、入力とm次元パラメータの線形結合を出力する。
パラメータ
・モデルに含まれる推定すべき未知のベクトル
・特徴量が予測値に対してどのように影響を与えるかを決定する重みの集合
・切片w0はxの0番目の要素が1であると、パラメータwの要素ととらえることができる
・最適なパラメータwを求めることが線形回帰モデルによる機械学習の目的といえる
\boldsymbol{w}=(w_{1}, w_{2}, \cdots, w_{m})^{T} \in \mathbb{R}^{m}
線形結合
・入力とパラメータの内積
・入力のベクトルが多次元（重回帰）の場合でも出力は1次元（スカラ）となる
\hat{y}=\boldsymbol{w}^{T}\boldsymbol{x}+w_{0}=\sum^{m}_{j=1} w_{j}x_{j}+w_{0}
幾何学的意味
説明変数が1次元（単回帰）の場合、下記の線形回帰モデルは下記の式となる。幾何学の観点から解釈すると、図のように表現できる。
y=w_{0}+w_{1}x_{1}+\epsilon
2. データの分割とモデルの汎化性能測定
機械学習ではモデルの汎化性能を評価するためにデータを学習用データと、検証用データに分割し、学習用データでモデルを学習させた後、検証用データで精度を測定する。汎化性能とは、未知のデータに対する予測性能のことであり、検証用データはこの性能を評価するために利用するため、モデルの学習時には使用してはならない。
3. パラメータの推定
線形回帰モデルでは、パラメータの推定に最小二乗法を用いる。最小二乗法は、データとモデル出力の二乗誤差の和（MSE）が最小となるようなパラメータを求める推定手法である。具体的な計算方法としては、MSEのパラメータwによる偏微分を取り、これが0となるパラメータwの値を求める。
2. 非線形回帰モデル
1. 概要
データ構造を線形で捉えられる場合は限られており、特に複雑な非線形構造を内在する現象を説明するには、非線形回帰モデリングを適用する必要がある.
基底展開法
回帰関数として、基底関数と呼ばれる既知の非線形関数とパラメータベクトルの線型結合を使用するモデリング法。下記の基底関数がよく利用される。
・多項式関数
・ガウス型基底関数
・スプライン関数/Bスプライン関数
y_{i}=f(\boldsymbol{w}_{i})+\epsilon_{i} \qquad y_{i}=w_{0}+\sum^{m}_{j=1}w_{j} \phi_{j}(\boldsymbol{w}_{i})+\epsilon_{i}
2. パラメータの推定
線形回帰モデルでは説明変数として、xiを用いていたが、非線形回帰モデルでは、例えばxikやsinxiなどの非線形関数に置き換えて、予測値を計算する。一方、パラメータwnについては非線形回帰モデルでも、線形回帰モデルでも、線形である点は共通しているため、同じ枠組みでパラメータの推定を行うことができる。
\begin{align}
説明変数&：\boldsymbol{x}_{i}=(x_{i1},x_{i2},\cdots,x_{im}) \in \mathbb{R}^{m} \\
非線形関数ベクトル&：\boldsymbol{\phi}(\boldsymbol{x}_{i})=
(\phi_{1}(\boldsymbol{x}_{i}),\phi_{2}(\boldsymbol{x}_{i}),\cdots,\phi_{k}(\boldsymbol{x}_{i}))^{T}
\in \mathbb{R}^{k} \\
非線形関数の計画行列&：\boldsymbol{\Phi}^(train)=
(\boldsymbol{\phi}(\boldsymbol{x}_{1}),\boldsymbol{\phi}(\boldsymbol{x}_{2}),\cdots,\boldsymbol{\phi}(\boldsymbol{x}_{n}))^{T}
\in \mathbb{R}^{n \times k} \\
最尤法による予測値&：\hat{y}=\boldsymbol{\Phi}(\boldsymbol{\Phi}^{(train)T}\boldsymbol{\Phi}^(train))^{-1}
\boldsymbol{\Phi}^{(train)T}y^{(train)}
\end{align}
3. 未学習と過学習
未学習
学習データに対して、十分に小さな誤差が得られない状態を未学習という。
モデルの表現力の低さが原因と考えられるため、より表現力の高いモデルに変更することが対策となる。
過学習
学習データに対して小さな誤差は得られたが、検証データに対しての誤差が大きい状態を過学習という。
過学習への対策
過学習への対策は下記がある。
・学習データの数を増やす
・不要な基底関数(変数)を削除して表現力を抑止
・正則化法を利用して表現力を抑止
不要な基底関数の削除
モデルの複雑さは基底関数の数、位置やバンド幅に依存する。解きたい問題の複雑さに対して、より多くの基底関数を用いると過学習となりやすいため、cvなどを利用して適切な基底関数に絞る方法。
正則化
モデルのパラメータ推定では、モデルによる推定値と訓練データの差が小さくなるよう最小化問題を解くが、正則化法は最小化したい関数に、正則化項を加える方法。正則化項はモデルの複雑さに伴って値が大きくなるため、最小化するにあたってペナルティーとして働く。正則化法の例として下記2つがある。正則化パラメータはモデルの曲線の滑らかさに寄与し、この値が大きくなると制約面が小さくなる。
Ridge
・l2ノルムを利用
・パラメータを0に近づけるよう推定
・縮小推定
Lasso
・l1ノルムを利用
・いくつかのパラメータが正確に0となるよう推定
・スパース推定
未学習/過学習の判定方法
・訓練誤差もテスト誤差もどちらも小さい → 汎化しているモデルの可能性
・訓練誤差は小さいがテスト誤差が大きい → 過学習
・訓練誤差もテスト誤差もどちらも小さくならない → 未学習
4. モデルの評価
ホールドアウト法
有限のデータを学習用とテスト用の2つに分割し、「予測精度」や「誤り率」を推定する方法。手元にデータが大量にある場合を除いて、良い性能評価を与えないという欠点がある。
・学習用データを多くすればテスト用データが減り学習精度は良くなるが、性能評価の精度は悪くなる
・逆にテスト用データを多くすれば学習用データが減少するので、学習そのものの精度が悪くなることになる
クロスバリデーション(交差検証)
有限のデータを、k個に分割し、「予測精度」や「誤り率」を推定する方法（k-分割交差検証）。k個に分割されたデータのうち、1群を検証用データとし、残りを学習用データとして、「予測精度」や「誤り率」を求める。この操作を、検証用のデータ群を変えながらk回繰り返し、それぞれのケースから得られた結果を平均し、全体の予測精度とする。
3. ロジスティック回帰モデル
1. 概要
分類問題を解くための教師あり機械学習モデルの1つで、入力とm次元パラメータの線形結合（上記の線形回帰モデルの式）をシグモイド関数の入力とし、y=1のクラスに分類される確率を出力する。
シグモイド関数
下記の関数をシグモイド関数という。シグモイド関数には次の性質がある。
\sigma(x)=\frac{1}{1+exp(-ax)}
・実数を入力としてとり、0～1の値を出力する
・単調増加関数
・シグモイド関数の微分は、シグモイド関数自身で表現することが可能
\frac{\partial{\sigma(x)}}{\partial{x}}=a\sigma(x)(1-\sigma(x))
2. パラメータの推定
ロジスティック回帰モデルのパラメータ推定では、最尤法を用いる。最尤法は得られたデータから、そのデータを生成したであろう尤もらしいパラメータを推定する方法である。得られたデータを固定し、求めたいロジスティック回帰モデルのパラメータを変化（変数とする）させて、パラメータの尤もらしさを表現する関数を尤度関数と呼ぶ。最尤法では、この関数の最大値を与える変数を、回帰モデルの最適パラメータとする。ロジスティック回帰モデルの尤度関数は下記となる。
y1～ynのデータが得られたとする \\
P(y_{1},y_{2},\cdots,y_{n}|w_{0},w_{1},\cdots,w_{m})=L(\boldsymbol{w})=
\prod_{i=1}^{n}
\sigma(\boldsymbol{w}^{T}\boldsymbol{x}_{i})^{y_{i}}(1-\sigma(\boldsymbol{w}^{T}\boldsymbol{x}_{i}))^{1-y_{i}}
対数尤度関数
尤度関数を最大にする変数を求める問題は、尤度関数の対数を取り、かつ符号を反転したもの（対数尤度関数）を最小化する変数を求める問題に読み替えて計算する。尤度関数は同時確率からなり、積の計算が多いため、対数をとることでそれらを和に変換して表現し、計算を簡単にする。対数関数は単調増加関数であるため、最大値を与える変数の値は変化しない。また、符号の反転は上記の線形回帰モデルにおける最小二乗法ではMSEを最小化したが、そこに合わせて最小化問題にするためにとられる措置である。
\begin{align}
E(w_{0},w_{1},\cdots,w_{m})&=-logL(w_{0},w_{1},\cdots,w_{m}) \\
&=-\sum^{n}_{i=1}(y_{i}logp_{i}+(1-y_{i})log(1-p_{i}))
\end{align}
勾配降下法
線形回帰モデルではMSEの最小化を行う際に、MSEを求めたいパラメータで偏微分して、微分係数＝0の方程式を解いていた。ロジスティック回帰モデルの最尤法でも、同様に対数尤度関数の最小化を計算するが、MSEの場合のように解析的に解くことは困難である。そこで、反復学習によりパラメータを逐次的に更新するアプローチの一つである勾配降下法を用いて計算を行う。
\boldsymbol{w}^{(k+1)}=\boldsymbol{w}^{(k)}-\eta \frac{\partial E(\boldsymbol{w})}{\partial \boldsymbol{w}}
確率的勾配降下法（SGD）
しかし、通常の勾配降下法（最急降下法）は、パラメータを更新するのにN個全てのデータに対する和を求める必要があるため、大きなメモリリソースを必要としたり、計算に膨大な時間を要するなどの欠点がある。そこで、データを一つずつランダムに(「確率的」に)選んでパラメータを更新する、確率的勾配降下法が用いられることがある。最急降下法でパラメータを1回更新するのと同じ計算量でパラメータをn回更新できるので、効率よく最適な解を探索できる。
\boldsymbol{w}^{(k+1)}=\boldsymbol{w}^{(k)}-\eta (y_{i}-p_{i})\boldsymbol{x_{i}}
3. 分類タスクの評価
分類タスクの予測精度は、次の4つの指標を用いて計算される。
・True Positive ：モデルの予測結果がPositiveで、検証用データもPositiveの場合（TP）
・False Positive：モデルの予測結果がPositiveで、検証用データがNegativeの場合（FP）
・True Negative ：モデルの予測結果がNegativeで、検証用データがPositiveの場合（TN）
・False Negative：モデルの予測結果がNegativeで、検証用データもNegativeの場合（FN）
正解率（Accuracy）
・正解した数/予測対象となった全データ数
・単純な正解率はあまり意味をなさないことが多い
適合率（Precision）
・モデルが「Positiveと予測」したものの中で本当にPositiveである割合
・本当にPositiveなものをNegativeとしてしまうことについては考えていない
・「見逃し(False Negative)が多くてもより確実な」予測をしたい際に利用
\frac{TP}{TP+FP}
再現率（Recall）
・「本当にPositiveなもの」の中からPositiveと予測できる割合
・本当はNegativeなものをPositiveとしてしまう事象については考えていない
・「誤り(False Positive)が多少多くても抜け漏れは少ない」予測をしたい際に利
\frac{TP}{TP+FN}
F値
・PrecisionとRecallの調和平均
\frac{2 \times Precision \times Recall}{Precision + Recall}
4. 主成分分析
1. 概要
多変量データの持つ構造を、情報の損失をなるべく小さくしつつ、より少数個の指標に圧縮する手法。少数変数を利用した分析や、2・3次元の場合にデータを可視化するために用いられる。分散の大きさを情報の量と捉えて、線形変換後の変数の分散が最大となる射影軸を探索することで求める。
\begin{align}
学習データ&：\boldsymbol{x}_{i}=(x_{i1},x_{i2},\cdots,x_{im}) \in \mathbb{R}^{m} \\
平均（ベクトル）&：\bar{\boldsymbol{x}}=\frac{1}{n}\sum^{n}_{i=1}\boldsymbol{x}_{i} \\
データ行列&：\bar{\boldsymbol{X}}=(\boldsymbol{x}_{1}-\bar{\boldsymbol{x}},
\boldsymbol{x}_{2}-\bar{\boldsymbol{x}},\cdots,\boldsymbol{x}_{2}-\bar{\boldsymbol{x}})^{T}
\in \mathbb{R}^{n \times m} \\
分散共分散行列&：\sum=Var(\bar{X})=\frac{1}{n}\bar{X}^{T} \bar{X} \\
線形変換後のベクトル&：\boldsymbol{s}_{j}=(s_{1j},s_{2j},\cdots,s_{nj})^{T}=\bar{X}\boldsymbol{a}_{j}
\qquad \boldsymbol{a}_{j} \in \mathbb{R}^{m}
\end{align}
射影軸の探索は、制限付き最適化問題であり、目的関数を最大にする係数ベクトルを探索することで解くことができる。
目的関数：arg \quad \underset{\boldsymbol{a} \in \mathbb{R}^{m}}{\textrm{max}} \quad
\boldsymbol{a}^{T}_{j} Var(\bar{X}) \boldsymbol{a}_{j} \\
制約条件：\boldsymbol{a}^{T}_{j} \boldsymbol{a}_{j}=1
2. 主成分と寄与率
元データがk次元であったとき、第1～第k主成分までの分散の和は、元データの分散と一致する。第k主成分の分散の全分散に対する割合を寄与率と呼ぶ。第1～第k主成分まで圧縮した際の情報損失量の割合を累積寄与率と呼ぶ。
寄与率
c_{k}=\frac{\lambda_{k}}{\sum^{m}_{i=1} \lambda_{i}}
累積寄与率
r_{k}=\frac{\sum^{k}_{j=1} \lambda_{j}}{\sum^{m}_{i=1} \lambda_{i}}
5. アルゴリズム
1. k近傍法
分類問題のための機械学習手法。クラス分類したいデータの最近傍のデータをk個取得し、それらがもっとも多く所属するクラスに識別する。
性質
・kの値を変化させると、それに伴い識別の結果も変化する
・kを大きくすると決定境界は滑らかになる
2. k-means
教師なし学習によるクラスタリング手法で、与えられたデータをk個のクラスタに分類する。
アルゴリズム
k-meansでは下記のアルゴリズムでクラスタリングを行う。
- 各クラスタ中心の初期値を設定する
- 各データ点に対して、各クラスタ中心との距離を計算し、最も距離が近いクラスタを割り当てる
- 各クラスタの平均ベクトル（中心）を計算する
- 収束するまで2, 3の処理を繰り返す
性質
・アルゴリズムの手順1のクラスタ中心の初期値を変えると、クラスタリングの結果も変わることがある
・kの値を変更すると、クラスタリングの結果も変わる
6. サポートベクターマシン
1. 概要
線形モデルの正負の2値による、2クラス分類のための機械学習手法。線形判別関数に最も近いデータ点との距離（マージン）が最大となるように、線形判別関数を求める。
2. 目的関数
線形判別関数は
\boldsymbol{w}^{T}\boldsymbol{x}+b=0
の式で表現される。マージンが最大となるよう、wとbを求めたい。
マージンをκと置くと、すべての点は線形判別式からκ以上離れたところに存在するので、
|\boldsymbol{w}^{T}\boldsymbol{x}_{i}+b|\geqq \kappa
となる。マージンκで正規化し、wとbを改めて置きなおせば、クラスti=+1のとき、
\boldsymbol{w}^{T}\boldsymbol{x}_{i}+b\geqq +1
クラスti=-1のとき、
\boldsymbol{w}^{T}\boldsymbol{x}_{i}+b\leqq -1
となるので、下記の式にまとめることができる。
t_{i}(\boldsymbol{w}^{T}\boldsymbol{x}_{i}+b)\geqq 1
ここで、2クラス間のマージンρ(w, b)は、クラスti=+1のデータのうち、線形判別式との距離が最小のものと、クラスti=-1のデータのうち、線形判別式との距離が最小のものを用いて、
\begin{align}
\rho(\boldsymbol{w},b)&=\underset{x_{i} \in C_{ti}=+1}{\textrm{min}} \frac{\boldsymbol{w}^{T}\boldsymbol{x}_{i}}{||\boldsymbol{w}||} - \underset{x_{i} \in C_{ti}=-1}{\textrm{max}} \frac{\boldsymbol{w}^{T}\boldsymbol{x}_{i}}{||\boldsymbol{w}||} \\
&=\frac{1-b}{||\boldsymbol{w}||} - \frac{-1-b}{||\boldsymbol{w}||} \\
&=\frac{2}{||\boldsymbol{w}||} \\
\end{align}
とあらわすことができる。マージンρ(w, b)が最大となるような線形判別式を求めたいので、まとめると下記の不等式制約条件最適化問題の主問題を解くことで、最適な線形判別式を求めることができることになる。
評価関数（最小化）：L_{p}(\boldsymbol{w})=\frac{1}{2} \boldsymbol{w}^{T} \boldsymbol{w} \\
不等式制約条件：t_{i}(\boldsymbol{w}^{T} \boldsymbol{x}_{i} + b) \geqq 1
上記をラグランジュ未定乗数法によって解くが、不等式制約条件であるため、その解はKKT条件を満たす。
求める最適なwをw0、bをb0とおく。
ラグランジュ関数:
L(\boldsymbol{w},b,\boldsymbol{a})=\frac{1}{2} \boldsymbol{w}^{T} \boldsymbol{w}-\sum_{i=1}^{n}a_{i}(t_{i}(\boldsymbol{w}^{T} \boldsymbol{x}_{i} + b) - 1)
\begin{align}
KKT条件:
&\frac{\partial L(\boldsymbol{w},b,\boldsymbol{a})}
{\partial \boldsymbol{w}}
|_{\boldsymbol{w}=\boldsymbol{w}_{0}}
= 0 \tag{1}\\
&\frac{\partial L(\boldsymbol{w},b,\boldsymbol{a})}
{\partial b}
= 0 \tag{2}\\
&t_{i}(\boldsymbol{w}^{T} \boldsymbol{x}_{i} + b) -1 \geqq 0 \tag{3}\\
&a_{i} \geqq 0 \tag{4}\\
&a_{i}(t_{i}(\boldsymbol{w}^{T} \boldsymbol{x}_{i} + b) - 1) = 0 \tag{5}
\end{align}
KKT条件(1)、(2)から、
\boldsymbol{w}_{0}=\sum_{i=1}^{n}a_{i}t_{i}\boldsymbol{x}_{i}, \quad \sum_{i=1}^{n}a_{i}t_{i}=0
これらを用いて、ラグランジュ関数からwとbを消去し、ラグランジュ未定乗数のみの関数とする。
L(\boldsymbol{a})=\sum_{i=1}^{n}a_{i}-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}a_{i}a_{j}t_{i}t_{j}\boldsymbol{x}^{T}_{i} \boldsymbol{x}_{j}
すると、解くべき問題は、Lを最大にするaを求める問題に置き換えることができ、これを主問題に対する双対問題という。
評価関数（最大化）：L_{d}(\boldsymbol{a})=\sum_{i=1}^{n}a_{i}-\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}a_{i}a_{j}t_{i}t_{j}\boldsymbol{x}^{T}_{i} \boldsymbol{x}_{j} \\
制約条件：\sum_{i=1}^{n}a_{i}t_{i}=0
よって、最大マージンは最適なラグランジュ未定乗数a*を用いて下記となる。
\frac{1}{||\boldsymbol{w}||}=\frac{1}{\sqrt{\sum_{i=1}^{n}a^{*}_{i}}}
3. 非線形分離
ソフトマージンSVM
サンプルを線形分離できないときに、誤差や、マージン内にデータが入ることをを許容し、対してペナルティを与える方法。パラメータCの大小で決定境界が変化する。
カーネルトリック
特徴空間に写像し、その空間で線形に分離する方法。カーネル関数を用いて、高次元ベクトルの内積をスカラー関数で表現する。特徴空間が高次元でも計算コストを抑えられる。
B. 実装演習結果
1. 線形回帰モデル
必要モジュールとデータのインポート
# from モジュール名 import クラス名（もしくは関数名や変数名）
from sklearn.datasets import load_boston
from pandas import DataFrame
import numpy as np
# ボストンデータを"boston"というインスタンスにインポート
boston = load_boston()
# インポートしたデータを確認(data / target / feature_names / DESCR)
print(boston)
# 見慣れない形式だったので型を確認
print(type(boston))
scikit-learn固有の型でした。
# DESCR変数の中身を確認
print(boston['DESCR'])
データの内容が記載されています。
データフレームの作成
# 説明変数らをDataFrameへ変換
df = DataFrame(data=boston.data, columns = boston.feature_names)
# 目的変数をDataFrameへ追加
df['PRICE'] = np.array(boston.target)
# 最初の5行を表示
df.head(5)
線形単回帰分析
# 説明変数
data = df.loc[:, ['RM']].values
# 目的変数
target = df.loc[:, 'PRICE'].values
## sklearnモジュールからLinearRegressionをインポート
from sklearn.linear_model import LinearRegression
# オブジェクト生成
model = LinearRegression()
# model.get_params()
# model = LinearRegression(fit_intercept = True, normalize = False, copy_X = True, n_jobs = 1)
# fit関数でパラメータ推定
model.fit(data, target)
# 予測
model.predict([[1]])
部屋数のみから物件の価格を推定すると、値が負になってしまい、説明変数が適切でないことが考えられます。
重回帰分析
# 説明変数
data2 = df.loc[:, ['CRIM', 'RM']].values
# 目的変数
target2 = df.loc[:, 'PRICE'].values
# オブジェクト生成
model2 = LinearRegression()
# fit関数でパラメータ推定
model2.fit(data2, target2)
model2.predict([[0.2, 7]])
上記より、部屋数だけでは良い予測が得られなかったので、物件周辺の犯罪率も加味してみると、物件価格の値は正となり予測精度は向上したといえます。
回帰係数と切片の値を確認
# 単回帰の回帰係数と切片を出力
print('推定された回帰係数: %.3f, 推定された切片 : %.3f' % (model.coef_, model.intercept_))
線形回帰モデルで最適化した、具体的な一次関数のパラメータは上記であったことがわかります。
モデルの検証
# train_test_splitをインポート
from sklearn.model_selection import train_test_split
# 70%を学習用、30%を検証用データにするよう分割
X_train, X_test, y_train, y_test = train_test_split(data, target,
test_size = 0.3, random_state = 666)
# 学習用データでパラメータ推定
model.fit(X_train, y_train)
# 作成したモデルから予測（学習用、検証用モデル使用）
y_train_pred = model.predict(X_train)
y_test_pred = model.predict(X_test)
# matplotlibをインポート
import matplotlib.pyplot as plt
# Jupyterを利用していたら、以下のおまじないを書くとnotebook上に図が表示
%matplotlib inline
# 学習用、検証用それぞれで残差をプロット
plt.scatter(y_train_pred, y_train_pred - y_train, c = 'blue', marker = 'o', label = 'Train Data')
plt.scatter(y_test_pred, y_test_pred - y_test, c = 'lightgreen', marker = 's', label = 'Test Data')
plt.xlabel('Predicted Values')
plt.ylabel('Residuals')
# 凡例を左上に表示
plt.legend(loc = 'upper left')
# y = 0に直線を引く
plt.hlines(y = 0, xmin = -10, xmax = 50, lw = 2, color = 'red')
plt.xlim([10, 50])
plt.show()
# 平均二乗誤差を評価するためのメソッドを呼び出し
from sklearn.metrics import mean_squared_error
# 学習用、検証用データに関して平均二乗誤差を出力
print('MSE Train : %.3f, Test : %.3f' % (mean_squared_error(y_train, y_train_pred), mean_squared_error(y_test, y_test_pred)))
# 学習用、検証用データに関してR^2を出力
print('R^2 Train : %.3f, Test : %.3f' % (model.score(X_train, y_train), model.score(X_test, y_test)))
学習データ/検証データに対する、平均二乗誤差/決定係数はそれぞれ上記でした。
2. 非線形回帰モデル
モジュールのインポートと設定
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
%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.linear_model import LinearRegression
clf = LinearRegression()
data = data.reshape(-1,1)
target = target.reshape(-1,1)
clf.fit(data, target)
p_lin = clf.predict(data)
plt.scatter(data, target, label='data')
plt.plot(data, p_lin, color='darkorange', marker='', linestyle='-', linewidth=1, markersize=6, label='linear regression')
plt.legend()
print(clf.score(data, target))
線形下記では表現力が不足していることが見受けられます。
カーネルモデル
from sklearn.kernel_ridge import KernelRidge
clf = KernelRidge(alpha=0.0002, kernel='rbf')
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')
plt.legend()
# plt.plot(data, p, color='orange', marker='o', linestyle='-', linewidth=1, markersize=6)
非線形モデルでは、精度よく予測できています。
3. ロジスティック回帰モデル
データ読み込み・クレンジング
# titanic data csvファイルの読み込み
titanic_df = pd.read_csv('/content/drive/My Drive/study_ai_ml_google/data/titanic_train.csv')
# 予測に不要と考えるからうをドロップ (本当はここの情報もしっかり使うべきだと思っています)
titanic_df.drop(['PassengerId', 'Name', 'Ticket', 'Cabin'], axis=1, inplace=True)
# nullを含んでいる行を表示
titanic_df[titanic_df.isnull().any(1)].head(10)
# Ageカラムのnullを中央値で補完
titanic_df['AgeFill'] = titanic_df['Age'].fillna(titanic_df['Age'].mean())
# 冒頭を表示
titanic_df.head()
ロジスティック回帰
# 運賃だけのリストを作成
data1 = titanic_df.loc[:, ["Fare"]].values
# 生死フラグのみのリストを作成
label1 = titanic_df.loc[:,["Survived"]].values
from sklearn.linear_model import LogisticRegression
model=LogisticRegression()
model.fit(data1, label1)
# 予測させる
model.predict([[61]])
# 予測の自信
model.predict_proba([[62]])
0の自信と1の自信は足して1となることが確認できます。
w_0 = model.intercept_[0]
w_1 = model.coef_[0,0]
# def normal_sigmoid(x):
# return 1 / (1+np.exp(-x))
def sigmoid(x):
return 1 / (1+np.exp(-(w_1*x+w_0)))
x_range = np.linspace(-1, 500, 3000)
plt.figure(figsize=(9,5))
# plt.xkcd()
plt.legend(loc=2)
# plt.ylim(-0.1, 1.1)
# plt.xlim(-10, 10)
# plt.plot([-10,10],[0,0], "k", lw=1)
# plt.plot([0,0],[-1,1.5], "k", lw=1)
plt.plot(data1,np.zeros(len(data1)), 'o')
plt.plot(data1, model.predict_proba(data1), 'o')
plt.plot(x_range, sigmoid(x_range), '-')
# plt.plot(x_range, normal_sigmoid(x_range), '-')
#
右に行くほど値が1に近づくシグモイド曲線が確認できます。運賃が高いほど生存率が高いことを示しています。
実装(2変数から生死を判別)
titanic_df['Gender'] = titanic_df['Sex'].map({'female': 0, 'male': 1}).astype(int)
titanic_df['Pclass_Gender'] = titanic_df['Pclass'] + titanic_df['Gender']
np.random.seed = 0
xmin, xmax = -5, 85
ymin, ymax = 0.5, 4.5
index_survived = titanic_df[titanic_df["Survived"]==0].index
index_notsurvived = titanic_df[titanic_df["Survived"]==1].index
from matplotlib.colors import ListedColormap
fig, ax = plt.subplots()
cm = plt.cm.RdBu
cm_bright = ListedColormap(['#FF0000', '#0000FF'])
sc = ax.scatter(titanic_df.loc[index_survived, 'AgeFill'],
titanic_df.loc[index_survived, 'Pclass_Gender']+(np.random.rand(len(index_survived))-0.5)*0.1,
color='r', label='Not Survived', alpha=0.3)
sc = ax.scatter(titanic_df.loc[index_notsurvived, 'AgeFill'],
titanic_df.loc[index_notsurvived, 'Pclass_Gender']+(np.random.rand(len(index_notsurvived))-0.5)*0.1,
color='b', label='Survived', alpha=0.3)
ax.set_xlabel('AgeFill')
ax.set_ylabel('Pclass_Gender')
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
ax.legend(bbox_to_anchor=(1.4, 1.03))
チケットクラスと性別（1, 0）の和をとり、新たな変数Pclass_Genderを作成し、生死を説明しています。縦方向に図を見ていくと、下のほうがプロットが青いことから、女性のほうが、階級が高いほうが生存率が高いことがわかります。右方向に図を見ると右に行くほど、赤いプロットが多くなっており、高齢であるほど生存率が低いことがわかります。
h = 0.02
xmin, xmax = -5, 85
ymin, ymax = 0.5, 4.5
xx, yy = np.meshgrid(np.arange(xmin, xmax, h), np.arange(ymin, ymax, h))
Z = model2.predict_proba(np.c_[xx.ravel(), yy.ravel()])[:, 1]
Z = Z.reshape(xx.shape)
fig, ax = plt.subplots()
levels = np.linspace(0, 1.0)
cm = plt.cm.RdBu
cm_bright = ListedColormap(['#FF0000', '#0000FF'])
# contour = ax.contourf(xx, yy, Z, cmap=cm, levels=levels, alpha=0.5)
sc = ax.scatter(titanic_df.loc[index_survived, 'AgeFill'],
titanic_df.loc[index_survived, 'Pclass_Gender']+(np.random.rand(len(index_survived))-0.5)*0.1,
color='r', label='Not Survived', alpha=0.3)
sc = ax.scatter(titanic_df.loc[index_notsurvived, 'AgeFill'],
titanic_df.loc[index_notsurvived, 'Pclass_Gender']+(np.random.rand(len(index_notsurvived))-0.5)*0.1,
color='b', label='Survived', alpha=0.3)
ax.set_xlabel('AgeFill')
ax.set_ylabel('Pclass_Gender')
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
# fig.colorbar(contour)
x1 = xmin
x2 = xmax
y1 = -1*(model2.intercept_[0]+model2.coef_[0][0]*xmin)/model2.coef_[0][1]
y2 = -1*(model2.intercept_[0]+model2.coef_[0][0]*xmax)/model2.coef_[0][1]
ax.plot([x1, x2] ,[y1, y2], 'k--')
モデルを用いて境界線を引くと図のようになります。
4. 主成分分析
データの読み込み
cancer_df = pd.read_csv('/content/drive/My Drive/study_ai_ml_google/data/cancer.csv')
print('cancer df shape: {}'.format(cancer_df.shape))
# 欠損値は列ごと削除
cancer_df.drop('Unnamed: 32', axis=1, inplace=True)
# 目的変数の抽出
y = cancer_df.diagnosis.apply(lambda d: 1 if d == 'M' else 0)
# 説明変数の抽出
X = cancer_df.loc[:, 'radius_mean':]
# 学習用とテスト用でデータを分離
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)
# 標準化
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)
# ロジスティック回帰で学習
logistic = LogisticRegressionCV(cv=10, random_state=0)
logistic.fit(X_train_scaled, y_train)
# 検証
print('Train score: {:.3f}'.format(logistic.score(X_train_scaled, y_train)))
print('Test score: {:.3f}'.format(logistic.score(X_test_scaled, y_test)))
print('Confustion matrix:\n{}'.format(confusion_matrix(y_true=y_test, y_pred=logistic.predict(X_test_scaled))))
PCAによる次元削減前の32次元のデータで、ロジスティック回帰による精度を測定したものです。
97%を分類できています。
寄与率
pca = PCA(n_components=30)
pca.fit(X_train_scaled)
plt.bar([n for n in range(1, len(pca.explained_variance_ratio_)+1)], pca.explained_variance_ratio_)
各主成分の寄与率を棒グラフにしたものです。
第一主成分と、第二主成分でおよそ60%は説明ができそうです。
可視化
# PCA
# 次元数2まで圧縮
pca = PCA(n_components=2)
X_train_pca = pca.fit_transform(X_train_scaled)
print('X_train_pca shape: {}'.format(X_train_pca.shape))
# X_train_pca shape: (426, 2)
# 寄与率
print('explained variance ratio: {}'.format(pca.explained_variance_ratio_))
# explained variance ratio: [ 0.43315126 0.19586506]
# 散布図にプロット
temp = pd.DataFrame(X_train_pca)
temp['Outcome'] = y_train.values
b = temp[temp['Outcome'] == 0]
m = temp[temp['Outcome'] == 1]
plt.scatter(x=b[0], y=b[1], marker='o') # 良性は○でマーク
plt.scatter(x=m[0], y=m[1], marker='^') # 悪性は△でマーク
plt.xlabel('PC 1') # 第1主成分をx軸
plt.ylabel('PC 2') # 第2主成分をy軸
PCAによりデータを二次元まで圧縮し、可視化したものです。
直感的には第一主成分と、第二主成分でおよそ60%説明できるというのが、正しそうです。
5. アルゴリズム
k近傍法
def gen_data():
x0 = np.random.normal(size=50).reshape(-1, 2) - 1
x1 = np.random.normal(size=50).reshape(-1, 2) + 1.
x_train = np.concatenate([x0, x1])
y_train = np.concatenate([np.zeros(25), np.ones(25)]).astype(np.int)
return x_train, y_train
X_train, ys_train = gen_data()
plt.scatter(X_train[:, 0], X_train[:, 1], c=ys_train)
訓練データを生成しました。
k=3の場合です。
k=10の場合です。3の場合と境界が異なります。
k-means
def gen_data():
x1 = np.random.normal(size=(100, 2)) + np.array([-5, -5])
x2 = np.random.normal(size=(100, 2)) + np.array([5, -5])
x3 = np.random.normal(size=(100, 2)) + np.array([0, 5])
return np.vstack((x1, x2, x3))
# データ作成
X_train = gen_data()
# データ描画
plt.scatter(X_train[:, 0], X_train[:, 1])
データ作成からプロットまでです。
k-meansによるクラスタリング結果です。初期値が近すぎたために、うまく分離できていなかったようです。
再度実行し、うまく分離できました。
6. サポートベクターマシン
実験データを生成しました。
SVMによって完全に分離できました。サポートベクターが４つ存在しています。
線形分離できないパターンの実験データを生成しました。
RBFカーネルを用いたカーネル法により、非線形で分類できました。
完全分離できないパターンの実験データを生成しました。
ソフトマージンSVMによって分類できました。マージン内にもデータが存在しています。
C. 関連記事レポート
1. 線形回帰モデル
2. 非線形回帰モデル
よく使われる既定関数に挙げられていたB-スプライン曲線について調べました。
・スプラインは、製図などに用いられる一種の自在定規のこと
・B-スプライン曲線は与えられた複数の制御点とノットベクトルから定義される滑らかな曲線
・区分多項式（区分的に定義された多項式）により表現されている
・一部を変更しても曲線全体に影響は及ばない性質がある
・ベジェ曲線とともに、コンピュータグラフィックスの世界で広く利用されている
・式は下記となる
\begin{align}
&S(t)=\sum^{m-n-2}_{i=0} \boldsymbol{P}_{i}b_{i,n}(t),\qquad t\in[t_{n},t_{m-n-1}] \\
&b_{j,0}(t)=\left\{
\begin{array}{ll}
1 & (if\quad t_{j} \leq t \lt t_{j+1}) \\
0 & (otherwise)
\end{array}
\right.
,\qquad j=0,\cdots ,m-2 \\
&b_{j,k}(t)=
\frac{t-t_{j}}{t_{j+k}-t_{j}}b_{j,k-1}(t)+\frac{t_{j+k+1}-t}{t_{j+k+1}-t_{j+1}}b_{j+1,k-1}(t)
,\qquad j=0,\cdots ,m-k-2
\end{align}
参考：
https://tajimarobotics.com/basis-spline-interpolation-example-2/
https://ja.wikipedia.org/wiki/B-%E3%82%B9%E3%83%97%E3%83%A9%E3%82%A4%E3%83%B3%E6%9B%B2%E7%B7%9A