可視化パッケージYellowbrickのインストール
以下を参照
https://www.scikit-yb.org/en/latest/
Google Colab.上で古いバージョンがインストールされている場合には，以下のコードを実行後に， RUNTIMEのリスタートが必要である．
#!pip install -U yellowbrick
import yellowbrick
yellowbrick.__version__
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()
独立変数（特徴ベクトル）$X$ は TV, Radio, Newspaperの列，従属変数（ターゲット） $y$ は Salesの列
y = df["Sales"]
X = df[["TV", "Radio", "Newspaper"]] # df.drop("Sales",axis=1) でも同じ
X.head()
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_)
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
print(reg.score(X, y)) # 決定係数の別計算
from yellowbrick.regressor import PredictionError
visualizer = PredictionError(reg)
visualizer.fit(X, y)
visualizer.score(X, y)
visualizer.show()
from yellowbrick.regressor import ResidualsPlot
visualizer = ResidualsPlot(reg)
visualizer.fit(X, y)
visualizer.score(X, y)
visualizer.show()
問題 (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()
concrete = pd.read_csv("http://logopt.com/data/concrete.csv")
concrete.head()
bikeshare = pd.read_csv("http://logopt.com/data/bikeshare.csv")
bikeshare.head()
ダミー変数を用いた回帰
ダイヤモンドの価格の予測（ダミー変数）
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()
diamond = pd.get_dummies(diamond, drop_first=True) # ダミー変数の最初のものを除く
# diamond = pd.get_dummies(diamond) # 除かなくても結果は同じ
diamond.head()
y = diamond.price # 従属変数(price)の抽出
X = diamond.drop("price", axis=1) # 独立変数（特徴ベクトル）をpriceの列を除くことによって生成
X.head()
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)) # 決定係数の別計算
visualizer = PredictionError(reg)
visualizer.fit(X, y)
visualizer.score(X, y)
visualizer.show()
visualizer = ResidualsPlot(reg)
visualizer.fit(X, y)
visualizer.score(X, y)
visualizer.show()
問題（車の価格）
http://logopt.com/data/carprice.csv から車の価格データを読み込み，線形回帰による予測を行え． また，結果を可視化せよ．
データの詳細は https://vincentarelbundock.github.io/Rdatasets/doc/DAAG/carprice.html にある．
車種(Type)，100マイル走る際のガロン数（gpm100），都市部での1ガロンあたりの走行距離（MPGcity），高速道路での１ガロン当たりの走行距離（MPGhighway）から，価格(Price)を予測せよ．
carprice = pd.read_csv("http://logopt.com/data/carprice.csv", index_col=0)
carprice.head()
import seaborn as sns
tips = sns.load_dataset("tips")
tips.head()
ロジスティック回帰による分類
ロジスティック回帰(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()
is_spam列が従属変数（ターゲット）$y$ になり，それ以外の列が独立変数（特徴ベクトル）$X$ になる
ロジスティック回帰を用いる。
X = spam.drop("is_spam", axis=1)
y = spam.is_spam
X
from sklearn.linear_model import LogisticRegression # ロジスティック回帰クラスの読み込み
logreg = LogisticRegression() # インスタンスの生成
logreg.fit(X, y) # 訓練
混合行列
- positive（陽性） : 予測(prediction)が １（スパム）；(注意）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))
from yellowbrick.classifier import ClassificationReport
visualizer = ClassificationReport(logreg, classes=["not spam", "is spam"])
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();
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();
credit = pd.read_csv("http://logopt.com/data/credit.csv")
credit.head()
occupancy = pd.read_csv("http://logopt.com/data/occupancy.csv")
occupancy.head()
多クラス分類
上ではロジスティック回帰を用いた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()
X = iris[["sepal_length", "sepal_width", "petal_length", "petal_width"]]
# 従属変数 y
y = iris["species_id"]
y.head()
手順１：分類するためのクラスをインポートして，インスタンスを生成する．
手順２：fitメソッドを用いて，訓練する．
手順３：predictメソッドを用いて予測を行う．
from sklearn.linear_model import LogisticRegression # ロジスティック回帰クラスの読み込み
logreg = LogisticRegression() # インスタンスの生成
logreg.fit(X, y) # 訓練
logreg.predict([[3, 5, 4, 2]]) # 試しに予測
y_pred = logreg.predict(X)
from sklearn import metrics
print(metrics.accuracy_score(y, y_pred))
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();
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();
mashroom = pd.read_csv("http://logopt.com/data/mashroom.csv")
mashroom.head()
mashroom.color.unique()
X = mashroom.drop("target", axis=1)
y = mashroom.target
X
from sklearn.preprocessing import LabelEncoder
y = LabelEncoder().fit_transform(y)
y
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)
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();