機械学習とは
機械学習(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~(=1,2,\ldots,m)$ に対して，仮説関数によって「予測」された値 $\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モジュールを準備する．
# csvファイルからデータ読み込み
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
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に設定すれば良い．
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$ を容易に求めることができるのである．
2値分類 （スパム判定）
メイルがスパム（spam；迷惑メイル）か否かを判定する例題を用いる．
https://archive.ics.uci.edu/ml/datasets/spambase
様々な数値情報から，is_spam列が $1$ （スパム）か， $0$ （スパムでない）かを判定する．
spam = pd.read_csv("http://logopt.com/data/spam.csv")
spam.head().T
is_spam列が従属変数（ターゲット）$y$ になり，それ以外の列が独立変数（特徴ベクトル）$X$ になる．
ロジスティック回帰を用いる。
X = spam.drop("is_spam", axis=1)
y = spam.is_spam
LogisticRegressionクラスの引数
最適化をするための手法をsolver引数で変えることができる． 以下のものが準備されている．
- "newton-cg": Newton共役方向法
- "lbfgs": 記憶制限付きBFGS(Broyden–Fletcher–Goldfarb–Shanno)法 （既定値）
- "liblinear": 座標降下法（小さなデータのとき有効）
- "sag" : 確率的平均勾配法 (Stochastic Average Gradient)； 大規模なデータで高速
- "saga" : 確率的平均勾配法の変形； 大規模なデータで高速
- "newton-cholesky": データ（サンプル）数が特徴の数と比較して非常に大きいときに有効
from sklearn.linear_model import LogisticRegression # ロジスティック回帰クラスの読み込み
logreg = LogisticRegression(solver="liblinear") # インスタンスの生成
logreg.fit(X, y)
# 訓練
混同行列
混同行列 (confusion matrix)とは， 行を正解，列を予測とし，あたり（true）が対角線になるようにした行列である． 普通は $1$ を左上にするが，yellowbrickだと $0,1$ の順に並ぶので注意する必要がある．
| True\Prediction | not spam | is spam |
|---|---|---|
| not spam (0) | TN | FP |
| is spam (1) | FN | TP |
- positive（陽性）: 予測(prediction)が $１$ （スパム）；(注意）positiveか否かは相対的なものであるので，どちらでも良い． $0$ （スパムでない）をpositiveとすることも可能である．
negative（陰性）: 予測が $0$ （スパムでない）
true （真）: あたり
false（偽）: はずれ
正解率 (accuracy) は以下のように定義される評価尺度である．
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など) が準備されている．
$$ \mathrm{true~positive~rate} = \mathrm{recall} = \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}} $$$$ \mathrm{precision} = \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}} $$$$ \mathrm{false~positive~rate} = \frac{\mathrm{FP}}{\mathrm{TN}+\mathrm{FP}} $$$$ \mathrm{f1~score} = \frac{2 (\mathrm{precision} \cdot \mathrm{recall}) }{\mathrm{precision} + \mathrm{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()
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$ を乗じることによって， ベクトル $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()
titanic = pd.read_csv("http://logopt.com/data/titanic.csv")
titanic.head()
問題（胸部癌）
http://logopt.com/data/cancer.csv にある胸部癌か否かを判定するデータセットを用いて分類を行え．
最初の列diagnosisが癌か否かを表すものであり，"M"が悪性（malignant），"B"が良性（benign）を表す．
cancer = pd.read_csv("http://logopt.com/data/cancer.csv", index_col=0)
cancer.head().T # 横長なので転置して表示する．
$K$ 近傍法による分類と可視化
4種類のフルーツを重量と色で分類する．
フルーツは数値だが，順に "apple", "mandarin","orange","lemon" を表す．
$K$ 近傍法とは，($X$の空間で）データに近い $K$ 個のデータの値($y$)の多数決で分類を行う最も簡単な分類手法である．
以下のようにしてクラスをインポートできる．
from sklearn.neighbors import KNeighborsClassifier
クラス KNeighborsClassifier は引数 n_neighbors でパラメータ $K$ を設定できる．
$K$ を色々変えて実験せよ（以下では$K=3$の結果を示す）．
fruit = pd.read_csv("http://logopt.com/data/fruit_simple.csv")
fruit.head()
from sklearn.preprocessing import StandardScaler
from sklearn.neighbors import KNeighborsClassifier
from yellowbrick.contrib.classifier import DecisionViz
X = fruit.drop("fruit_label", axis=1)
y = fruit.fruit_label
X = StandardScaler().fit_transform(X) # 可視化のためにスケーリングしておく
knn = KNeighborsClassifier(3)
knn.fit(X, y) # 訓練
y_pred = knn.predict(X)
print("正解率=", metrics.accuracy_score(y, y_pred))
viz = DecisionViz(
knn,
title="Nearest Neighbors",
features=["mass", "color score"],
classes=["apple", "mandarin", "orange", "lemon"],
)
viz.fit(X, y)
viz.draw(X, y)
viz.show()
多項式回帰
回帰や分類を行う際に予測誤差が小さい方が望ましいことは言うまでもない． しかし，予測誤差を小さくするために変数（特性）を増やしすぎると， トレーニングデータに適合しすぎて，実際のデータに対して良い結果を出さない 過剰適合(overfitting; 過学習)が発生してしまう．
たとえば，線形回帰の仮説関数として以下のような $x$ に対する多項式関数を考える． $$ h_w (x)= w_0 +w_1 x_1 + w_2 x_2 + \cdots + w_n x_n + w_{n+1} x_1^2 + w_{n+2} x_2^2 + \cdots + w_{2n} x_n^2 + w_{2n+1} x_1^3 + w_{2n+2} x_2^3 + \cdots + w_{3n} x_n^3 $$ この関数は重み $w$ に対する線形関数であるので，（特性の数は $3$倍になるが） 線形回帰と同じように最適な重みを計算することができる． ちなみに，このような回帰を多項式回帰(polynomial regression)とよぶ． 一般に高次の多項式回帰でトレーニングデータの数が少ない場合には，過剰適合が発生する．
宣伝の効果のデータを用いて多項式回帰を行う．
pandasで2次の項（たとえばテレビとラジオの相乗効果の列なら TV*Radio）を生成してから，線形回帰を行う．
data = pd.read_csv("http://logopt.com/data/Advertising.csv", index_col=0)
data["TV*Radio"] = data.TV * data.Radio
data["TV*Newspaper"] = data.TV * data.Newspaper
data["Radio*Newspaper"] = data.Radio * data.Newspaper
data.head()
X = data[["TV", "Radio", "Newspaper", "TV*Radio", "TV*Newspaper", "Radio*Newspaper"]]
y = data["Sales"]
from sklearn.linear_model import LinearRegression # クラスのインポート
lin_reg = LinearRegression() # 線形回帰クラスのインスタンス生成
lin_reg.fit(X, y) # 訓練
yhat = lin_reg.predict(X) # 予測
print(lin_reg.score(X, y)) # 決定係数 （線形回帰だと R^2は 0.897210638179
print("y-切片= ", lin_reg.intercept_)
print("係数 = ", lin_reg.coef_)
例としてBostonの住宅データの予測を行う．
medvが住宅の価格で，他のデータ（犯罪率や人口など）から予測する．
boston = pd.read_csv("http://logopt.com/data/Boston.csv", index_col=0)
boston.head()
X = boston.drop("medv", axis=1) # 最後の列以外のデータを独立変数（特徴ベクトル）として抽出
y = boston.medv # 最後の列(medv)を従属変数として抽出
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.3
) # 30%のデータをテスト用，それ以外をトレーニング用に分離
from sklearn.linear_model import LinearRegression # クラスのインポート
reg = LinearRegression() # 線形回帰クラスのインスタンス生成
reg.fit(X_train, y_train) # 訓練
yhat = reg.predict(X_test) # 予測
print(reg.score(X_test, y_test)) # 決定係数
残差プロットを描画すると，トレーニングデータとテストデータに色分けして表示される．
from yellowbrick.regressor import ResidualsPlot
visualizer = ResidualsPlot(reg)
visualizer.fit(X_train, y_train)
visualizer.score(X_test, y_test)
visualizer.show()
リッジ回帰
過剰適合を避けるための方法として正則化 (regulalization) がある。
正則化を追加した線形回帰は，リッジ回帰(ridge regression)とよばれる．
正則化のアイディアは単純であり，重み $w$ を小さくするような項を費用関数に追加するだけである．
線形回帰の場合には，以下の費用関数を用いる． $$ J(w)= \frac{1}{2m} \left[ \sum_{i=1}^m \left( h_{w}(x^{(i)}) -y^{(i)} \right)^2 + \lambda \sum_{j=1}^n w_j^2 \right] $$
追加された2項目が重みを小さくし，これによって過剰適合が軽減される． パラメータ $\lambda$ は正則化パラメータ(regularization parameter)とよばれ， 誤差と過剰適合のトレードオフをとるために用いられる． $\lambda$ が小さいときには，誤差は少なくなるが過剰適合の危険性が増大する． 逆に $\lambda$ を大きくすると，誤差は大きくなり，過剰適合は回避される傾向になる．
パラメータ $\lambda$ はリッジ回帰クラス Ridge では引数 alpha で与える．
from sklearn.linear_model import Ridge # クラスのインポート
reg = Ridge(alpha=10.0) # リッジ回帰クラスのインスタンス生成
reg.fit(X_train, y_train) # 訓練
yhat = reg.predict(X_test) # 予測
print(reg.score(X_test, y_test)) # 決定係数
最も良い正則化パラメータalphaを求める．
後述する交差検証を使い，yellowbrickのAlphaSelectionで可視化する．
from sklearn.linear_model import RidgeCV
from yellowbrick.regressor import AlphaSelection
import numpy as np
alphas = np.logspace(-100, -1, 400)
model = RidgeCV(alphas=alphas)
visualizer = AlphaSelection(model)
visualizer.fit(X, y)
visualizer.show()
reg = Ridge(alpha=0.032) # リッジ回帰クラスのインスタンス生成
reg.fit(X_train, y_train) # 訓練
yhat = reg.predict(X_test) # 予測
print(reg.score(X_test, y_test)) # 決定係数
問題（アヤメ，タイタニック，ダイヤモンド）
- irisデータをトレーニング用とテスト用に分けてから，ロジスティック回帰による分類を行い，正解率を計算せよ．
- titanicデータトレーニング用とテスト用に分けてから，ロジスティック回帰による分類を行い，正解率を計算せよ．
- 例題2で用いた ダイアモンドの価格データhttp://logopt.com/data/Diamond.csv に対して線形回帰とリッジ回帰による予測を行え．また，トレーニング用とテスト用に分けて決定係数 $R^2$ を計算し，評価せよ． ただしリッジ回帰の正則化パラメータは $0.03$ と設定せよ．
カーネルとSVM (Support Vector Machine）
分類のための手法の1つである SVM(Support Vector Machine；サポートベクトルマシン) について述べる．
SVMはロジスティック回帰を計算しやすいように近似したものと考えることができる．
ロジスティック回帰における費用関数 $$ 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. $$ は $$ Cost (h_w(x),y)= -y \log (h_{w}(x) ) - (1-y) \log (1- h_{w}(x) ) $$ と書き直すことができる．
ロジスティック回帰の仮説関数 $h_w (x)$ は， $z=w_0 +w_1 x_1 + w_2 x_2 + \cdots + w_n x_n$ とおいたときに $$ h_w (x)= g(z) =\frac{1}{1+e^{-z}} $$ であったことを思い起こすと，費用関数は $$ Cost (h_w(x),y)= -y \log \frac{1}{1+e^{-z}} - (1-y) \log (1- \frac{1}{1+e^{-z}} ) $$ となる．
$y=1$ のときの費用関数は第1項目だけになるので， $$ -\log \frac{1}{1+e^{-z}} $$ となる．これは $z$ が大きくなると $0$ に近づく曲線である （グラフを見たい場合には，Googleで-log(1/(1+e^-x))と検索する）． これを $z \geq 1$ のとき $0$ となる区分的線形関数 $cost_1(z)$ で近似する．
$y=0$ のときも同様に第2項 $$ -\log (1- \frac{1}{1+e^{-z}} ) $$ は，$z$ が小さくなると $0$ に近づく曲線となるので， $z \leq -1$ のときに $0$ になる区分的線形関数 $cost_0(z)$ で近似する． これに正則化パラメータを加えたものがSVMの費用関数 $J(w)$ になる． $$ J(w)= C \displaystyle\sum_{i=1}^m \left[ y^{(i)} cost_1(w^T x^{(i)}) + (1-y^{(i)}) cost_0(w^T x^{(i)}) \right] +\frac{1}{2} \displaystyle\sum_{j=1}^m w_j^2 $$ ここで $C$ は正則化パラメータを $\lambda$ としたとき $1/2 \lambda$ に相当するパラメータである． $\lambda$ のかわりに $C$ を用いるのは歴史的な理由に基づく慣例であり，特に意味はない．
パラメータ $C$ を大きく設定することは正則化パラメータ $\lambda$ を小さく設定することに相当するので， 誤差は小さくなるが過剰適合の危険性が高まり，$C$ を小さくするとその逆になる．
仮説関数を用いて実際の分類を行うためには， $w_0 +w_1 x_1 + w_2 x_2 + \cdots + w_n x_n \geq 0$ のときには $y=1$ と予測し， それ以外のときには $y=0$ と予測すれば良い．
費用関数 $J(w)$ の第1項が $0$ の場合には，$y^{(i)}=1$ なら $w^T x^{(i)} \geq 1$ かつ $y^{(i)}=0$ なら $w^T x^{(i)} \leq -1$ になる必要がある． このことから，データが直線で2つのクラスに分類可能な場合には， SVMの決定境界は，$1,-1$ の幅をもたせた境界になることが言える．
次に，カーネル(kernel)の概念を用いたSVMを解説する． これは，決定境界が非線形になる場合にも適用可能な手法である． もちろん，$x$ に対する多項式を作成してからロジスティック回帰 を適用しても良いのだが，非線形性が強い場合には高次の多項式が必要になり， それは計算量の増大をもたらす． 以下では，カーネルの中でも非線形な場合に適用可能で計算量的にも優秀な Gaussカーネル関数(Gaussian kernel function)を紹介する．
トレーニングデータの空間にランドマークとよばれる点列 $\ell^{(1)}, \ell^{(2)},\ldots,\ell^{(k)}$ が与えられているものとする（設定法については後述する）． このとき特性ベクトル $x=(x_1, x_2, \ldots, x_n)$ とランドマーク $\ell^{(i)}$ の「近さ」を表す尺度として 以下のGaussカーネル関数を用いる． $$ f_i = \exp \left( - \frac{ \sum_{j=1}^n (x_j-\ell_j^{(i)})^2 }{2 \sigma^2} \right) = \exp \left( - \frac{ \| x-\ell^{(i)} \|}{2 \sigma^2} \right) $$ ここで $\| \cdot \|$ はベクトルのノルムを表す． この関数は $x$ がランドマーク $\ell^{(i)}$ に近づくと $1$ になり， 遠ざかると $0$ になる多次元正規(Gauss)分布のような形状をとる． これがGaussカーネルとよばれる所以である． $\sigma$ は正規分布の標準偏差に対応するパラメータであり， $\sigma$ が小さくなると分布が $\ell^{(i)}$ の周りに集中し， 逆に$\sigma$ が大きくなると裾野が広がる．
仮説関数は，カーネルを用いて計算された $f_i$ と重み $w$ の線形関数として，以下のように定義する． $$ h_w (x)= w_0 +w_1 f_1 + w_2 f_2 + \cdots + w_k f_k $$ つまり，元の変数 $x$ をそのまま特性ベクトルとして用いるのではなく， $f= (f_1,f_2,\ldots,f_k)$ を特性ベクトルとして予測を行う訳である． これは多次元正規分布を重ね合わせたような分布になるので，非線形な決定境界を うまく表現できることが期待できる．
ランドマーク $\ell^{(i)}$ としては，トレーニング集合の要素 $x^{(i)}$ を用いる． すなわち，トレーニングデータ $x^{(i)},y^{(i)} (i=1,2,\ldots,m)$ に対して $\ell^{(i)} =x^{(i)}$ と設定するのである． このとき，トレーニングデータ $x^{(i)}$ に対する特性ベクトル $f^{(i)} =(f_0^{(i)}, f_1^{(i)}, \ldots, f_m^{(i)})^T$ は， 各ランドマークへの近さを表すベクトルとして， $$ \begin{array}{c} f_0^{(i)} =1 \\ f_1^{(i)} = \exp \left( - \frac{ \| x^{(i)}-\ell^{(1)} \|}{2 \sigma^2} \right) \\ \vdots \\ f_i^{(i)} = \exp \left( - \frac{ \| x^{(i)}-\ell^{(i)} \|}{2 \sigma^2} \right) = \exp (0) = 1 \\ \vdots \\ f_m^{(i)} = \exp \left( - \frac{ \| x^{(i)}-\ell^{(m)} \|}{2 \sigma^2} \right) \end{array} $$ と計算される $m+1$次元ベクトルである．ここで $f_0^{(i)}$ はランドマークに依存しない定数項を表す （線形回帰における $y$-切片を表す）パラメータであり，常に $1$ に設定しておく． これは元のトレーニングデータである $n$次元ベクトル $x^{(i)}$ の空間から $m$次元の特性ベクトル $f^{(i)}$ の空間に射影したものと考えられる．
費用関数 $J(w)$ は， ランドマークに対する重みベクトル $w=(w_0, w_1, \ldots,w_m)$ の関数として， $$ J(w)= C \displaystyle\sum_{i=1}^m \left[ y^{(i)} cost_1(w^T f^{(i)}) + (1-y^{(i)}) cost_0(w^T f^{(i)}) \right] +\frac{1}{2} \displaystyle\sum_{j=0}^m w_j^2 $$ となる． 2項目の$\sum_{j=0}^m w_j^2$ は正則化項を表し， ベクトル表記すると $w^T w$ となる． これはカーネルによって定まる行列 $M$ を用いて $w^T M w$ とする場合もある． 仮説関数を用いて実際の分類を行うためには， $w^T f \geq 0$ のときには $y=1$ と予測し，それ以外のときには $y=0$ と予測すれば良い．
前節と同様に， パラメータ $C$ を大きく設定すると誤差は小さくなるが過剰適合の危険性が高まり， $C$ を小さくするとその逆になる． パラメータ $\sigma$ を小さくするとカーネルはランドマークの周辺に集中するので， 誤差は小さくなるが，一方では過剰適合の危険性が高まり，$\sigma$ を大きくするとその逆になる．
元の特性ベクトルの数 $n$ がトレーニングデータの数 $m$ に対して小さいときには， 決定境界の非線形性が強いことが予想されるので， Gaussカーネルを用いたSVMが有効になる．ただし計算量はそこそこかかるので， $m$ が非常に大きいときには，特性ベクトルの数 $n$ を増やしてから比較的高速な ロジスティクス回帰や線形カーネルを用いたSVMを用いた方が良い．
カーネルはクラス SVC (Support Vector Classifier）の引数 kernel で与える．
Gaussカーネル（引数はrbf）を用いると非線形な境界をもつ問題に対しても、分類が可能になる．
irisデータで試してみよう．
import plotly.express as px
iris = px.data.iris()
# 独立変数（特徴ベクトル） X
X = iris[["sepal_length", "sepal_width", "petal_length", "petal_width"]]
# 従属変数 y
y = iris["species_id"]
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3)
from sklearn.svm import SVC # クラスのインポート
svc = SVC(kernel="rbf", gamma="scale") # インスタンス生成
svc.fit(X_train, y_train) # 訓練
yhat = svc.predict(X_test) # 予測
from sklearn import metrics
print(metrics.accuracy_score(y_test, yhat))
from yellowbrick.classifier import ClassificationReport
visualizer = ClassificationReport(svc)
visualizer.fit(X_train, y_train)
visualizer.score(X_test, y_test)
visualizer.show()
from yellowbrick.classifier import ConfusionMatrix
cm = ConfusionMatrix(svc)
cm.fit(X_train, y_train)
cm.score(X_test, y_test)
cm.show()
from yellowbrick.classifier import ROCAUC
visualizer = ROCAUC(svc, size=(600, 400))
visualizer.fit(X_train, y_train)
visualizer.score(X_test, y_test)
visualizer.show()
from sklearn.preprocessing import StandardScaler
from sklearn.neighbors import KNeighborsClassifier
from yellowbrick.contrib.classifier import DecisionViz
X = iris[["sepal_width", "petal_width"]] # 2次元を切り出す
X = StandardScaler().fit_transform(X) # 可視化のためにスケーリングしておく
viz = DecisionViz(
svc,
features=["sepal_width", "petal_width"],
classes=["Iris-setosa", "Iris-versicolor", "Iris-virginica"],
)
viz.fit(X, y)
viz.draw(X, y)
viz.show()
仮説の評価（交差検証）
トレーニングデータだけを用いて訓練すると，他のデータに対して役に立たない危険性がある． このような過剰適合を避けるためには，特徴の数を減らしたり，正則化パラメータを調整したりする方法があるが， ここではどのモデルを採用すれば良いかを評価するための方法論について述べる．
学習した仮説関数が正しいかどうかを検証するための基本は，データを分けることである． 最も簡単なのは，トレーニングデータとテストデータの2つに分け，訓練（費用関数を最小化する重み $w$ の決定）はトレーニングデータで行い， 評価（正解データとの誤差を計算）はテストデータで行う方法である．
データの数が足りない場合には，トレーニングデータとテストデータを入れ替えて行う方法が行われる． これを 交差検証 (cross validation)とよぶ．
上で用いたサポートベクトル分類 svcのインスタンスで試してみる． 引数の cv は交差検証の分割数であり， 引数の scoring="accuracy" は評価尺度を正解率に設定している （回帰分析のときには，引数を決定変数 "r2" や最大誤差 "max_error" などに設定する）．
以下ではデータを $10$ に分割し，そのうちの１つをテストデータ（残りをトレーニングデータ）として $10$ 個の正解率を計算し，その平均値を出している．
from sklearn.model_selection import cross_val_score
scores = cross_val_score(svc, X, y, cv=10, scoring="accuracy")
print(scores)
print(scores.mean())
ニューラルネット
ここでは，古典的なニューラルネットを機械学習の観点から解説する．
ニューラルネットワークは脳の動きを模した極めて単純な構造から構成される． 基本となるのは， ロジスティック回帰の仮説関数 $h_w(x)$ である． これは，入力ベクトル $x=(x_1,x_2,\ldots,x_n)$ を出力 $h_w(x)$ に変換する． $$ h_w(x)= g(w_0 +w_1 x_1 + w_2 x_2 + \cdots + w_n x_n ) $$ これは，入力ベクトルに行列を乗じた後に，非線形変換 $g(z)$ をしているだけである． ここで $g(z)$ は，活性化関数 (activation function)と呼ばれる関数であり，通常は非線形関数を用いる．
代表的な活性化関数を以下に示す．
- ロジスティック（シグモイド）関数 $\sigma$： 値域は $(0,1)$ $$ \frac{1}{1+e^{-z}} $$
- ReLU(rectified linear unit) $(x)^+$: 値域は $[0,\infty]$ $$ \max(0,z) $$
- $\tanh$ (hyperbolic tangent): 値域は $(-1,1)$ $$ \frac{e^z- e^{-z}}{e^z + e^{-z}} $$
ニューラルネットワークでは，上の機構で入力を出力に変換する点を，複数の層に並べることによって学習を行う． 線形関数の重み $w$ は枝上に定義されており， 左端の1層目に入力 $x$ を与えることで計算された出力が2層目の入力となり， 2層目の出力が3層目の入力になる． 最後の（右端の）層の出力が $y$ となる．このように計算を行う関数が， ニューラルネットワークの仮説関数となる．
分類を例としてニューラルネットワークの学習（重みの最適化）を考える． トレーニングデータの集合を $x^{(i)}, y^{(i)} (i=1,2,\ldots,m)$ とする． 各層 $\ell~(=1,2,\ldots,L)$ における点の数を $s_{\ell}$ と記す．
$y=0,1$ に分類する場合には最後の層 $L$ は1つの点を含む（すなわち $s_{L}=1$ になる）． $K$個のクラスに分類する場合には，最後の層 $L$ は $K$個の点を含む（すなわち $s_{L}=K$ になる）．
2つのクラスへの分類（$y=0,1$）のときの費用関数は，ロジスティック回帰と同様に， $$ J(w)= \frac{1}{m} [ \sum_{i=1}^m y^{(i)} \log h_{w}(x^{(i)}) + (1-y^{(i)})\log (1- h_{w}(x^{(i)}) )]+ \frac{\lambda}{m} \sum_{\ell=1}^{L-1} \sum_{i=1}^{s_{\ell}} \sum_{j=1}^{s_{\ell+1}} (w_{ji}^{(\ell)})^2 $$ と書くことができる． ここで $h_{w}(x^{(i)})$ はトレーニングデータ $x^{(i)}$ に対する最後の層の出力であり， $w_{ji}^{(\ell)}$ は第 $\ell$ 層の点 $i$ から 第 $\ell+1$ 層の点 $j$ への枝上の重みである． 上の費用関数では，正則化のために，重みの自乗平均に 正則化パラメータ $\lambda$ を乗じたものを加えている．
費用関数 $J(w)$ を最小化するためには，費用関数の勾配ベクトルを計算する必要がある． ニューラルネットワークでは，逆伝播(backpropagation；バックプロパゲーション)とよばれる方法を用いて勾配を計算する． 勾配を用いて費用関数の最小化を行うことができるが， ニューラルネットワークの場合には費用関数 $J(w)$ は非凸であり， 局所的最適解に停留する可能性があることに注意する必要がある．
最近では，ニューラルネットは深層学習と進化し， ロジスティック関数以外の活性化関数が使われるようになり，実用化が進んでいる． 深層学習は， http://playground.tensorflow.org/ で体験することができる．
from sklearn.neural_network import MLPClassifier # ニューラルネット
neural = MLPClassifier()
neural.fit(X_train, y_train) # 訓練
yhat = neural.predict(X_test) # 予測
print(metrics.accuracy_score(y_test, yhat))
scores = cross_val_score(neural, X, y, cv=10, scoring="accuracy")
print(scores)
print(scores.mean())
単純Bayes
単純Bayes(naive Bayes)は，名前の通り，高校で習うBayesの定理（規則）(Bayes' theorem; Bayes' rule)を用いた手法である． $P(X)$ を事象 $X$ が起きる確率，$P(Y|X)$ を事象 $X$ が起きたもとで事象 $Y$ が起きる確率 としたとき，Bayesの定理は以下の式で表される． $$ P(X) P(Y|X) =P(Y)P(X|Y) $$
分類の目的は，特性ベクトル $x$ に対して出力が $y$ になる確率 $P(y|x)$ を推定することであった． これはBayesの定理を用いると， $$ P(y|x) = \frac{ P(y) P(x|y) }{P(x)} $$ と計算できるので，我々の目的は $P(x|y)$ を推定することに置き換えられる． つまり，$y=0,1$ の各々に対して， トレーニングデータのしたがう分布 $p(x|y)$ を推定すれば良い．
いま，データは $n$次元実数ベクトルであり，各成分 $j=1,2,\ldots,n$ は独立な 正規(Gauss)分布にしたがっているものとする「単純な」仮定をする （これが単純Bayesの名前の由来である）．
平均 $\mu$，標準偏差 $\sigma$ の正規分布 $N(\mu, \sigma^2)$ の密度関数を $$ p(x; \mu, \sigma^2)= \frac{1}{\sqrt{2\pi\sigma^{2}}} \exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2} \right) $$ と記す．
$y$ が $1$ であるトレーニング集合を $\{ x^{(1)},x^{(2)}, \ldots, x^{(m)} \}$ とし， 各$j=1,2,\ldots,n$ に対して平均と分散を推定する． $$ \mu_j=\frac{1}{m} \displaystyle\sum_{i=1}^{m} x_j^{(i)} $$ $$ \sigma_j^2 = \frac{1}{m} \displaystyle\sum_{i=1}^{m} (x_j^{(i)}-\mu_j)^2 $$ 1式目は平均，2式目は分散を表す． 統計学者は不偏分散（$1/m$ のかわりに $1/(m-1)$）を用いる場合が多いが， 機械学習では上の定義を用いる場合が多い．
成分ごとに独立な正規分布という仮定の下では，$y=1$ に分類されるデータの密度関数は， 正規分布の密度関数の積として $$ P(x|y=1) = \displaystyle\prod_{j=1}^n p(x_j; \mu_j, \sigma_j^2) $$ と書くことができる． 同様に $y=0$ に対しても $P(x|y=0)$ が計算できる．
これらの情報を用いることによって， （Bayesの定理を用いると）与えられたデータ $x$ が $y=1$ に分類される確率は， $$ P(y=1|x) = \frac{ P(y=1) P(x|y=1) }{P(x)} $$ と計算できる．ここで分母の $P(x)$ も $$ P(x)= P(x|y=1) P(y=1) + P(x|y=0) P(y=0) $$ と計算できることに注意されたい．
単純Bayes法はデータがしたがう確率分布についての情報をもつ場合に有効であり， 分類だけでなく不適合値の検出などの応用にも用いることができる．
from sklearn.naive_bayes import GaussianNB # 単純Bayes
bayes = GaussianNB()
bayes.fit(X_train, y_train) # 訓練
yhat = bayes.predict(X_test) # 予測
print(metrics.accuracy_score(y_test, yhat))
scores = cross_val_score(bayes, X, y, cv=10, scoring="accuracy")
print(scores)
print(scores.mean())
from sklearn import tree
tree_class = tree.DecisionTreeClassifier() # 決定木
tree_class.fit(X_train, y_train) # 訓練
yhat = tree_class.predict(X_test) # 予測
print(metrics.accuracy_score(y_test, yhat))
scores = cross_val_score(tree_class, X, y, cv=10, scoring="accuracy")
print(scores)
print(scores.mean())
アンサンブル法
アンサンブル法(ensemble method)とは，複数のアルゴリズムを合わせることによって， より高精度の予測もしくは分類を得るためのメタアルゴリズムである． その基本原理は「ばらつきをもったデータを集約するとばらつきが減少する」という 統計の公理であり，様々なデータや手法で得られた結果を， 平均したり多数決したりすることによって，よりばらつきの少ない結果を得ようというものである．
バギング(bagging; bootstrap aggregatingの略）は，ブートストラップ(bootstrap)とよばれるリサンプリング法によって複数のトレーニング集合を作成し， それらのデータを用いて作成したモデルを平均（回帰の場合）もしくは多数決（分類の場合）によって統合する方法である．
ブースティング(boosting)は，複数の手法の重み付けを学習によって調整する方法である．
ランダム森(random forest)は，決定木ベースのバギングである．
from sklearn.ensemble import RandomForestClassifier # ランダム森
forest = RandomForestClassifier()
forest.fit(X_train, y_train) # 訓練
yhat = forest.predict(X_test) # 予測
print(metrics.accuracy_score(y_test, yhat))
scores = cross_val_score(forest, X, y, cv=10, scoring="accuracy")
print(scores)
print(scores.mean())
クラスタリング
UCI機械学習レポジトリのワインに関するデータセットを用いてクラスタリングを解説する．
使用するのはKMeansクラスで実装されている $k$-平均法である．
$k$-平均法($k$-means method)は，与えられたデータを $k$個の部分集合 （クラスター）に分けるための手法である． $k$-平均法は， クラスター数 $k$ と$n$次元実数ベクトル $x^{(i)} \in \mathbf{R}^n$ を要素とする トレーニング集合 $\{ x^{(1)}, x^{(2)}, \ldots, x^{(m)} \}$ を 与えたとき，$k$個のクラスター（データの部分集合）の 重心(centroid) $\mu_1,\mu_2, \ldots ,\mu_k$ を（近似的に）求めることを目的とする． 重心はクラスター内のデータの平均値をとることによって得られるので，この名前がついた．
アルゴリズムは単純である．まずランダムに重心 $\mu_1, \mu_2, \ldots, \mu_k \in \mathbf{R}^n$ を決める． 次に，トレーニング集合内のデータを最も近い重心に割り当てることによってクラスターを決める． そして，各クラスターの重心を新しい重心として，重心が移動しなくなるまで，この操作を繰り返す．
$k$-平均法で対象とする問題は$NP$-困難であるので，このアルゴリズムは局所的最適解を 求めるだけである．初期重心をランダムに変更して繰り返すことによって， より良い解を求めることができるが，大域的最適解を得られる保証はない．
データは http://logopt.com/data/wine.data にある．
列名は https://archive.ics.uci.edu/ml/datasets/Wine で解説されている．
以下では，$k$-平均法で4つのクラスターに分類する．
L = [
"Alcohol",
"Malic",
"Ash",
"Alcalinity",
"Magnesium",
"Phenols",
"Flavanoids",
"Nonflavanoid",
"Proanthocyanins",
"Color",
"Hue",
"OD280",
"OD315",
"Proline"
]
wine = pd.read_csv("http://logopt.com/data/wine.data", names=L)
wine.head()
from sklearn.cluster import KMeans # クラスをインポート
kmeans = KMeans(n_clusters=4) # インスタンス生成
kmeans.fit(wine) # 訓練
wine["label"] = kmeans.labels_
wine.head()
パラメータ $k$ の適正化 （エルボー法）
上では4つのクラスターに分けたが， 幾つに分けるかはデータに依存する． 目安になる方法としてエルボー法がある．
yellowbrickパッケージのKElbowVisualizerクラスを用いて可視化する．
クラスの引数のkで調べる $k$ の範囲を指定する．以下では，$2$から$12$ の範囲で $k$-平均法を適用し，結果をプロットしている．
評価尺度（重心への距離の和）が $k$ の増加とともに減少するが， それが減少が緩やかになった $k$ が適正なクラスター数であると考える．
この図が肘（エルボー）を曲げたように見えることからエルボー法とよばれ．
from yellowbrick.cluster import KElbowVisualizer
visualizer = KElbowVisualizer(kmeans, k=(2, 12))
visualizer.fit(wine)
visualizer.show()
from yellowbrick.cluster import InterclusterDistance
kmeans = KMeans(n_clusters=4, random_state=123)
visualizer = InterclusterDistance(kmeans)
visualizer.fit(wine)
visualizer.show()
次元削減
ここでは、高次元のデータを低次元に変換するための方法について学ぶ。
多くの機械学習の実際問題は多次元データである。これを人間が見てわかる程度（通常は2次元か3次元）に（本質を失うことなしに）変換できれば、 問題に対する洞察を得ることができ、大変便利である。また、低次元に変換してから、回帰や分類を行う方が、良い結果を出す場合もある。
主成分分析
最初に紹介するのは、古典的な主成分分析 (principal component analysis； PCA)である。 これは、低次元に射影変換したときに、なるべくデータがばらつくような軸（主成分）を選択するものである。
主成分分析は， $n$次元実数ベクトル $x^{(i)} \in \mathbf{R}^n$ を要素とする トレーニング集合 $\{ x^{(1)}, x^{(2)}, \ldots, x^{(m)} \}$ を， $k(<n)$次元平面上に射影したときの距離の自乗和を最小にするような $k$次元平面を求める ことによって次元削減を行う．
行列を用いると，トレーニング集合を表す $m \times n$ 行列 $X$ は， $$ X = \left( \begin{array}{cccc} x^{(1)}_{1} & \ldots & x^{(1)}_{n} \\ x^{(2)}_{1} & \ldots & x^{(2)}_{n} \\ \vdots & \ddots & \vdots \\ x^{(m)}_{1} & \ldots & x^{(m)}_{n} \end{array} \right) $$ となる．これを $m \times k$ 行列 $Z$ に，$n \times k$ の行列 $P$ で $$ Z = X P $$ と射影することによって次元削減が行われる．
アルゴリズムを適用する前にデータの正規化を行っておく． つまり，各特性の平均値 $\mu_j (j=1,2,\ldots,n)$ を $$ \mu_j =\frac{1}{m} \displaystyle\sum_{i=1}^m x_j^{(i)} $$ と計算し，新たな $x_j^{(i)}$ を $x_j^{(i)}-\mu_j$ もしくは平均値で正規化して $$ \frac{x_j^{(i)}-\mu_j}{\mu_j} $$ と設定する． これはデータの平均を原点に移動させたことに相当し， 特性スケーリング(feature scaling)もしくは平均正規化(mean normalization)とよばれる．
次に，データの共分散行列 $$ \Sigma =\frac{1}{m} \displaystyle\sum_{i=1}^m (x^{(i)}) (x^{(i)})^T $$ の固有値分解 $\Sigma = U S U^T$ を求める．ここで $U$ は直交行列であり， $S$ は固有値を対角成分にもつ対角行列である． ここで，正方行列 $\Sigma$ の固有値と固有ベクトルとは， $$ \Sigma u = \lambda u $$ を満たすスカラー $\lambda$ と（長さが $1$ で $u^T u=0$ を満たす）正規直交ベクトル $u$ である．
以下では，$k=1$ の場合のみを解説する．このとき，$Z$ と $P$ はベクトルになるので， $$ z = X p $$ が射影変換になる．$z$ の分散が最大になる射影 $p$ が，射影したときの距離の和を最小にする． $$ z^T z = (Xp)^T (Xp) =p^T X^T X p = m p^T \Sigma p $$ であり，$\Sigma$ の最大固有値を $\lambda$，そのときの固有ベクトルを $u$ とすると， $$ u^T \Sigma u = u^T \lambda u =\lambda $$ であるので，分散を最大にするには，$p=u$ にとれば良いことが分かる． $k>1$ の場合も同様に計算ができ， $n \times n$ 行列 $U$ の（固有値が大きい順に並べた）最初の $k$個の列が， 求めたい平面への射影を表すベクトルになる．
irisデータセットを用いて主成分分析と可視化の方法を説明する．
iris = px.data.iris()
X = iris[["sepal_length", "sepal_width", "petal_length", "petal_width"]]
from sklearn.decomposition import PCA
pca = PCA(n_components=2) # 2次元に射影
pca.fit(X)
pca.components_ # 射影行列
Z = pca.transform(X) # 射影したデータ
# 射影した2次元データを元のデータフレームに追加
iris["X"] = Z[:, 0]
iris["Y"] = Z[:, 1]
iris.head()
import seaborn as sns
sns.lmplot(x="X", y="Y", hue="species", fit_reg=False, data=iris);
from sklearn.manifold import TSNE
Z = TSNE(n_components=2, random_state=0).fit_transform(X)
# 射影した2次元データを元のデータフレームに追加
iris["X"] = Z[:, 0]
iris["Y"] = Z[:, 1]
sns.lmplot(x="X", y="Y", hue="species", fit_reg=False, data=iris)
問題（アルコール摂取量）
http://logopt.com/data/drinks.csv にある国別のアルコール摂取量データを用いて主成分分析を行え．
4次元の数値データを2次元に射影し，2次元座標で表示せよ．その際，色調としては，大陸(continent)列を用いよ．
from sklearn.ensemble import RandomForestRegressor
data = pd.read_csv(
"http://logopt.com/data/Advertising.csv", index_col=0
) # 0行目をインデックスにする．
X = data[["TV", "Radio", "Newspaper"]]
y = data["Sales"]
rf = RandomForestRegressor()
rf.fit(X, y);
from sklearn.inspection import permutation_importance
result = permutation_importance(rf, X, y)
result.importances_mean
結果から， TVが最も重要で，次いで Radio であり， Newspaperは広告効果がほとんどないことが分かる．
yellowbrickパッケージのFeatureImportancesクラスを用いて可視化する．
from yellowbrick.model_selection import FeatureImportances
viz = FeatureImportances(rf)
viz.fit(X, y)
viz.show()
import plotly.express as px
iris = px.data.iris()
X = iris[["sepal_length", "sepal_width", "petal_length", "petal_width"]]
y = iris["species_id"]
問題（胸部癌）
http://logopt.com/data/cancer.csv にある胸部癌か否かを判定するデータセットに対してランダム森で分類し， 特徴重要度を計算し，可視化せよ．
最初の列diagnosisが癌か否かを表すものであり，Mが悪性（malignant），Bが良性（benign）を表す．
決定木の可視化
決定木 (decision tree) を可視化することができる．
以下は、linuxなら簡単だが，MacやWindowsだとインストールが面倒． Google Colab.の場合には，以下のセルのコメントを外して実行するだけで良い．
MacやWindowsの場合には，以下のサイトを参照して， graphvizとdtreevizをインストールする。
https://github.com/parrt/dtreeviz
dtreeviz関数に，引数として決定木回帰 DecisionTreeRegressor のインスタンス，X, y, 特徴名のリスト feature_names，ターゲット名 target_name を与えると，回帰の決定木を可視化してくれる．
また，引数Xに特徴ベクトルをいれると，それが決定木のどのように辿って予測を行うかを示してくれる．
#!sudo apt install graphviz
#!pip install dtreeviz
from sklearn import tree
data = pd.read_csv(
"http://logopt.com/data/Advertising.csv", index_col=0
) # 0行目をインデックスにする．
X = data[["TV", "Radio", "Newspaper"]]
y = data["Sales"]
reg = tree.DecisionTreeRegressor(max_depth=2) # limit depth of tree
reg.fit(X, y)
from dtreeviz.trees import dtreeviz
viz = dtreeviz(reg, X, y, feature_names=data.columns, target_name="Sales")
display(viz)
viz = dtreeviz(reg, X, y, feature_names=data.columns, target_name="Sales", X=X.iloc[0])
print(y.iloc[0])
display(viz)
irisデータセットに対する決定木による分類も同様に可視化できる．
dtreeviz関数に，引数として決定分類 DecisionTreeClassifier のインスタンス，X, y, 特徴名のリスト feature_names，ターゲット名 target_nameならびにクラス名class_namesを与えると，分類の決定木を可視化してくれる．
また，引数 X に特徴ベクトルをいれると，それが決定木のどのように辿って予測を行うかを示してくれる．
from sklearn import tree
from dtreeviz.trees import dtreeviz
import plotly.express as px
iris = px.data.iris()
X = iris[["sepal_length", "sepal_width", "petal_length", "petal_width"]]
y = iris["species_id"]
classifier = tree.DecisionTreeClassifier(max_depth=2) # 木の深さを2に設定
classifier.fit(X, y)
viz = dtreeviz(
classifier,
X,
y,
target_name="species_id",
feature_names=["sepal_length", "sepal_width", "petal_length", "petal_width"],
class_names=["setosa", "versicolor", "virginica"],
)
display(viz)
viz = dtreeviz(
classifier,
X,
y,
target_name="species_id",
feature_names=["sepal_length", "sepal_width", "petal_length", "petal_width"],
class_names=["setosa", "versicolor", "virginica"],
X=X.iloc[100],
)
display(viz)
問題（胸部癌）
"http://logopt.com/data/cancer.csv" にある胸部癌か否かを判定するデータセットに対する決定木（深さは4）を可視化せよ．
最初の列diagnosisが癌か否かを表すものであり，"M"が悪性（malignant），"B"が良性（benign）を表す．
なお，dtreevizではターゲットは数値である必要があるので，LabelEncoderで数値化しておく．
cancer = pd.read_csv("http://logopt.com/data/cancer.csv", index_col=0)
y = cancer["diagnosis"]
X = cancer.drop("diagnosis", axis=1)
from sklearn.preprocessing import LabelEncoder
y = LabelEncoder().fit_transform(y)