1 人工知能と機械学習

1.1 人工知能とは

人工知能（artificial intelligence）とは、人間の知的なふるまいをプログラムで再現したもののことである。

例えば「画像に何が写っているのかを識別する」というふるまいをコンピュータで実行するプログラムは人工知能といえる。

人工知能を作るための方法はいくつかあり、機械学習が登場する前のルールベースと呼ばれるアプローチでは、画像に何が写っているのかを推論する部分のプログラム（推論のルール）を人間が考えて作るものであった。

ルールベースの人工知能は、医者など専門家の知的なふるまいをシステム化することを目的としたエキスパートシステムという試みとして1970~80年代に発展したものの、推論のルール（人工知能の中身）を人間が作るため、システムの開発・保守運用に多数のプログラマーを必要として莫大な人件費がかかるものであった。

1.2 機械学習とは

他方で機械学習（machine learning）は人工知能のプログラムを人間が考えるのではなくデータから統計学的に推定させる（機械に学習させる）タイプのアプローチである。

このアプローチは大量のデータと高い処理能力のコンピュータが必要になるが、人間の労働力は少なく済む。現在「人工知能」と呼ばれるものはほとんどが機械学習に基づくもので、これまでに触れてきた「線形回帰（回帰分析）」や近年ブームになっている「ディープラーニング」も機械学習の手法のひとつである。

1.3 学習の種類

学習の種類は大別して次の３種類がある。

教師あり学習（supervised learning）：予測（推論結果）の正解データ（教師データ）があるもの
- 回帰分析のように、入力データ（説明変数）\(x\)と出力データ（被説明変数）\(y\)の相関関係を学習する
教師なし学習（unsupervised learning）：予測の正解データが与えられないもの
- 入力データの規則性から「同じ系統のもの」「異なる系統のもの」を分けるような推論を行う
強化学習（reinforcement learning）：予測の正解データは無いが、利得関数という間接的な正解データがあるもの
- 囲碁や将棋のAIなどで使われる

本講義では主に教師あり学習に関して述べていく。

1.4 機械学習のタスク

教師あり学習が行う主なタスクは次の２種類の予測問題になる。

回帰（regression）：入力データから連続変数の出力を予測する
- 例：過去の気象データから明日の気温を予測する
分類（classification）：入力データから離散変数の出力を予測する
- 例：顧客情報から、その顧客が商品を「買う」か「買わない」かを予測する

今回は回帰を行う機械学習の例として線形回帰に再度触れ、分類を行う機械学習の例としてロジスティック回帰を紹介する。

2 線形回帰

線形回帰（linear regression）は、予測の目的変数（objective variable, target variable）\(y_i\)と特徴量（feature）\(x_{i1}, x_{i2}, ..., x_{im}\)の間に次のような線形関係を仮定したモデルを建てて予測を行う¹。

\[ y_i = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \cdots + \beta_m x_{im} + \varepsilon_i \]

ここで\(\varepsilon_i\)は誤差項である。

2.1 パラメータの推定方法

2.1.1 二乗誤差の最小化

教師あり学習では予測の悪さの指標として誤差関数（error function）あるいは損失関数（loss function）と呼ばれる関数\(L(\boldsymbol{\theta})\)を定義し、誤差を最小化することを主な基準としてパラメータを推定していく。

線形回帰のパラメータ\(\boldsymbol{\beta} = (\beta_1, \beta_2, \dots, \beta_m)^\top\)の推定では、実測値\(y_i\)と予測値\(\hat{y}_i\)の二乗誤差\((y_i - \hat{y}_i)^2\)の和

\[ L(\boldsymbol{\beta}) = \sum_{i=1}^n (y_i - \hat{y}_i)^2 \]

を誤差関数に用いて、誤差\(L(\boldsymbol{\beta})\)を最小にするパラメータ

\[ % 定義 \newcommand{\argmin}{\mathop{\rm arg~min}\limits} \hat{\boldsymbol{\beta}}^{OLS} = \argmin_{\boldsymbol{\beta}} L(\boldsymbol{\beta}) \]

採用する（なお、\(\argmin\)という記号は「関数の最小値を与える引数」という意味である）。

二乗和誤差\(L(\boldsymbol{\beta})\)は（パラメータが1つの場合）U字型のカーブを描く関数になる。

最小値が１つであり微分可能なので、推定値\(\hat{\boldsymbol{\beta}}^{OLS}\)を得るための計算が容易であるという便利な特性をもっている。

2.1.2 推定量

以下では説明の簡単のために、特徴量が１つの単回帰モデル（simple regression model）

\[ y_i = \beta_0 + \beta_1 x_{i} + \varepsilon_i \]

の場合で説明していく。

最小二乗法では二乗和誤差\(L(\boldsymbol{\beta})\)を最小化するようなパラメータ\(\beta_j^*\)の値を推定するのが目的である。\(L(\boldsymbol{\beta})\)は二次元で見たときU字型の関数であるため、最小値は次の条件（接線の傾きが0であること）を満たすパラメータとなる。

\[ \begin{cases} \frac{\partial L(\boldsymbol{\beta})}{\partial \beta_0} = -2 \sum_{i=1}^n (y_i - \beta_0 - \beta_1 x_i) = 0\\ \frac{\partial L(\boldsymbol{\beta})}{\partial \beta_1} = -2 \sum_{i=1}^n (y_i - \beta_0 - \beta_1 x_i)x_i = 0 \end{cases} \]

上記の式を整理すると、最小二乗推定量を求める式

\[ \begin{align} \hat{\beta}_0 &= \bar{y} - \beta_1 \bar{x}\\ \hat{\beta}_1 &= \frac{\sum_{i=1}^n(x_i-\bar{x})(y_i-\bar{y_i})} {\sum_{i=1}^n (x_i - \bar{x})^2} \end{align} \]

が得られる。ここで\(\bar{y}\)と\(\bar{x}\)はそれぞれ\(y_i\)や\(x_i\)の平均値である。

2.2 Rでの実践

2.2.1 パッケージの読み込み

今回の分析で使うパッケージを読み込む。

pacman::p_load(tidyverse, 
               AER, # HousePricesデータを使用
               stargazer, # lmやplmのときの結果表
               MLmetrics, # 機械学習で用いる関数が実装
               carData) # タイタニック号データを使用

2.2.2 データの用意

{AER}パッケージのHousePricesデータを使用し、住宅の価格を予測する問題に取り組むことにする。

data("HousePrices") #AERパッケージのデータ
head(HousePrices)

##   price lotsize bedrooms bathrooms stories driveway recreation fullbase gasheat
## 1 42000    5850        3         1       2      yes         no      yes      no
## 2 38500    4000        2         1       1      yes         no       no      no
## 3 49500    3060        3         1       1      yes         no       no      no
## 4 60500    6650        3         1       2      yes        yes       no      no
## 5 61000    6360        2         1       1      yes         no       no      no
## 6 66000    4160        3         1       1      yes        yes      yes      no
##   aircon garage prefer
## 1     no      1     no
## 2     no      0     no
## 3     no      0     no
## 4     no      0     no
## 5     no      0     no
## 6    yes      0     no

このデータセットには次の変数が含まれている

price：住宅の販売価格
lotsize：敷地面積（平方フィート）
bedrooms：寝室の数
bathrooms：浴室（full bathroom）の数
stories：地下室を除いた階数
driveway：家に私道があるかどうか
recreation：レクリエーションルームがあるか
fullbase：地下室があるか
gasheat：給湯にガスを使用しているか
aircon：セントラルエアコンがあるか
garage：ガレージの数
prefer：都市の近くに立地しているか

機械学習の目的は未知のデータに対する予測性能を高めることである。したがって、データをモデルの学習（パラメータの推定）に使う「訓練用データ」と、学習に使わず予測性能の試験に使う「テスト用データ」に分割しておく。

# ID列を追加
df = HousePrices %>% rownames_to_column("ID")

# 80%を訓練用データに
set.seed(0)
train <- df %>% sample_frac(size = 0.8)

# 訓練用データに使っていないIDの行をテスト用データに
test <- anti_join(df, train, by = "ID")

ID列は予測に使わないため削除しておく

# ID列の削除
train <- train %>% select(-ID)
test <- test %>% select(-ID)

2.2.3 線形回帰の学習

線形回帰はlm()関数で行うことができる。

reg <- lm(log(price) ~ . , data = train)
stargazer(reg, type = "text")

## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                             log(price)         
## -----------------------------------------------
## lotsize                      0.0001***         
##                              (0.00001)         
##                                                
## bedrooms                       0.024           
##                               (0.016)          
##                                                
## bathrooms                    0.177***          
##                               (0.023)          
##                                                
## stories                      0.097***          
##                               (0.014)          
##                                                
## drivewayyes                  0.116***          
##                               (0.032)          
##                                                
## recreationyes                0.090***          
##                               (0.031)          
##                                                
## fullbaseyes                  0.068***          
##                               (0.026)          
##                                                
## gasheatyes                   0.215***          
##                               (0.056)          
##                                                
## airconyes                    0.178***          
##                               (0.024)          
##                                                
## garage                       0.056***          
##                               (0.013)          
##                                                
## preferyes                    0.143***          
##                               (0.026)          
##                                                
## Constant                     10.041***         
##                               (0.053)          
##                                                
## -----------------------------------------------
## Observations                    437            
## R2                             0.683           
## Adjusted R2                    0.675           
## Residual Std. Error      0.215 (df = 425)      
## F Statistic          83.139*** (df = 11; 425)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

2.2.4 予測精度の評価

lmオブジェクトから予測値を計算させるにはpredict()関数にlmオブジェクトとデータフレームを入れて実行すればよい。

# 既知のデータ（train）に対する予測値
y_train_pred = predict(reg, train) %>% exp()  # 学習時に対数変換していたため予測値を指数変換して戻す
# 実測値
y_train_true = train$price

予測精度の評価は、例えば平均二乗誤差（Mean of Squared Error: MSE）のような関数で評価を行う。

\[ MSE = \frac{1}{n} \sum_{i=1}^n (y_i - \hat{y}_i)^2 \]

なお、ここでの「誤差」とは、\(y_i - \hat{y}_i\)すなわち観測値と予測値の差であって、統計学や計量経済学でいうところの誤差項(\(\varepsilon\))とは異なる点に注意されたい。

{MLmetrics}パッケージにはMSEのような機械学習で用いる関数が実装されているため、このパッケージをインストールして読み込む

# 平均二乗誤差（MSE）による予測精度の評価
MSE(y_train_pred, y_train_true)

## [1] 229164893

MSEは二乗しているため、非常に値が大きくなってしまう。

MSEの平方根をとることでこの点を修正したものが\(RMSE\)(root mean squared error)である。

\[ RMSE = \sqrt{\frac{1}{n} \sum_{i=1}^n (y_i - \hat{y}_i)^2} \]

# 平方根平均二乗誤差（RMSE）による予測精度の評価
RMSE(y_train_pred, y_train_true)

## [1] 15138.19

訓練用データのpriceの予測に際して約1.5万の誤差があることを示している。

テスト用データに対する予測誤差

回帰モデルにとって未知のデータであるtestデータセットに対する予測値をみてみよう。

# 未知のデータ（test）に対する予測値
y_pred <- predict(reg, test) %>% exp()
# 実測値
y_test = test$price

# 平方根平均二乗誤差（RMSE）による予測精度の評価
RMSE(y_pred, y_test)

## [1] 15036

データを訓練用とテスト用に任意の割合で分割してテスト用データに対する予測誤差を測る方法をホールドアウト（hold out）法といい、その誤差の推定量をホールドアウト誤差（hold-out error）という。

一方で最初に行ったような訓練用データに対する予測誤差を訓練誤差（training error）という。訓練誤差は一般に真の誤差（真の分布で誤差関数の期待値をとった誤差）以下になる推定量であることが知られているため、機械学習の分野ではあまり重視されない。

一方でホールドアウト誤差はデータを分割しなければならないためサンプル数が小さい状況と相性が悪い。この点を改善する方法として交差検証（cross validation）という方法がある。詳細は次回に述べる。

（参考）その他の評価指標

平均絶対誤差（mean absolute error: MAE）

\[ MAE = \frac{1}{n} \sum_{i=1}^n |y_i - \hat{y}_i| \]

MAE(y_test, y_pred)

## [1] 10568.73

絶対誤差は二乗誤差に比べると外れ値の影響を受けにくい特徴を持つため、予測を大きく外した個体の評価への影響を抑えたい状況ではMAEのほうが使いやすい場合がある。

また、二乗して平方根をとっているRMSEに比べると直感的で多くの人にとって比較的わかりやすい指標でもある。

決定係数（coefficient of determination）

回帰分析の当てはまりの指標として習う決定係数\(R^2\)も回帰モデルの予測性能の指標として使うことができる。

\[ R^2 = 1 - \frac{ \sum^n_{i=1} (y_i - \hat{y})^2 }{ \sum^n_{i=1} (y_i - \bar{y})^2 } = 1 - \frac{ 残差の変動 }{ データの変動 } \]

R2_Score(y_test, y_pred)

## [1] 0.5130499

3 ロジスティック回帰

ロジスティック回帰（logistic regression）は、目的変数に\(\{0,1\}\)のような二値変数や確率を使う回帰モデルである。

特徴量\(x_{i1}, x_{i2}, ..., x_{im}\)の下で目的変数\(y_i\)が1の値をとる確率\(P(y_i = 1| x_{i1}, x_{i2}, ..., x_{im})\)を\(p\)、\(P(y_i = 0| x_{i1}, x_{i2}, ..., x_{im})\)を\(1-p\)とおくと、ロジスティック回帰モデルは

\[ \log \left( \frac{p}{1-p} \right) = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_m x_{im} \]

あるいは両辺の指数をとって整理した

\[ P(y_i = 1| x_{i1}, x_{i2}, ..., x_{im}) = \frac{\exp(\beta_0 + \beta_1 x_{i1} + \cdots + \beta_m x_{im})} {1 + \exp(\beta_0 + \beta_1 x_{i1} + \cdots + \beta_m x_{im})} \] と表される。

後者の式の右辺はロジスティック・シグモイド関数（logistic sigmoid function）あるいは単にシグモイド関数（sigmoid function）と呼ばれる関数

\[ \sigma(z) = \frac{\exp(z)}{1 + \exp(z)} = \frac{1}{1 + \exp(-z)} \]

で、任意の数を0から1の範囲に変換する関数である。

3.1 パラメータの推定方法

ロジスティック回帰では最尤法（さいゆうほう、method of maximum likelihood）によりパラメータを推定する。

最尤法は「得られたデータはデータが従う確率分布でもっとも出現確率が高い値が生じたもの」という仮定の下で、その確率分布のパラメータとして最も尤もらしい（最も確率が高くなる）パラメータを推定する。

\(n\)個のサンプル\(x_1, x_2, \dots, x_n\)が同一の分布から独立に得られたと仮定する。データが従うと想定する確率分布の確率関数\(P(x|\theta)\)（\(\theta\)は確率分布のパラメータ）の同時確率関数

\[ P(x_1|\theta) \times P(x_2|\theta) \times \cdots \times P(x_n|\theta) = \prod_{i=1}^n P(x_i|\theta) \]

を尤もらしさの指標として用いる。これを尤度関数（likelihood function）といい、\(L(\theta)\)と表す。

3.1.1 例：コイン投げ

例えば、あるコインを4回投げて表が1回出たとする。これを\(p\)の確率で表（\(x=1\)）が、\(1-p\)の確率で裏（\(x=0\)）が出る試行（ベルヌーイ試行）と捉えて、コインの表が出るかどうかは確率関数

\[ P(X = x | p) = \begin{cases} p & x = 1 \text{のとき}\\ 1 - p & x = 0 \text{のとき} \end{cases} \]

のベルヌーイ分布に従うと仮定する（モデリングする）。

4回投げて表が1回、裏が3回という結果（データ）を踏まえると、尤度関数\(L(p)\)は

\[ \begin{align} L(p) &= P(X = 1| p) \times P(X = 0| p) \times P(X = 0| p) \times P(X = 0| p)\\ &= p^1(1-p)^3 \end{align} \]

となる。

この尤度関数を描画してみると、尤度\(L(p)\)が最大になっている点でパラメータが\(p = 0.25 = \frac{1}{4}\)になっていることがわかる。

3.1.2 対数尤度と交差エントロピー誤差

確率の積は非常に小さな値になってしまうため、実際に計算する際は対数で変換した対数尤度（log-likelihood）

\[ \ell(\theta) = \ln L(\theta) = \sum_{i=1}^n \ln P(x_i|\theta) \]

を使用する。

また統計学の分野では尤度という言葉を使うが、機械学習の分野ではベルヌーイ分布の対数尤度の符号を反転させたものを交差エントロピー誤差（cross entropy error）と呼ぶ。

\[ \mathcal{L}(\theta) = - \ell(\theta) = - \ln L(\theta) \]

符号を反転させることで「誤差を最小化するようなパラメータを推定する」という教師あり学習の標準的な枠組みに持っていくことができる。

3.1.3 ニュートン法

最尤法では一般に解析的には解が求まらないため、数値的に解いていく。尤度が高くなる方向（交差エントロピー誤差が下がる方向）へ少しずつ推定値を変える計算を何度も反復することで推定する。

この方法の大まかな流れは次のようになる。まず、初期値\(\hat{\theta}_{0}\)を乱数などで適当に決める。そして\(t\)回目の反復では、現状の推定値\(\hat{\theta}_{t}\)と対数尤度の1次と2次の導関数を用いて次のように推定値を更新する。

\[ \hat{\theta}_{t+1} = \hat{\theta}_{t} - \hat{\theta}_{t}の地点における\mathcal{L}の傾き \]

そして、更新してもあまり変化しなくなったら（ある小さな正の数\(\varepsilon\)について\(|\hat{\theta}_{t+1} - \hat{\theta}_{t}| < \varepsilon\)を満たすなら）反復を終了して\(\hat{\theta}_{t+1}\)を最終的な解とする。

ロジスティック回帰の推定ではニュートン法（Newton’s method）と呼ばれるという方法を用いることが多い。ニュートン法では傾きの計算に\(\mathcal{L}\)の2次の導関数も用いて次のように推定値を更新する。

\[ \hat{\theta}_{t+1} = \hat{\theta}_{t} - \frac { \mathcal{L}'(\hat{\theta}_{t}) } { \mathcal{L}''(\hat{\theta}_{t}) } \]

コイン投げの例で使用した交差エントロピー誤差において実際にこの反復計算を行うと、次の図のようになる（初期値は0.9とした）

3.2 Rでの実践

3.2.1 パッケージの読み込み

今回の分析で使うパッケージを読み込む（すでに一度やっているので、飛ばしても良い）。

pacman::p_load(tidyverse, 
               AER, # HousePricesデータを使用
               stargazer, # lmやplmのときの結果表
               MLmetrics, # 機械学習で用いる関数が実装
               carData) # タイタニック号データを使用

3.2.2 データの用意

{carData}パッケージのTitanicSurvivalデータセットを使う

data("TitanicSurvival")　#carDataに入っている。
head(TitanicSurvival)

##                                 survived    sex     age passengerClass
## Allen, Miss. Elisabeth Walton        yes female 29.0000            1st
## Allison, Master. Hudson Trevor       yes   male  0.9167            1st
## Allison, Miss. Helen Loraine          no female  2.0000            1st
## Allison, Mr. Hudson Joshua Crei       no   male 30.0000            1st
## Allison, Mrs. Hudson J C (Bessi       no female 25.0000            1st
## Anderson, Mr. Harry                  yes   male 48.0000            1st

##                                 survived    sex     age passengerClass
## Allen, Miss. Elisabeth Walton        yes female 29.0000            1st
## Allison, Master. Hudson Trevor       yes   male  0.9167            1st
## Allison, Miss. Helen Loraine          no female  2.0000            1st
## Allison, Mr. Hudson Joshua Crei       no   male 30.0000            1st
## Allison, Mrs. Hudson J C (Bessi       no female 25.0000            1st
## Anderson, Mr. Harry                  yes   male 48.0000            1st

このデータセットは1912年のタイタニック号の沈没事故の乗客の生死に関するデータで、次の変数が含まれている

survived：生存したかどうか
sex：性別
age：年齢（1歳に満たない幼児は小数）。263の欠損値を含む。
passengerClass：船室の等級

このデータセットには欠損値が含まれているため、まず欠損値を除去する

# NA（欠損値）を含む行を削除
titanic <- na.omit(TitanicSurvival)

そしてデータを訓練用・テスト用に分割する。

# ID列を追加
df = titanic %>% rownames_to_column("ID")

# 80%を訓練用データに
set.seed(0)
train <- df %>% sample_frac(size = 0.8)

# 訓練用データに使っていないIDの行をテスト用データに
test <- anti_join(df, train, by = "ID")

# ID列は予測に使わないため削除しておく
train <- train %>% select(-ID)
test <- test %>% select(-ID)

3.2.3 ロジスティック回帰の学習

ロジスティック回帰はglm()関数で実行できる（なお、glm()はパラメータ推定にニュートン・ラフソン法を使用している）。

clf <- glm(survived ~ . , family = binomial(link = "logit"), data = train)
stargazer(clf, type = "text")

## 
## =============================================
##                       Dependent variable:    
##                   ---------------------------
##                            survived          
## ---------------------------------------------
## sexmale                    -2.492***         
##                             (0.186)          
##                                              
## age                        -0.031***         
##                             (0.007)          
##                                              
## passengerClass2nd          -1.187***         
##                             (0.250)          
##                                              
## passengerClass3rd          -2.335***         
##                             (0.253)          
##                                              
## Constant                   3.370***          
##                             (0.359)          
##                                              
## ---------------------------------------------
## Observations                  837            
## Log Likelihood             -387.704          
## Akaike Inf. Crit.           785.408          
## =============================================
## Note:             *p<0.1; **p<0.05; ***p<0.01

推定された係数を見ると、性別が男性であったり、年齢が高かったり、客室の等級が低い場合は生存する確率が下がる関係が見られ、女性や子供を守るように努めたのであろうことが推測できる。

3.2.4 予測精度の評価

glm()関数の予測値もpredict()関数で出力できる。ただし、ロジスティック回帰はデフォルトではオッズ比\(\log \left( \frac{p}{1-p} \right)\)が出される。type = "response"を指定してやることで確率値\(P(y_i = 1| x_{i1}, x_{i2}, ..., x_{ip})\)のほうを出力してくれる。

# 未知のデータ（test）に対する予測値（確率）
p_pred <- predict(clf, test, type = "response")

predict(type = "response")で出力される値は確率値であり、そのままでは使いにくいのでifelse()関数で二値変数に変換する

# 予測した確率が「0.5未満ならno, それ以外はyes」となるよう変換
y_pred <- ifelse(p_pred < 0.5, "no", "yes")

# 実測値
y_test = test$survived

混同行列

分類（被説明変数が離散変数の予測）の評価指標を考える際に基本となるのが混同行列（confusion matrix）である。

これは予測値と実測値のクロス集計表であり、上の図のような構造になっている。

表のPositive, Negativeが予測の被説明変数がとる値の種類で、「実測値がPositiveなものを予測値でもPositiveと予測できていればTrue Positive（真陽性）」という具合に、予測の正誤と予測値の値から４つの区分に分類される。

正しく予測できたものはTrue PositiveやTrue Negativeになり、予測を誤ったものはFalse PositiveかFalse Negativeに入る。

Rではtable()関数を実行すればクロス集計表を作成できる。

# confusion matrix
table(y_test, y_pred)

##       y_pred
## y_test no yes
##    no  97  18
##    yes 30  64

誤った分類をしているものもそれなりにある。

誤分類の程度を測る指標の代表的なものに次の３つの指標がある。

正解率（Accuracy）：正しく予測できた数÷予測した全ての数

\[ Accuracy = \frac{TP + TN}{TP + FP + FN + TN} \]

適合率（precision）：予測したPositiveにおいて正しく予測できたものの割合

\[ Precision = \frac{TP}{TP+FP} \]

再現率（recall）：実際にPositiveであるデータのうち正しく予測できたものの割合

\[ Recall = \frac{TP}{TP+FN} \]

もっともシンプルな考え方で作られていて直感的にもわかりやすいのは正解率である。ただし、世の中の分類問題では、それぞれのカテゴリーの出現頻度が大きく異なることが少なくない（例えばローンを返せなくなる人を予測する場合、そうした人は全体の中のごくわずかである）。そのため適合率・再現率など正解率以外の様々な評価指標が作られている²。

Rでは{MLmetrics}パッケージにそれぞれ用意されている。

# 正解率
Accuracy(y_pred, y_test)

## [1] 0.7703349

# 適合率
Precision(y_test, y_pred)

## [1] 0.7637795

# 再現率
Recall(y_test, y_pred)

## [1] 0.8434783

いずれの指標も7~8割程度には予測できている。

複数の閾値にわたって評価する指標

分類においては「このクラスに属すると思われる確率」を算出し、何らかの閾値を決めて確率値を分類のカテゴリの値に変換して評価を行うことになる。今回の場合、0.5を閾値としていた。

この閾値の取り方によってPrecisionやRecallなどの値も変わってくる。

# 予測確率と実際のラベル
tibble(y_test, p_pred) %>% ggplot(aes(x = p_pred, fill = y_test)) + 
  geom_histogram(bins = 20, alpha = .5)

# 閾値=0.7
y_pred <- ifelse(p_pred < 0.7, "no", "yes")

tibble(
  precision = Precision(y_test, y_pred),
  recall = Recall(y_test, y_pred)
)

## # A tibble: 1 × 2
##   precision recall
##       <dbl>  <dbl>
## 1     0.696  0.974

# 閾値=0.05
y_pred <- ifelse(p_pred < 0.05, "no", "yes")

tibble(
  precision = Precision(y_test, y_pred),
  recall = Recall(y_test, y_pred)
)

## # A tibble: 1 × 2
##   precision recall
##       <dbl>  <dbl>
## 1         1 0.0261

Precision-Recall curve / average precision

閾値を様々な値に変えていった場合におけるPrecisionとRecallをプロットした曲線のことをPrecision-Recall curveという。

# （参考） Precision-Recall curve
# さまざまな閾値のもとでのPrecisionとRecallを計算する
metrics <- tibble() # 空のデータフレーム
for (threshold in seq(0, 1, 0.01)) {
  # ある閾値（threshold）のもとでの分類結果
  y_pred <- ifelse(p_pred < threshold, "no", "yes")
  # 1行だけのデータフレーム
  metric <- tibble(
    threshold,
    precision = Precision(y_test, y_pred),
    recall = Recall(y_test, y_pred)
  )
  metrics <- bind_rows(metrics, metric) # データフレームを結合して行を追加
}

# plot
metrics %>% filter(!is.na(precision)) %>% ggplot(aes(x = recall, y = precision)) +
  geom_line() + 
  labs(title="Precision-Recall curve")

分類モデルの性能が高いほど曲線の下側の面積（area under the curve: AUC）が大きくなる。この面積は\(K\)個の閾値の下でのrecall \(R_k\)とprecision \(P_k\)を用いて

\[ AP = \sum^K_{k=1} (R_k - R_{k-1}) P_k \]

のように計算されるもので、この指標はaverage precision (AP)やarea under the precision-recall curve (AUC-PR; PR-AUC)と呼ばれる。

# Area Under the Precision-Recall Curve
MLmetrics::PRAUC(p_pred, y_test)

## [1] 0.8331927

ROC-AUC

true positive rate\(TPR\)（recallの別名、陽性者を正しく陽性だと当てる率、sensitivityとも）とfalse positive rate\(FPR\)（陰性者のうち偽陽性になる率）

\[ TPR = \frac{TP}{P} = \frac{TP}{TP + FN} = \frac{Positiveを当てたもの}{Positiveのもの}\\ FPR = \frac{FP}{N} = \frac{FP}{FN + TN} = \frac{Positiveを外したもの}{Negativeのもの} \]

を用いて閾値を変えながら描いた曲線をreceiver operating characteristic（ROC; 受信者操作特性）曲線という。

# （参考） ROC curve
metrics <- tibble() # 空のデータフレーム
for (threshold in seq(0, 1, 0.01)) {
  # ある閾値（threshold）のもとでの分類結果
  y_pred <- ifelse(p_pred < threshold, "no", "yes")
  # 1行だけのデータフレーム
  metric <- tibble(
    threshold,
    tpr = MLmetrics::Sensitivity(y_test, y_pred),
    fpr = 1 - MLmetrics::Specificity(y_test, y_pred) # FPR = 1 - 特異度
  )
  metrics <- bind_rows(metrics, metric) # データフレームを結合して行を追加
}

# plot
metrics %>% filter(!is.na(tpr) & !is.na(fpr)) %>% ggplot(aes(x = fpr, y = tpr)) +
  geom_line() + 
  geom_abline(slope=1, linetype = "dashed") +
  labs(title="ROC curve", x="False Positive Rate", y="True Positive Rate")

この曲線の下側の面積をROC-AUCと呼ぶ（単にAUCと呼ばれることもある）。ランダムな予測（図中の点線）の下ではROC-AUCは0.5になるため、0.5との比較によって予測モデルの性能を評価できる。

y_test_int <- ifelse(y_test == "yes", 1, 0) # {0, 1}の整数型に変換
MLmetrics::AUC(p_pred, y_test_int) # ROC-AUC

## [1] 0.8259944

tibble(y_test, p_pred) %>% arrange(desc(p_pred))

## # A tibble: 209 × 2
##    y_test p_pred
##    <fct>   <dbl>
##  1 yes     0.948
##  2 yes     0.946
##  3 yes     0.945
##  4 yes     0.943
##  5 yes     0.943
##  6 yes     0.936
##  7 yes     0.934
##  8 yes     0.932
##  9 yes     0.920
## 10 yes     0.920
## # ℹ 199 more rows

4 参考文献

機械学習の教科書はたくさん出版されており、ウェブ上の解説も多い。各自いろいろ調べて、自分の関心やニーズに合うものを探してほしい。

計量経済学と機械学習の関係については、やや高度な内容だが、下記の文章に詳しい。

Katagiri(2020)計量経済学と機械学習の関係‒AI はさだめ, さだめは反事実‒

統計学や計量経済学では\(y\)は被説明変数、\(x\)は説明変数のように呼ばれるが、機械学習の文脈ではそれぞれ目的変数、特徴量と呼ばれる。同様に残差（residual）を誤差（error）や損失（loss）と呼ぶ。↩︎
なお、そのようなクラスごとの出現頻度の異なるデータは不均衡データ（imbalanced data）と呼ばれ、研究も多い（大崎、2022）↩︎

Rで学ぶ計量経済学と機械学習 10

機械学習１：線形回帰・ロジスティック回帰

安藤道人（立教大学）　三田匡能 (株式会社 GA technologies)

2023-12-15

1 人工知能と機械学習

1.1 人工知能とは

1.2 機械学習とは

1.3 学習の種類

1.4 機械学習のタスク

2 線形回帰

2.1 パラメータの推定方法

2.1.1 二乗誤差の最小化

2.1.2 推定量

2.2 Rでの実践

2.2.1 パッケージの読み込み

2.2.2 データの用意

2.2.3 線形回帰の学習

2.2.4 予測精度の評価

テスト用データに対する予測誤差

（参考）その他の評価指標

平均絶対誤差（mean absolute error: MAE）

決定係数（coefficient of determination）

3 ロジスティック回帰

3.1 パラメータの推定方法

3.1.1 例：コイン投げ

3.1.2 対数尤度と交差エントロピー誤差

3.1.3 ニュートン法

3.2 Rでの実践

3.2.1 パッケージの読み込み

3.2.2 データの用意

3.2.3 ロジスティック回帰の学習

3.2.4 予測精度の評価

混同行列

複数の閾値にわたって評価する指標

Precision-Recall curve / average precision

ROC-AUC

4 参考文献

Rで学ぶ計量経済学と機械学習 10 機械学習１：線形回帰・ロジスティック回帰

安藤道人（立教大学） 三田匡能 (株式会社 GA technologies)

2023-12-15

1 人工知能と機械学習

1.1 人工知能とは

1.2 機械学習とは

1.3 学習の種類

1.4 機械学習のタスク

2 線形回帰

2.1 パラメータの推定方法

2.1.1 二乗誤差の最小化

2.1.2 推定量

2.2 Rでの実践

2.2.1 パッケージの読み込み

2.2.2 データの用意

2.2.3 線形回帰の学習

2.2.4 予測精度の評価

テスト用データに対する予測誤差

（参考）その他の評価指標

平均絶対誤差（mean absolute error: MAE）

決定係数（coefficient of determination）

3 ロジスティック回帰

3.1 パラメータの推定方法

3.1.1 例：コイン投げ

3.1.2 対数尤度と交差エントロピー誤差

3.1.3 ニュートン法

3.2 Rでの実践

3.2.1 パッケージの読み込み

3.2.2 データの用意

3.2.3 ロジスティック回帰の学習

3.2.4 予測精度の評価

混同行列

複数の閾値にわたって評価する指標

Precision-Recall curve / average precision

ROC-AUC

4 参考文献

Rで学ぶ計量経済学と機械学習 10

機械学習１：線形回帰・ロジスティック回帰

安藤道人（立教大学）　三田匡能 (株式会社 GA technologies)