「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。目下、データ分析・語学に力点を置いています。

MENU

統計的機械学習の数理100問(09/20)

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

を用いることにする。

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のときを考える。判定のための基準である


\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}

および


\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}

において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")


サンプルデータと境界

*1:ただしf_{\pm}正規分布でその母数(ベクトル)が既知、更に事前確率が既知であるという前提の下で成り立つ。実際には母数(ベクトル)としてその推定値を用いる。

*2:もし\pi_{\pm1},f_{\pm1}が未知ならば、観測データから推定する。

プライバシーポリシー お問い合わせ