統計的機械学習の数理100問（09/20） - 「大人の教養・知識・気付き」を伸ばすブログ

　いい加減時代の潮流に乗ろうということで機械学習を学びたいと思う。またRはともかくとしてPythonは未経験であるため、丁度良い書籍として

統計的機械学習の数理100問 with R (with R) (機械学習の数理100問シリーズ)

作者:讓, 鈴木
共立出版

Amazon

統計的機械学習の数理100問 with Python (機械学習の数理100問シリーズ)

作者:讓, 鈴木
共立出版

Amazon

を用いることにする。

前回
2.　分類
- 2.3.　ロジスティック回帰のシミュレーション例
- 2.4　線形判別と二次判別
  - 2.4.1　Rでのシミュレーション
  - 2.4.2　Pythonでのシミュレーション
次回

前回

power-of-awareness.com

2.　分類

　有限個の数値を目的変数とするような分類の問題を取り扱う。

2.3.　ロジスティック回帰のシミュレーション例

　目的変数

$\begin{aligned} Y_i=\begin{cases}1,\\-1\end{cases} \end{aligned}$

に対して、標本数 $n=2000,\$ 説明変数は

$\begin{aligned} \boldsymbol{X}&={}^{t}(X_1,\cdots,X_n),\\ X_i&=\begin{cases}\varepsilon_i+1,&Y_i=1\\\varepsilon_i-1,&Y_i=-1\end{cases},\\ \varepsilon_i&\sim N(0,1),\\ X_i&\in\mathbb{R}^{p},p=1 \end{aligned}$

としてデータセット $\boldsymbol{Z}=\begin{bmatrix}\boldsymbol{1}_{n}&\boldsymbol{X}\end{bmatrix}$ を作成した。そして全データセットのうち半分に相当する $\displaystyle{\frac{n}{2}}=1000$ 個をランダムに選び、それを学習データ、選ばれなかったデータセットをテストデータとする。

############################################
### ロジスティック回帰のシミュレーション ###
############################################

# 説明変数を1つとするロジスティック回帰

### サンプルデータの作成
nm_sample <- 1000
vc_x <- c(rnorm(nm_sample)+1,rnorm(nm_sample)-1) # 説明変数に差異を与えるべく、±1する
vc_y <- c(rep(1,nm_sample),rep(-1,nm_sample))

# 学習データは全データの半分とする
vc_train <- sample(1:(2*nm_sample),nm_sample,replace = F) # 学習データとする行番号
df <- data.frame(vc_x,vc_y) # 全データセット

# 学習データを抽出する
mt_x_train <- as.matrix(df[vc_train,1])
mt_y_train <- as.vector(df[vc_train,2])

### モデルの構築
p <- 1 # 説明変数の数
mt_X_train <- cbind(1,mt_x_train) # 切片項のための1を追加
vc_beta <- 0
vc_gamma <- rnorm(p+1)

# Newton-Raphson法にて母数を最尤法で推定
while(sum((vc_beta - vc_gamma)^2)>0.0001){
  # スコアおよび
  vc_beta <- vc_gamma
  vc_s <- as.vector(mt_X_train %*% vc_beta) # s_i = X_iβ
  vc_v <- exp(-vc_s * mt_y_train) # exp(-y_i*s_i)
  
  # Jacobi行列を計算
  vc_u <- mt_y_train * vc_v /(1+vc_v)
  vc_w <- vc_v/(1+vc_v)^2
  mt_W <- diag(vc_w)
  vc_z <- vc_s + vc_u/vc_w
  
  vc_gamma <- as.vector(solve(t(mt_X_train) %*% mt_W %*% mt_X_train) %*% t(mt_X_train) %*% mt_W %*% vc_z)
  
  print(vc_gamma)
}

mt_X_test <- as.matrix(df[-vc_train,1])
mt_y_test <- as.matrix(df[-vc_train,2])

### 予測性能の評価
# 予測結果
vc_z_test <- 2 * as.integer((vc_beta[1] + mt_X_test * vc_beta[2])>0) - 1

# データと比較
tbl_results <- table(mt_y_test,vc_z_test)
print(tbl_results)
# 正答率
sum(diag(tbl_results))/sum(tbl_results) # 0.829

2.4　線形判別と二次判別

　観測データ $(\boldsymbol{X}_1,Y_1),\cdots,(\boldsymbol{X}_n,Y_n)\in\mathbb{R}^{p}\times\{-1,1\},p\in\mathbb{N}$ が存在し、それらから対応

$\begin{aligned} f:&\mathbb{R}^{p}\rightarrow\{-1,1\},\\ f:&\boldsymbol{X}_i\mapsto Y_i \end{aligned}$

を考えることとする。
　まず各 $y=\pm1$ の下での各 $\boldsymbol{X}_i\in\mathbb{R}^p$ の分布が $\mathcal{N}_p\left(\boldsymbol{\mu}_{j},\Sigma_{j}\right),j=\pm1$ であることが既知であると仮定する。このとき

$\begin{aligned} f_j(\boldsymbol{x})=\displaystyle{\frac{1}{\sqrt{(2\pi)^p \mathrm{det}\left(\Sigma_{j}\right)}}\exp\left\{-\frac{1}{2}{}^{t}(\boldsymbol{x}-\boldsymbol{\mu}_j)\Sigma_{j}^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_j)\right\}} \end{aligned}$

である。
　また $Y_i=j$ の生起する事前確率を $\pi_{j}$ としこれもまた既知であるとする。このときBayesの定理から

$\begin{aligned} P\{Y_i=j|\boldsymbol{X}_i=\boldsymbol{x}_i\}=\displaystyle{\frac{\pi_{j} f_j(\boldsymbol{x}_i)}{\pi_1 f_1(\boldsymbol{x}_i)+\pi_{-1}f_{-1}(\boldsymbol{x}_i)}} \end{aligned}$

と事後確率を評価できる。そして

$\begin{aligned} \displaystyle{\frac{\pi_{1} f_1(\boldsymbol{x}_i)}{\pi_1 f_1(\boldsymbol{x}_i)+\pi_{-1}f_{-1}(\boldsymbol{x}_i)}}\geq\displaystyle{\frac{\pi_{j} f_j(\boldsymbol{x}_i)}{\pi_{-1}f_{-1}(\boldsymbol{x}_i)+\pi_{-1}f_{-1}(\boldsymbol{x}_i)}}&\Longrightarrow Y_i=1,\\ \displaystyle{\frac{\pi_{1} f_1(\boldsymbol{x}_i)}{\pi_1 f_1(\boldsymbol{x}_i)+\pi_{-1}f_{-1}(\boldsymbol{x}_i)}}\lt\displaystyle{\frac{\pi_{j} f_j(\boldsymbol{x}_i)}{\pi_{-1}f_{-1}(\boldsymbol{x}_i)+\pi_{-1}f_{-1}(\boldsymbol{x}_i)}}&\Longrightarrow Y_i=-1, \end{aligned}$

もしくはこれと同値な

$\begin{aligned} \pi_{1} f_1(\boldsymbol{x}_i)\geq \pi_{-1} f_{-1}(\boldsymbol{x}_i)&\Longrightarrow Y_i=1,\\ \pi_{1} f_1(\boldsymbol{x}_i)\lt\pi_{-1} f_{-1}(\boldsymbol{x}_i)&\Longrightarrow Y_i=-1,\\ \end{aligned}$

を用いて推計すると誤り率を最小にできる*1。このような原理による判別は事後確率最大の原理という。
　2値でなくとも $K\geq2$ 値でも適用できる。説明変数の観測値 $\boldsymbol{x}$ があるときに、各 $k=1,\cdots,K$ について $y=k$ の確率が $P\{y=k|\boldsymbol{X}=\boldsymbol{x}\}$ であれば、推定値 $y=\hat{k}$ に対して $y=\hat{k}=k$ であるような確率は

$\begin{aligned} 1-\displaystyle{\sum_{k\neq \hat{k}}P\{y=k|\boldsymbol{X}=\boldsymbol{x}\}}=P\{y=\hat{k}|\boldsymbol{X}=\boldsymbol{x}\} \end{aligned}$

であるから、事後確率 $P\{y=k|\boldsymbol{X}=\boldsymbol{x}\}$ を最大にするような $k=\hat{k}$ とすることが誤り率最小の意味で最適である。
　以下では、 $K=2$ のときを考える。判定のための基準である

および

において $y=\pm1$ の境界線

$\begin{aligned} &-{}^{t}(\boldsymbol{x}-\boldsymbol{\mu}_{1})\Sigma_{1}^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_{1})+{}^{t}(\boldsymbol{x}-\boldsymbol{\mu}_{-1})\Sigma_{-1}^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_{-1})\\ =&\log\left(\displaystyle{\frac{\mathrm{det}\Sigma_{1}}{\mathrm{det}\Sigma_{-1}}}\right)-2\log\left(\displaystyle{\frac{\pi_{1}}{\pi_{-1}}}\right) \end{aligned}$

の性質を考える。まず左辺は $\boldsymbol{x}$ に関する二次式である。これを二次判別という。
　他方でもし $\Sigma_1=\Sigma_{-1}=\Sigma$ ならば、境界線が平面になる。これを線形判別という。実際、

$\begin{aligned} &-{}^{t}(\boldsymbol{x}-\boldsymbol{\mu}_{1})\Sigma^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_{1})+{}^{t}(\boldsymbol{x}-\boldsymbol{\mu}_{-1})\Sigma^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_{-1})\\ =&2(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{-1})\Sigma^{-1}\boldsymbol{x})\left({}^{t}\mu_1\Sigma^{-1}\mu_1-{}^{t}\mu_{-1}\Sigma^{-1}\mu_{-1}\right)\\ =&-2\log\left(\displaystyle{\frac{\pi_{1}}{\pi_{-1}}}\right) \end{aligned}$

が成り立つ。したがって $\pi_{1}=\pi_{-1}$ ならば $x=\displaystyle{\frac{\mu_1+\mu_{-1}}{2}}$ が境界になる*2。

2.4.1　Rでのシミュレーション

############
### 判別 ###
############

library("ggplot2")

### 標本を作成する
# 集団1の母数
mu_1 <- c(2,2)
sigma_1 <- diag(c(2,2))
rho_1 <- 0

# 集団2の母数
mu_2 <- c(-3,-3)
sigma_2 <- diag(c(1,1))
rho_2 <- 0.8

nm_sample <- 1000

### 標本の作成

# 集団1の標本を作成
vc_u1 <- rnorm(nm_sample)
vc_v1 <- rnorm(nm_sample)
vc_x1 <- sigma_1[1,1] * vc_u1 + mu_1[1]
vc_y1 <- (rho_1 * vc_u1 + sqrt(1 - rho_1^2) * vc_v1) * sigma_1[2,2] + mu_1[2]

# 集団2の標本を作成
vc_u2 <- rnorm(nm_sample)
vc_v2 <- rnorm(nm_sample)
vc_x2 <- sigma_2[1,1] * vc_u2 + mu_2[1]
vc_y2 <- (rho_2 * vc_u2 + sqrt(1 - rho_2^2) * vc_v2) * sigma_2[2,2] + mu_2[2]

# 判別分析(Linear discriminant analysis)
LDA <- function(vc_x, vc_mu, mt_inv, mt_de) {
  return(drop(-0.5 * t(vc_x - vc_mu) %*% mt_inv %*% (vc_x - vc_mu) - 0.5 * log(mt_de)))
}

### 実際のデータから母数を推定
mu_1_est <- mean(c(vc_x1, vc_y1))
mu_2_est <- mean(c(vc_x2, vc_y2))

# 
df1 <- data.frame(vc_x1, vc_y1)
mt <- cov(df1)
inv1 <- solve(mt)
de1 <- det(mt)

# 
df2 <- data.frame(vc_x2, vc_y2)
mt <- cov(df2)
inv2 <- solve(mt)
de2 <- det(mt)

# 母数の推定値を代入した判別境界値
LDA_1 <- function(u, v) {
  return(LDA(c(u, v), mu_1_est, inv1, de1))
}

LDA_2 <- function(u, v) {
  return(LDA(c(u, v), mu_2_est, inv2, de2))
}

# 事前確率の設定
nm_pi1 <- 0.5
nm_pi2 <- 0.5

# 標本とその評価値を図示
df_samples <- data.frame(group = rep("group1",length(vc_x1)),"x" = vc_x1, "y" = vc_y1)
df_samples <- rbind(df_samples,
                    data.frame(group = rep("group2",length(vc_x2)),"x" = vc_x2, "y" = vc_y2))

vc_u <- vc_v <- seq(-6,6,length = 100)
m <- length(vc_u)

# df_Wが正か負かで判断
df_W <- array(dim = c(m,m))

for(i in 1:m){
  for(j in 1:m){
    df_W[i,j] <- log(nm_pi1) + LDA_1(vc_u[i],vc_v[j]) - log(nm_pi2) - LDA_2(vc_u[i],vc_v[j])
  }
}

# 図示
contour(vc_u,vc_v,df_W,level = 1)
points(df_samples[df_samples[,1]=="group1",2],
       df_samples[df_samples[,1]=="group1",3],col = "red")
points(df_samples[df_samples[,1]=="group2",2],
       df_samples[df_samples[,1]=="group2",3],col = "blue")

サンプルデータと境界

2.4.2　Pythonでのシミュレーション

##############################
### 判別を人工データで行う ###
##############################

import numpy as np
from numpy import random as rd
import matplotlib.pyplot as plt
%matplotlib inline
import japanize_matplotlib

### 真の母数を設定
mu_1 = np.array([2,2])
sigma_1_1,sigma_1_2 = 2,2
rho_1 = 0

mu_2 = np.array([-3,-3])
sigma_2_1, sigma_2_2 = 1,1
rho_2 = 0.8

### 真の母数に基づきデータを生成
num_sample = 1000

# グループ1
lst_u = rd.randn(num_sample); lst_v = rd.randn(num_sample)
x_1 = sigma_1_1 * lst_u + mu_1[0]; y_1 = (rho_1 * lst_u + np.sqrt(1 - rho_1 **2) * lst_v) * sigma_1_2 + mu_1[1]

# グループ2
lst_u = rd.randn(num_sample); lst_v = rd.randn(num_sample)
x_2 = sigma_2_1 * lst_u + mu_2[0]; y_2 = (rho_2 * lst_u + np.sqrt(1 - rho_2 **2) * lst_v) * sigma_2_2 + mu_2[1]

### 
mu_1_est = np.average((x_1,y_1),1); mu_2_est = np.average((x_2,y_2),1)
df = np.array([x_1, y_1])
mt = np.cov(df, rowvar=1)
inv_1 = np.linalg.inv(mt)
de_1 = np.linalg.det(mt)

df = np.array([x_2, y_2])
mt = np.cov(df, rowvar=1)
inv_2 = np.linalg.inv(mt)
de_2 = np.linalg.det(mt)

### 母数の推定値を分布式に代入
def dst(x, mu, inv, de):
    return(-0.5 * (x - mu).T @ inv @ (x - mu) - 0.5 * np.log(de))
def dst_1(u, v):
    return dst(np.array([u,v]), mu_1, inv_1, de_1)
def dst_2(u, v):
    return dst(np.array([u,v]), mu_2, inv_2, de_2)

# 等高線を図示
pi_1, pi_2 = 0.5, 0.5

u = v = np.linspace(-8, 8, 100)
m = len(u)
w = np.zeros([m,m])

for i in range(m):
    for j in range(m):
        w[i,j] = np.log(pi_1) + dst_1(u[i], v[j]) - np.log(pi_2) - dst_2(u[i], v[j])

# 図示
plt.contour(u, v, w, levels = 1, colors = ['green'])
plt.scatter(x_1, y_1, c = "red")
plt.scatter(x_2, y_2, c = "blue")