機械学習の定番であるscikit-learnを可視化パッケージyellowbrickを用いて解説する。

可視化パッケージYellowbrickのインストール

以下を参照

https://www.scikit-yb.org/en/latest/

Google Colab.上で古いバージョンがインストールされている場合には,以下のコードを実行後に, RUNTIMEのリスタートが必要である.

#!pip install -U yellowbrick
import yellowbrick

yellowbrick.__version__
'1.3.post1'

Warningについて

Yellowbrickとscikit-learnのバージョン違いで注意 (Warning) が出ることがある. また,最適化が収束しないと注意が出ることがある.

鬱陶しい場合には,以下の「おまじない」コードを実行しておく.

from warnings import simplefilter
from sklearn.exceptions import ConvergenceWarning

simplefilter("ignore", category=ConvergenceWarning)
simplefilter(action="ignore", category=FutureWarning)

機械学習とは

機械学習(machine learning)は,Googleの検索エンジン,Facebookの友人紹介,Amazonのおすすめ商品の紹介,スパムメールの分別, ドローンの宙返り,自然言語処理,手書き文字の認識,データマイニングなど様々な応用に使われている. 機械学習は,コンピュータに新しい能力を与えるための基本学問体系であり, コンピュータに知能を与えるための総合学問である人工知能の一分野と位置づけられる. 機械学習の定義には様々なものがあるが,チェッカーゲームに対する初めての学習で有名な Arthur Samuelによれば「コンピュータに明示的なプログラムを書くことなしに学習する能力を与える学問分野」 (``the field of study that gives computers the ability to learn without being explicitly programmed") と定義される.

教師あり学習

機械学習は,大きく分けて教師あり学習(supervised learning)と 教師なし学習(unsupervised learning)に分類される. 教師あり学習は,入力と出力の組から成るデータを与えて, その関係を学習するものであり,教師なし学習は, 入力だけから成るデータから学習を行うものである.

教師あり学習は,大きく回帰(regression)と分類(classification)に分けられる. 回帰は出力が連続な値をもつ場合であり, 分類は離散的な値をもつ場合に対応する. 日にちを入力として各日の商品の需要量を出力としたデータが与えられたとき, 需要量を予測するためには回帰を用いる.需要量は連続な値をもつからだ.一方,顧客の過去の購入履歴から, ある商品を購入するか否かを判定するためには分類を用いる.購入するか否かは $0,1$ の離散的な値で表されるからだ. 回帰や分類を行うための入力を特徴(feature)とよぶ.特徴は1つでも複数でも(場合によっては)無限個でもよい.

教師なし学習

教師なし学習は,出力(正解)が与えられていないデータを用いる機械学習であり, 代表的なものとしてクラスタリング(clustering)や次元削減(dimension reduction)があげられる. クラスタリングは,入力されたデータ(変数)間の関係をもとにデータのグループ化を行う際に用いられ, 次元削減は,高次元データを低次元に落とす際に用いられる. これらの手法は,教師あり学習の前処理として有用である.

線形回帰

線形回帰の仮説関数(hypothesis function)は,

$$ h_w (x)=w_0 +w_1 x_1 + w_2 x_2 + \cdots + w_n x_n $$

と表すことができる.これは,入力 $x$ と出力 $y$ が重み $w$ を用いて $h_w(x)$ の ように表すことができるという仮説を表す.

学習アルゴリズムの目的は,トレーニング集合を用いて, 最も「良い」重みベクトル $w$ を求めることである. それでは,どのような重み $w$ が「良い」のであろうか? 線形回帰では,各データ $i$ に対して,仮説関数によって「予測」された値 $\hat{y}^{(i)} = h_{w}(x^{(i)})$ と本当の値 $y^{(i)}$ の 誤差の自乗(2乗)が小さい方が「良い」と考える.これを損出関数 (loss function)とよぶ. すべてのデータに対する損出関数の平均値(を $2$で割ったもの;これは微分したときに打ち消し合うためである)を費用関数 (cost function) $J(w)$ とよぶ.

$$ J(w)= \frac{1}{2m} \displaystyle\sum_{i=1}^m \left( h_{w}(x^{(i)}) -y^{(i)} \right)^2 $$

費用関数 $J(w)$ を最小にする重みベクトル $w$ を求めることが線形回帰の目的となる.

実は,費用関数 $J(w)$ は凸関数になるので, これを最小にする $w$ は凸関数に対する非線形最適化の手法を用いて簡単に求めることができる.

ベクトルと行列を用いて最適な $w$ を求める公式を導いておこう. トレーニングデータをベクトル $y$ と行列 $X$ で表しておく.

$$ y = \left( \begin{array}{c} y^{(1)} \\ y^{(2)} \\ \vdots \\ y^{(m)} \end{array} \right) $$$$ X = \left( \begin{array}{cccc} 1 & x^{(1)}_{1} & \ldots & x^{(1)}_{n} \\ 1 & x^{(2)}_{1} & \ldots & x^{(2)}_{n} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x^{(m)}_{1} & \ldots & x^{(m)}_{n} \end{array} \right) $$

行列 $X$ に対しては最初の列にすべて $1$ のベクトルを追加していることに注意されたい.

重みベクトル $w$ は $n+1$次元のベクトルとして以下のように表す. $$ w = \left( \begin{array}{c} w_0 \\ w_1 \\ \vdots \\ w_n \end{array} \right) $$

このとき費用関数は,上で定義したベクトルと行列で記述すると $$ J(w)= \frac{1}{2m} (y-Xw)^T (y-Xw) $$ となる.これを $w$ で偏微分したものを $0$ とおくと方程式 $$ X^T (y-X w) =0 $$ が得られるので,$w$ について解くことによって, $$ w =(X^T X)^{-1} X^T y $$ を得る.これが線形回帰の解ベクトルとなる.

線形回帰がどれだけ良いかを表す尺度として決定係数(coefficient of determination) $R^2$ がある. $R^2$ は予測値とトレーニングデータの誤差の自乗和から計算される.

まず,最適な重み $w$ の下での誤差の自乗和(SSE: Sum of Square Error)を以下のように計算する. $$ SSE = \displaystyle\sum_{i=1}^m \left( h_{w}(x^{(i)}) -y^{(i)} \right)^2 $$ 次に,基準となる誤差の自乗和として, 平均値 $\bar{y}$ とトレーニングデータ $y^{(i)}$ の差の 自乗和(SST: Total Sum of Square)を計算する. $$ SST = \displaystyle\sum_{i=1}^m \left( \bar{y} -y^{(i)} \right)^2 $$

決定係数 $R^2$ はこれらの自乗和の比 $SSE/SST$ を $1$ から減じた値と定義される. $$ R^2 = 1 - SSE/SST $$

定義から分かるように $R^2$ は $1$以下の値をもち, $1$に近いほど誤差が少ない予測であることが言える.

広告による売り上げの予測

広告のデータ http://logopt.com/data/Advertising.csv を用いる.

テレビ(TV),ラジオ(Radio),新聞(Newspaper)への広告から売り上げ(Sales)を予測する.

import pandas as pd  # まずはpandasモジュールを準備する.
df = pd.read_csv("http://logopt.com/data/Advertising.csv", index_col=0)  # 0行目をインデックスにする.
df.tail()
TV Radio Newspaper Sales
196 38.2 3.7 13.8 7.6
197 94.2 4.9 8.1 9.7
198 177.0 9.3 6.4 12.8
199 283.6 42.0 66.2 25.5
200 232.1 8.6 8.7 13.4

独立変数(特徴ベクトル)$X$ は TV, Radio, Newspaperの列,従属変数(ターゲット) $y$ は Salesの列

y = df["Sales"]
X = df[["TV", "Radio", "Newspaper"]]  # df.drop("Sales",axis=1) でも同じ
X.head()
TV Radio Newspaper
1 230.1 37.8 69.2
2 44.5 39.3 45.1
3 17.2 45.9 69.3
4 151.5 41.3 58.5
5 180.8 10.8 58.4

scikit-learnの基本手順

  • 手順1:クラスをインポートして,インスタンスを生成する.

  • 手順2:fitメソッドを用いて,データから訓練する.

  • 手順3:predictメソッドを用いて予測を行う.

from sklearn.linear_model import LinearRegression  # 線形回帰クラス LinearRegression のインポート

reg = LinearRegression()  # 線形回帰クラスのインスタンス reg を生成
reg.fit(X, y)  # fitによる訓練
yhat = reg.predict(X)  # predictによる予測
print("y-切片= ", reg.intercept_)
print("係数 = ", reg.coef_)
y-切片=  2.938889369459412
係数 =  [ 0.04576465  0.18853002 -0.00103749]
SSE = ((yhat - y) ** 2).sum()  # Sum of Square Error
SST = ((y.mean() - y) ** 2).sum()  # Total Sum of Square
print("R2 =", 1 - SSE / SST)  # 決定係数 R^2
R2 = 0.8972106381789522
print(reg.score(X, y))  # 決定係数の別計算
0.8972106381789522

可視化(回帰)

  • Yellowbrickパッケージを用いる.

  • 可視化の基本手順

  • 手順1:クラスをインポートして,可視化インスタンスを生成する.

  • 手順2:fitメソッドで,データを用いて訓練する.

  • 手順3:scoreメソッドを用いて評価尺度を計算する.

  • 手順4:showメソッドを用いて図を表示する.

  • 回帰に対しては,以下の2種類がある.

    • 予測誤差 (PredictionError)
    • 残差プロット (ResidualPlot)
from yellowbrick.regressor import PredictionError

visualizer = PredictionError(reg)

visualizer.fit(X, y)
visualizer.score(X, y)
visualizer.show()
<AxesSubplot:title={'center':'Prediction Error for LinearRegression'}, xlabel='$y$', ylabel='$\\hat{y}$'>
from yellowbrick.regressor import ResidualsPlot

visualizer = ResidualsPlot(reg)

visualizer.fit(X, y)
visualizer.score(X, y)
visualizer.show()
<AxesSubplot:title={'center':'Residuals for LinearRegression Model'}, xlabel='Predicted Value', ylabel='Residuals'>

問題 (SAT,GPA)

http://logopt.com/data/SATGPA.csv データを用いて,2種類のSATの成績からGPAを予測せよ. さらに結果の可視化を行え.

問題(住宅価格)

http://logopt.com/data/Boston.csv のBostonの住宅データを用いて回帰分析を行え. さらに結果の可視化を行え.

medvが住宅の価格で,他のデータ(犯罪率や人口など)から予測する.

問題(車の燃費)

http://logopt.com/data/Auto.csv の車の燃費のデータを用いて回帰分析を行え. さらに結果の可視化を行え.

データの詳細については,

https://vincentarelbundock.github.io/Rdatasets/doc/ISLR/Auto.html

を参照せよ.

最初の列が燃費(mpg: Mile Per Gallon)であり,これを他の列の情報を用いて予測する.最後の列は車名なので無視して良い.

car = pd.read_csv("http://logopt.com/data/Auto.csv", index_col=0)
car.head()
mpg cylinders displacement horsepower weight acceleration year origin name
0 18.0 8 307.0 130 3504 12.0 70 1 chevrolet chevelle malibu
1 15.0 8 350.0 165 3693 11.5 70 1 buick skylark 320
2 18.0 8 318.0 150 3436 11.0 70 1 plymouth satellite
3 16.0 8 304.0 150 3433 12.0 70 1 amc rebel sst
4 17.0 8 302.0 140 3449 10.5 70 1 ford torino

問題(コンクリートの強度)

以下のコンクリートの強度の例題に対して,strength列の強度を他の列の情報から,線形回帰を用いて推定せよ. さらに結果の可視化を行え.

concrete = pd.read_csv("http://logopt.com/data/concrete.csv")
concrete.head()
cement slag ash water splast coarse fine age strength
0 540.0 0.0 0.0 162.0 2.5 1040.0 676.0 28 79.986111
1 540.0 0.0 0.0 162.0 2.5 1055.0 676.0 28 61.887366
2 332.5 142.5 0.0 228.0 0.0 932.0 594.0 270 40.269535
3 332.5 142.5 0.0 228.0 0.0 932.0 594.0 365 41.052780
4 198.6 132.4 0.0 192.0 0.0 978.4 825.5 360 44.296075

問題(シェアバイク)

以下のシェアバイクのデータに対して,riders列が利用者数を,線形回帰を用いて推定せよ.ただし,date列とcasual列は除いてから回帰を行え. さらに結果の可視化を行え.

また,なぜcasual列を含めて推定をしないのか考察せよ.

bikeshare = pd.read_csv("http://logopt.com/data/bikeshare.csv")
bikeshare.head()
date season year month hour holiday weekday workingday weather temp feelslike humidity windspeed casual registered riders
0 2011-01-01 1 0 1 0 0 6 0 1 0.24 0.2879 0.81 0.0 3 13 16
1 2011-01-01 1 0 1 1 0 6 0 1 0.22 0.2727 0.80 0.0 8 32 40
2 2011-01-01 1 0 1 2 0 6 0 1 0.22 0.2727 0.80 0.0 5 27 32
3 2011-01-01 1 0 1 3 0 6 0 1 0.24 0.2879 0.75 0.0 3 10 13
4 2011-01-01 1 0 1 4 0 6 0 1 0.24 0.2879 0.75 0.0 0 1 1

ダミー変数を用いた回帰

ダイヤモンドの価格の予測(ダミー変数)

http://logopt.com/data/Diamond.csv からダイアモンドの価格データを読み込み,線形回帰による予測を行う.

列は ["carat","colour","clarity","certification","price"] であり,他の情報から価格(price)の予測を行え.

カラット(carat)以外の列は情報が文字列として保管されている.

これはカテゴリー変数とよばれ,sciki-learnで扱うには,数値に変換する必要がある.

pandasのget_dummies関数で数値情報(ダミー変数)に変換してから,線形回帰を行う.

たとえば,色を表すcolour列はD,E,F,G,H,Iの文字列が入っている.これを各値が入っているとき1,それ以外のとき0の数値に変換したものがダミー変数になる.

色はいずれかの値をとるので,ダミー変数は独立でない(1つが1になると,他のすべては0になる).

最初のダミー変数を除くには,get_dummies関数の引数のdrop_firstをTrueに設定すれば良い.

import pandas as pd

diamond = pd.read_csv("http://logopt.com/data/Diamond.csv", index_col=0)
diamond.head()
carat colour clarity certification price
1 0.30 D VS2 GIA 1302
2 0.30 E VS1 GIA 1510
3 0.30 G VVS1 GIA 1510
4 0.30 G VS1 GIA 1260
5 0.31 D VS1 GIA 1641
diamond = pd.get_dummies(diamond, drop_first=True)  # ダミー変数の最初のものを除く
# diamond = pd.get_dummies(diamond) # 除かなくても結果は同じ
diamond.head()
carat price colour_E colour_F colour_G colour_H colour_I clarity_VS1 clarity_VS2 clarity_VVS1 clarity_VVS2 certification_HRD certification_IGI
1 0.30 1302 0 0 0 0 0 0 1 0 0 0 0
2 0.30 1510 1 0 0 0 0 1 0 0 0 0 0
3 0.30 1510 0 0 1 0 0 0 0 1 0 0 0
4 0.30 1260 0 0 1 0 0 1 0 0 0 0 0
5 0.31 1641 0 0 0 0 0 1 0 0 0 0 0
y = diamond.price  # 従属変数(price)の抽出
X = diamond.drop("price", axis=1)  # 独立変数(特徴ベクトル)をpriceの列を除くことによって生成
X.head()
carat colour_E colour_F colour_G colour_H colour_I clarity_VS1 clarity_VS2 clarity_VVS1 clarity_VVS2 certification_HRD certification_IGI
1 0.30 0 0 0 0 0 0 1 0 0 0 0
2 0.30 1 0 0 0 0 1 0 0 0 0 0
3 0.30 0 0 1 0 0 0 0 1 0 0 0
4 0.30 0 0 1 0 0 1 0 0 0 0 0
5 0.31 0 0 0 0 0 1 0 0 0 0 0
from sklearn.linear_model import LinearRegression  # 線形回帰クラスのインポート

reg = LinearRegression()  # 線形回帰クラスのインスタンス生成
reg.fit(X, y)  # 訓練
yhat = reg.predict(X)  # 予測
print("y-切片= ", reg.intercept_)
print("係数 = ", reg.coef_)
print("決定変数= ", reg.score(X, y))  # 決定係数の別計算
y-切片=  169.17604383492744
係数 =  [12766.39597047 -1439.0853427  -1841.69054716 -2176.67218633
 -2747.14998002 -3313.1023993  -1474.56614749 -1792.01092358
  -689.29043537 -1191.16426364    15.22672874   141.2624469 ]
決定変数=  0.9581280577870392
visualizer = PredictionError(reg)

visualizer.fit(X, y)
visualizer.score(X, y)
visualizer.show()
<AxesSubplot:title={'center':'Prediction Error for LinearRegression'}, xlabel='$y$', ylabel='$\\hat{y}$'>
visualizer = ResidualsPlot(reg)

visualizer.fit(X, y)
visualizer.score(X, y)
visualizer.show()
<AxesSubplot:title={'center':'Residuals for LinearRegression Model'}, xlabel='Predicted Value', ylabel='Residuals'>

問題(車の価格)

http://logopt.com/data/carprice.csv から車の価格データを読み込み,線形回帰による予測を行え. また,結果を可視化せよ.

データの詳細は https://vincentarelbundock.github.io/Rdatasets/doc/DAAG/carprice.html にある.

車種(Type),100マイル走る際のガロン数(gpm100),都市部での1ガロンあたりの走行距離(MPGcity),高速道路での1ガロン当たりの走行距離(MPGhighway)から,価格(Price)を予測せよ.

carprice = pd.read_csv("http://logopt.com/data/carprice.csv", index_col=0)
carprice.head()
Type MinPrice Price MaxPrice RangePrice RoughRange gpm100 MPGcity MPGhighway
6 Midsize 14.2 15.7 17.3 3.1 3.09 3.8 22 31
7 Large 19.9 20.8 21.7 1.8 1.79 4.2 19 28
8 Large 22.6 23.7 24.9 2.3 2.31 4.9 16 25
9 Midsize 26.3 26.3 26.3 0.0 -0.01 4.3 19 27
10 Large 33.0 34.7 36.3 3.3 3.30 4.9 16 25

問題(チップ)

チップ(tips)データに対して線形回帰を用いてもらえるチップの額 (tip) を予測せよ. また,結果を可視化せよ.

import seaborn as sns

tips = sns.load_dataset("tips")
tips.head()
total_bill tip sex smoker day time size
0 16.99 1.01 Female No Sun Dinner 2
1 10.34 1.66 Male No Sun Dinner 3
2 21.01 3.50 Male No Sun Dinner 3
3 23.68 3.31 Male No Sun Dinner 2
4 24.59 3.61 Female No Sun Dinner 4

ロジスティック回帰による分類

ロジスティック回帰(logistic regression)は,その名前とは裏腹に分類のための手法である.

ここでは出力が $0$ もしくは $1$ の値をもつ分類問題を考える. ロジスティック回帰の仮説関数 $h_w (x)$ は,分類が $0$ になるのか $1$ になるのかを表す確率を表す. ロジスティック回帰の名前の由来は, 条件 $0 \leq h_w (x) \leq 1$ を満たすために, 線形回帰の仮説関数にロジスティック関数(logistic function; シグモイド関数)を用いることによる. ロジスティック関数 $g(z)$ は,以下のように定義される. (グラフを見たい場合には,Googleで 1/(1+e^-x)と検索する.)

$$ g(z)= \frac{1}{1+e^{-z}} $$

ここで $e$ は自然対数の底である.$z=0$ のときに $g(z)=0.5$ になり,$z \rightarrow \infty$ になると $g \rightarrow 1$,$z \rightarrow -\infty$ になると $g \rightarrow 0$ になることに注意されたい.

この関数を用いると,ロジスティック回帰の仮説関数は, $$ h_w (x)= g(w_0 +w_1 x_1 + w_2 x_2 + \cdots + w_n x_n ) $$ と表すことができる.

仮説関数を用いて実際の分類を行うためには,$h_w (x) \geq 0.5$ のときには $y=1$ と予測し, $h_w (x) < 0.5$ のときには $y=0$ と予測すれば良い. すなわち,$h_w (x) = 0.5$ を境界として, 領域を $y=1$ と $y=0$ の2つに分けるのである. このことから,$h_w (x) = 0.5$ は決定境界(decision boundary)とよばれる.

ロジスティック回帰における費用関数 は以下のように定義される.

$$ J(w)= \frac{1}{m} \displaystyle\sum_{i=1}^m Cost ( h_{w}(x^{(i)}),y^{(i)} ) $$$$ Cost (h_w(x),y)= \left\{ \begin{array}{l l } -\log (h_{w}(x) ) & y=1 \\ -\log (1- h_{w}(x) ) & y=0 \end{array} \right. $$

$y=1$ の場合には,この費用関数は $h_w(x)=1$ のとき $0$ になり, $h_w(x) \rightarrow 0$ になると無限大に近づく. $y=0$ の場合には,この費用関数は $h_w(x)=0$ のとき $0$ になり, $h_w(x) \rightarrow 1$ になると無限大に近づく. つまり,予測が当たる場合には費用が $0$ に,外れるにしたがって大きくなるように設定されている訳である. さらに,この費用関数は凸であることが示せるので,線形回帰と同様に非線形最適化の手法を用いることによって, 最適な重み $w$ を容易に求めることができるのである.

スパムの判定

メールがスパム(spam;迷惑メイル)か否かを判定する例題を用いる.

https://archive.ics.uci.edu/ml/datasets/spambase

様々な数値情報から,is_spam列が $1$ (スパム)か, $0$ (スパムでない)かを判定する.

import pandas as pd

spam = pd.read_csv("http://logopt.com/data/spam.csv")
spam.head()
word_freq_make word_freq_address word_freq_all word_freq_3d word_freq_our word_freq_over word_freq_remove word_freq_internet word_freq_order word_freq_mail ... char_freq_; char_freq_( char_freq_[ char_freq_! char_freq_$ char_freq_# capital_run_length_average capital_run_length_longest capital_run_length_total is_spam
0 0.21 0.28 0.50 0.0 0.14 0.28 0.21 0.07 0.00 0.94 ... 0.00 0.132 0.0 0.372 0.180 0.048 5.114 101 1028 1
1 0.06 0.00 0.71 0.0 1.23 0.19 0.19 0.12 0.64 0.25 ... 0.01 0.143 0.0 0.276 0.184 0.010 9.821 485 2259 1
2 0.00 0.00 0.00 0.0 0.63 0.00 0.31 0.63 0.31 0.63 ... 0.00 0.137 0.0 0.137 0.000 0.000 3.537 40 191 1
3 0.00 0.00 0.00 0.0 0.63 0.00 0.31 0.63 0.31 0.63 ... 0.00 0.135 0.0 0.135 0.000 0.000 3.537 40 191 1
4 0.00 0.00 0.00 0.0 1.85 0.00 0.00 1.85 0.00 0.00 ... 0.00 0.223 0.0 0.000 0.000 0.000 3.000 15 54 1

5 rows × 58 columns

is_spam列が従属変数(ターゲット)$y$ になり,それ以外の列が独立変数(特徴ベクトル)$X$ になる

ロジスティック回帰を用いる。

X = spam.drop("is_spam", axis=1)
y = spam.is_spam
X
word_freq_make word_freq_address word_freq_all word_freq_3d word_freq_our word_freq_over word_freq_remove word_freq_internet word_freq_order word_freq_mail ... word_freq_conference char_freq_; char_freq_( char_freq_[ char_freq_! char_freq_$ char_freq_# capital_run_length_average capital_run_length_longest capital_run_length_total
0 0.21 0.28 0.50 0.0 0.14 0.28 0.21 0.07 0.00 0.94 ... 0.0 0.000 0.132 0.0 0.372 0.180 0.048 5.114 101 1028
1 0.06 0.00 0.71 0.0 1.23 0.19 0.19 0.12 0.64 0.25 ... 0.0 0.010 0.143 0.0 0.276 0.184 0.010 9.821 485 2259
2 0.00 0.00 0.00 0.0 0.63 0.00 0.31 0.63 0.31 0.63 ... 0.0 0.000 0.137 0.0 0.137 0.000 0.000 3.537 40 191
3 0.00 0.00 0.00 0.0 0.63 0.00 0.31 0.63 0.31 0.63 ... 0.0 0.000 0.135 0.0 0.135 0.000 0.000 3.537 40 191
4 0.00 0.00 0.00 0.0 1.85 0.00 0.00 1.85 0.00 0.00 ... 0.0 0.000 0.223 0.0 0.000 0.000 0.000 3.000 15 54
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
4595 0.31 0.00 0.62 0.0 0.00 0.31 0.00 0.00 0.00 0.00 ... 0.0 0.000 0.232 0.0 0.000 0.000 0.000 1.142 3 88
4596 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00 0.00 ... 0.0 0.000 0.000 0.0 0.353 0.000 0.000 1.555 4 14
4597 0.30 0.00 0.30 0.0 0.00 0.00 0.00 0.00 0.00 0.00 ... 0.0 0.102 0.718 0.0 0.000 0.000 0.000 1.404 6 118
4598 0.96 0.00 0.00 0.0 0.32 0.00 0.00 0.00 0.00 0.00 ... 0.0 0.000 0.057 0.0 0.000 0.000 0.000 1.147 5 78
4599 0.00 0.00 0.65 0.0 0.00 0.00 0.00 0.00 0.00 0.00 ... 0.0 0.000 0.000 0.0 0.125 0.000 0.000 1.250 5 40

4600 rows × 57 columns

LogisticRegressionクラスの引数

最適化の引数をsolver引数で変えることができる. 引数には以下のものが準備されている.

  • ‘newton-cg’: Newton共役方向法
  • ‘lbfgs’: limited memory BFGS法 (既定値)
  • ‘liblinear’
  • ‘sag’
  • ‘saga’
from sklearn.linear_model import LogisticRegression  # ロジスティック回帰クラスの読み込み

logreg = LogisticRegression()  # インスタンスの生成
logreg.fit(X, y)  # 訓練
LogisticRegression()

分類に対する可視化

  • 混合行列 (ConfusionMatrix)
  • 分類レポート (ClassificationReport)
  • 2値分類に対する閾値変化図 (DiscriminationThreshold)
  • ROC曲線 (ROCAUC)

混合行列

  • positive(陽性) : 予測(prediction)が 1(スパム);(注意)positiveか否かは相対的なものであるので,どちらでも良い.0(スパムでない)をpositiveとすることも可能.
  • negative(陰性): 予測が0(スパムでない)

  • true (真): あたり

  • false(偽): はずれ

混合行列

行を正解,列を予測とし,Trueが対角線になるようにした行列(普通は1を左上にするが,yellowbrickだと0,1の順に並ぶので注意)

True\Prediction not spam is spam
not spam (0) TN FP
is spam (1) FN TP

正解率

$$ accuracy = (TP+TN)/(TP+FN+FP+TN) $$
from yellowbrick.classifier import ConfusionMatrix

cm = ConfusionMatrix(logreg, classes=["not spam", "is spam"])

cm.fit(X, y)
cm.score(X, y)
cm.show();

メトリクス(評価尺度)

spamでないメイルをspamと判断するのは(false positive),大事なメイルがスパムフォルダに入ってしまうので困る.

spamをspamでないと判断するのは(false negative),スパムを消せば良いのであまり困らない.

したがって,他の尺度 (recall, precision, f1 score) を計算する.

true positive rate = recall = TP/(TP+FN)

false positive rate = 1-precision = FP/(TN+FP)

F1 score = 2 * ((precision * recall) / (precision + recall))
TP, FN, FP, TN = 2665, 123, 190, 1622
precision = TP / (TP + FP)
recall = TP / (FN + TP)
print("precision=", precision)
print("recall=", recall)
print("F1 score=", 2 * (precision * recall) / (precision + recall))
precision= 0.9334500875656743
recall= 0.9558823529411765
F1 score= 0.9445330497962078
from yellowbrick.classifier import ClassificationReport

visualizer = ClassificationReport(logreg, classes=["not spam", "is spam"])

visualizer.fit(X, y)
visualizer.score(X, y)
visualizer.show();

閾値 (threshold)を変えてみる

ロジスティック回帰は,spamでない確率を推定し,それが閾値より大きいとspam,それ以外のときspamでないと判定する.

通常のロジスティック回帰の閾値は0.5

閾値を大きくすると,全部spamでない(negative)と判定(false positiveが減り,precisionが上がる)

閾値を小さくすると,全部spam(positive)と判定(false negativeが減り,recall が上がる)

f1 scoreはバランスをとる.

from yellowbrick.classifier import DiscriminationThreshold

visualizer = DiscriminationThreshold(logreg)

visualizer.fit(X, y)
visualizer.score(X, y)
visualizer.show();

ROC曲線

ROC: 受信者操作特性(receiver operating characteristic)

https://ja.wikipedia.org/wiki/受信者操作特性

閾値を1から0に変えると,recall (true positive rate) が増加し,precision (1-false positive rate) が減少する.

x軸に false positive rate,y軸にtrue positive rateをとると,両方とも1に近づく曲線になる.

曲線の下の面積 (AUC: area under the curve) が評価尺度(大きいほど良い)になる.

true positive rate = recall = TP/(TP+FN)

false positive rate = 1-precision = FP/(TN+FP)
from yellowbrick.classifier import ROCAUC

visualizer = ROCAUC(logreg, size=(600, 400))

visualizer.fit(X, y)
visualizer.score(X, y)
visualizer.show();

問題(クレジットカード)

以下のクレジットカードのデフォルトの判定データに対してロジスティック回帰を行い, 混合行列を描画せよ.

default列にデフォルトか否かの情報があり,他の列の情報を用いて分類せよ.

credit = pd.read_csv("http://logopt.com/data/credit.csv")
credit.head()
limit sex edu married age apr_delay may_delay jun_delay jul_delay aug_delay ... jul_bill aug_bill sep_bill apr_pay may_pay jun_pay jul_pay aug_pay sep_pay default
0 20000 2 2 1 24 2 2 -1 -1 -2 ... 0 0 0 0 689 0 0 0 0 1
1 120000 2 2 2 26 -1 2 0 0 0 ... 3272 3455 3261 0 1000 1000 1000 0 2000 1
2 90000 2 2 2 34 0 0 0 0 0 ... 14331 14948 15549 1518 1500 1000 1000 1000 5000 0
3 50000 2 2 1 37 0 0 0 0 0 ... 28314 28959 29547 2000 2019 1200 1100 1069 1000 0
4 50000 1 2 1 57 -1 0 -1 0 0 ... 20940 19146 19131 2000 36681 10000 9000 689 679 0

5 rows × 24 columns

問題(部屋)

以下の部屋が使われているか否かを判定するデータに対してロジスティック回帰を行い,混合行列とROC曲線を描画せよ.

occupancy列が部屋が使われているか否かを表す情報であり,これをdatetime列以外の情報から分類せよ.

occupancy = pd.read_csv("http://logopt.com/data/occupancy.csv")
occupancy.head()
datetime temperature relative humidity light CO2 humidity occupancy
0 2015-02-04 17:51:00 23.18 27.2720 426.0 721.25 0.004793 1
1 2015-02-04 17:51:59 23.15 27.2675 429.5 714.00 0.004783 1
2 2015-02-04 17:53:00 23.15 27.2450 426.0 713.50 0.004779 1
3 2015-02-04 17:54:00 23.15 27.2000 426.0 708.25 0.004772 1
4 2015-02-04 17:55:00 23.10 27.2000 426.0 704.50 0.004757 1

多クラス分類

上ではロジスティック回帰を用いた2値分類について考えた. ここでは,iris(アヤメ)のデータを用いて,3種類以上の分類(これを多クラス分類とよぶ)を解説する.

アヤメの例題では,4つの花の特徴を用いて3つの種類に分類する.

入力は長さ4のベクトル $x$ であり,$4 \times 3$ の重み行列 $W$ に $x$ を乗じることによって, 長さ3のベクトル $z =[z_1,z_2,z_3]$ を得る.重み $W$ は分類がうまくいくように決めるパラメータである. 

$z_1, z_2, z_3$ の各値を用いて, 各クラス(アヤメの種類)が選択される確率を計算する.これには,以下のソフトマックス関数を用いる. $$ p_j = \frac{e^{z_j}}{\sum_{k=1}^3 e^{z_k}} \ \ \ j=1,2,3 $$ ここで $e$ は自然対数の底である.

分類は $p_j$ が最も大きい $j$ を選択することによって行われる. これをソフトマックス回帰とよぶ.ただしscikit-learnでは,ロジスティック回帰と同じクラスを用いる. ソフトマックス回帰における費用関数は,以下のように定義される.

$$ J(w)= \frac{1}{m} \displaystyle\sum_{i=1}^m Cost ( h_{w}(x^{(i)}),y^{(i)} ) $$$$ Cost (h_w(x),y)= - \sum_{j=1}^3 [y=j] \log \frac{e^{z_j}}{\sum_{k=1}^3 e^{z_k}} $$

ここで,$[y=j]$ はデータ $y$ の正解が $j$ のとき $1$,そうでないとき $0$ の記号である.

以下ではアヤメのデータを用いてロジスティック(ソフトマックス)回帰を行い,3種類のアヤメを分類する.

import plotly.express as px

iris = px.data.iris()
iris.head()
sepal_length sepal_width petal_length petal_width species species_id
0 5.1 3.5 1.4 0.2 setosa 1
1 4.9 3.0 1.4 0.2 setosa 1
2 4.7 3.2 1.3 0.2 setosa 1
3 4.6 3.1 1.5 0.2 setosa 1
4 5.0 3.6 1.4 0.2 setosa 1
X = iris[["sepal_length", "sepal_width", "petal_length", "petal_width"]]
# 従属変数 y
y = iris["species_id"]
y.head()
0    1
1    1
2    1
3    1
4    1
Name: species_id, dtype: int64
  • 手順1:分類するためのクラスをインポートして,インスタンスを生成する.

  • 手順2:fitメソッドを用いて,訓練する.

  • 手順3:predictメソッドを用いて予測を行う.

from sklearn.linear_model import LogisticRegression  # ロジスティック回帰クラスの読み込み

logreg = LogisticRegression()  # インスタンスの生成
logreg.fit(X, y)  # 訓練
LogisticRegression()
logreg.predict([[3, 5, 4, 2]])  # 試しに予測
array([1])

予測と実際の誤差を検証

元データXを入れたときの予測y_predと本当の値yの比較

metricsにある正解率を計算する関数 accuracy_score を利用

y_pred = logreg.predict(X)
from sklearn import metrics

print(metrics.accuracy_score(y, y_pred))
0.9733333333333334

可視化

以下の2つの可視化を追加

  • どの特徴が重要かを可視化するRadViz

  • 重要と思われる2つの特徴に対して,ロジスティック回帰の分類がどのように行われたかを示すDecisionViz(決定の境界を示す)

スパムと同じ可視化

  • 混合行列 (ConfusionMatrix)
  • 分類レポート (ClassificationReport)
  • ROC曲線 (ROCAUC)

3値以上の場合には,閾値を変化させる可視化はできないことに注意

from yellowbrick.features import RadViz

visualizer = RadViz(
    classes=["Iris-setosa", "Iris-versicolor", "Iris-virginica"], size=(600, 400)
)

visualizer.fit(X, y)
visualizer.transform(X)
visualizer.show();
from sklearn.preprocessing import StandardScaler
from yellowbrick.contrib.classifier import DecisionViz

X = iris[["sepal_width", "petal_width"]]  # 2次元を切り出す
X = StandardScaler().fit_transform(X)  # 可視化のためにスケーリングしておく

viz = DecisionViz(
    logreg,
    features=["sepal_width", "petal_width"],
    classes=["Iris-setosa", "Iris-versicolor", "Iris-virginica"],
)
viz.fit(X, y)
viz.draw(X, y)
viz.show();
/Users/mikiokubo/Library/Caches/pypoetry/virtualenvs/analytics-v-sH3Dza-py3.8/lib/python3.8/site-packages/yellowbrick/contrib/classifier/boundaries.py:435: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3.  Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading'].  This will become an error two minor releases later.
  self.ax.pcolormesh(
from yellowbrick.classifier import ClassificationReport

visualizer = ClassificationReport(logreg)

visualizer.fit(X, y)
visualizer.score(X, y)
visualizer.show();
from yellowbrick.classifier import ConfusionMatrix

cm = ConfusionMatrix(logreg)

cm.fit(X, y)
cm.score(X, y)
cm.show();
from yellowbrick.classifier import ROCAUC

visualizer = ROCAUC(logreg, size=(600, 400))

visualizer.fit(X, y)
visualizer.score(X, y)
visualizer.show();

毒キノコの判定

ここでは,データから毒キノコか否かを判定する.

target列がターゲット(従属変数)であり,edibleが食用,poisonousが毒である.

他の列のデータもすべて数値ではない.ここでは,scikit-learnのprocessingを用いて,数値に変換してからロジスティック回帰により分類を行う.

mashroom = pd.read_csv("http://logopt.com/data/mashroom.csv")
mashroom.head()
target shape surface color
0 edible convex smooth yellow
1 edible bell smooth white
2 poisonous convex scaly white
3 edible convex smooth gray
4 edible convex scaly yellow
mashroom.color.unique()
array(['yellow', 'white', 'gray', 'brown', 'red', 'pink', 'buff',
       'purple', 'cinnamon', 'green'], dtype=object)
X = mashroom.drop("target", axis=1)
y = mashroom.target
X
shape surface color
0 convex smooth yellow
1 bell smooth white
2 convex scaly white
3 convex smooth gray
4 convex scaly yellow
... ... ... ...
8118 knobbed smooth brown
8119 convex smooth brown
8120 flat smooth brown
8121 knobbed scaly brown
8122 convex smooth brown

8123 rows × 3 columns

前処理(preprocessing)

sklearn.preprocessingモジュールにあるLabelEncoderクラスのfit_transformメソッドを用いてターゲット列を数値化する.

edibleは0に,poisonousは1に変換されていることが分かる.

独立変数 $X$ に対しては,OrdinalEncoderを用いる.これによって文字列データが数値に変換されていることが分かる.

なお,OneHotEncoderを使うと,ダイヤモンドの例題でダミー変数に変換したのと同じことができる.

from sklearn.preprocessing import LabelEncoder

y = LabelEncoder().fit_transform(y)
y
array([0, 0, 1, ..., 0, 1, 0])
from sklearn.preprocessing import OrdinalEncoder, OneHotEncoder
#X = OneHotEncoder().fit_transform(X)
#X = OrdinalEncoder().fit_transform(X)
#print(X)
from sklearn.linear_model import LogisticRegression

logreg = LogisticRegression()
logreg.fit(X, y)
LogisticRegression()
from yellowbrick.classifier import ConfusionMatrix

cm = ConfusionMatrix(logreg, classes=["edible", "poisonous"])

cm.fit(X, y)
cm.score(X, y)
cm.show();
from yellowbrick.classifier import ClassificationReport

visualizer = ClassificationReport(logreg)

visualizer.fit(X, y)
visualizer.score(X, y)
visualizer.show();
from yellowbrick.classifier import DiscriminationThreshold

visualizer = DiscriminationThreshold(logreg)

visualizer.fit(X, y)
visualizer.score(X, y)
visualizer.show();
from yellowbrick.classifier import ROCAUC

visualizer = ROCAUC(logreg, size=(600, 400))

visualizer.fit(X, y)
visualizer.score(X, y)
visualizer.show();