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

一流の大人（ビジネスマン、政治家、リーダー…）として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。データ分析・語学に力点を置いています。　→現在、コンサルタントの雛になるべく、少しずつ勉強中です(※2024年1月21日改訂)。

MENU

統計的機械学習の数理100問（17/20）

機械学習 統計学 Python R

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

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

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

Amazon

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

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

Amazon

を用いることにする。

前回

power-of-awareness.com

7.　決定木

7.1　回帰の決定木

　説明変数(p個)と目的変数の関係を、観測データ(x_1,y_1),\cdots,(x_n,y_n)\in\mathbb{R}^p\times\mathbb{R}から推定し決定木を推定する。
　決定木は頂点と枝からなり、枝が左右に分岐する頂点を分岐点(または内点)、分岐しない頂点を端点という。枝で隣接する2頂点のうち端点に近い頂点を子、遠い頂点を親という。親を持たない頂点を根という。
　決定木が構成されると、各\boldsymbol{x}\in\mathbb{R}^pは端点に対応する領域R_1,\cdots,R_mのいずれかに属することになる。すなわち同時密度関数f_{XY}(\boldsymbol{x},y)に対して


\begin{aligned} \bar{y}_j=E[Y|R_j]=\displaystyle{\frac{\displaystyle{\int_{-\infty}^{\infty}\int_{R_j}y f_{XY}(x,y)dxdy} }{\displaystyle{\int_{-\infty}^{\infty}\int_{R_j}f_{XY}(x,y)dxdy}}} \end{aligned}


として対応


\begin{aligned} x_i\in R_j\Rightarrow\hat{y}_i=\bar{y}_j \end{aligned}


を定め、


\begin{aligned} \displaystyle{\int_{-\infty}^{\infty}\sum_{j=1}^{m}\int_{R_j}(y-\bar{y}_j)^2 f(x,y)dxdy} \end{aligned}


を最小にするようにm\geq1および領域R_1,\cdots,R_mを決める。

7.1　回帰の決定木における論点

　同時確率密度関数f_{XY}は一般には未知で、標本から推定する必要がある。そのため、x_i\in R_jであるようなiの個数をn_jとすれば、


\begin{aligned} \bar{y}_j=\displaystyle{\frac{1}{n_j}\sum_{i\ \mathrm{s,t,}\ x_i\in R_j} y_i} \end{aligned}


に置き換えて


\begin{aligned} x\in R_j\Rightarrow \hat{y}_i=\bar{y}_j \end{aligned}


というルールを定め、


\begin{aligned} \displaystyle{\sum_{j=1}^{m}\sum_{i\ \mathrm{s,t,}\ x_i\in R_j}(y_i-\bar{y}_j)^2} \end{aligned}


を最小にするようにm\geq1および領域R_1,\cdots,R_mを決めることになる。
　しかし上式を最小にするように領域を定めると、領域数mを大きくすればするほどこの値は小さくなる。そのため極端にはm=Nとすれば0となるのであって、過学習を危惧しなければならず、クロス・バリデーションのように訓練データとテストデータとを分ける措置が必要である。
　また過学習を積極的に防ぐべく\alpha\gt0として、訓練データから


\begin{aligned} \displaystyle{\sum_{j=1}^{m}\sum_{i\ \mathrm{s,t,}\ x_i\in R_j}(y_i-\bar{y}_j)^2}+\alpha m \end{aligned}


と罰則項を含めて、これを最小にするm\geq1およびR_1,\cdots,R_mを決める方法がある。この\alphaはクロス・バリデーションで決める。たとえば

  1. 複数の\alpha\gt0の候補を用意する。
  2. 一定割合のデータを訓練データとして上式を最小にする決定木を求める。
  3. 訓練データにしなかったデータをテストデータとして性能を評価する。
  4. 2-3を複数回行いその算術平均を取ることで1つの\alphaに対する性能が求まる。
  5. 4.を繰り返すことですべての\alphaの性能を求めた後に性能が最も高い\alphaを決定する。

という流れを行うことになる。


　さらに各分岐点で変数を何個用いるか、分岐を何通りにするかなどの自由度がある。また本来は各分岐点で分岐に用いる変数をトップダウン的に決めても最適な決定木は得られない。決定木のすべての頂点を同時に見て、最適な変数の組み合わせを選択する必要がある。そこで最適解が得られないことが分かっていても、根から各端点までトップダウン的に分岐に用いる変数を1個(X_j)とし、その変数がある値以上か否かという閾値(しきいち)を定めることが多い。他には、ある頂点にサンプルが一定数（n_{\mathrm{min}}）以上無いとそれ以上分割しないというルールを定めることが多い。以下では標本x_1,\cdots x_nの部分集合x_i,i\in Sが残っており、これを2つの部分集合に更に分割する場合、X_jがx_{i,j}以上というルールの中で最適なi,jを定める。以下では


\begin{aligned} \displaystyle{\sum_{k\ \mathrm{s.t.}\ x_{k.j}\lt x_{i,j}}(y_i-\bar{y}_{i,j}^{L})^2}+\displaystyle{\sum_{k\ \mathrm{s.t.}\ x_{k.j}\geq x_{i,j}}(y_i-\bar{y}_{i,j}^{R})^2} \end{aligned}


を最小にするようなi,jを選ぶ。ただしn_L,n_Rをそれぞれx_{k,j}\lt x_{i,j}を満たすようなkの個数、x_{k,j}\geq x_{i,j}を満たすようなkの個数である。また


\begin{aligned} \bar{y}_{i,j}^{L}&=\displaystyle{\frac{1}{n_L}\sum_{k\ \mathrm{s.t.}\ x_{k.j}\lt x_{i,j}}y_i},\\ \bar{y}_{i,j}^{R}&=\displaystyle{\frac{1}{n_R}\sum_{k\ \mathrm{s.t.}\ x_{k.j}\geq x_{i,j}}y_i} \end{aligned}


と定義した。

7.2　分類の決定木

　分類の場合の決定木では、同じ端点に含まれるサンプルに対して同じクラスを割り当てる。このときその領域で最も蓋然性の高いクラスを選ぶことで誤り率が最小になる。そこでp個の説明変数の値からY=1,2,\cdots,Kのいずれかの値に対応させる場合に、その同時確率f_{XY}(x,k)が与えられている場合は、


\begin{aligned} \displaystyle{\frac{\displaystyle{\int_{R_j}f_{XY}(x,k)dx}}{\displaystyle{\sum_{h=1}^{K}\int_{R_j}f_{XY}(x,h)dx}}} \end{aligned}


を最大にするようなkを\bar{y}_jとして


\begin{aligned} x_i\in R_j\Rightarrow \hat{y}_i=\bar{y}_j \end{aligned}


というルールを設定することで平均誤り率


\begin{aligned} \displaystyle{\sum_{k=1}^{K}\sum_{j=1}^{m}\int_{R_j}\boldsymbol{1}_{\{\bar{y}_j\neq k\}}f_{XY}(x,k)dx} \end{aligned}


を最小にさせる。
　ただし実際にはx_i\in R_jかつy_i=kとなるようなiの中で頻度が最大のkを\bar{y}_jとおくとき、


\begin{aligned} x_i\in R_j\Rightarrow \hat{y}_i=\bar{y}_j \end{aligned}


という決定を行い、


\begin{aligned} \displaystyle{\sum_{j=1}^{m}\sum_{i\ \mathrm{s.t.}\ x_i\in R_j}\boldsymbol{1}(y_i\neq \bar{y}_k)} \end{aligned}


を最小にするようにm\geq1および領域R_1,\cdots,R_mを決定する。
　分類においても過学習に注意しなければならない。また分類の決定木を生成する際にどのような基準で分岐を行うかが重要である。


\begin{aligned} \displaystyle{\sum_{j=1}^{m}\sum_{i\ \mathrm{s.t.}\ x_i\in R_j}\boldsymbol{1}(y_i\neq \bar{y}_k)} \end{aligned}


は領域R_jにクラスY=kがn_{j,k}個あるとすれば、n_j=\displaystyle{\sum_{k=1}^{K}n_{j,k} }とおくと、この値は\displaystyle{\sum_{j=1}^{m}\left(n_j-\displaystyle{\max_{k}n_{j,k}}\right)}と書けるから、\hat{p}_{j,k}=\displaystyle{\frac{n_{j,k}}{n_j}}とすれば各領域R_jで誤り率


\begin{aligned} E_j=1-\displaystyle{\max_{k}\hat{p}_{j,k}} \end{aligned}


を最小にすればよいことになる、
　しかしサンプル集合の分割途中で端点まで遠い場合やKの値が大きい場合、誤り率Eではなく\mathrm{Gini}指標


\begin{aligned} G_j=\displaystyle{\sum_{k=1}^{K}\hat{p}_{j,k}(1-\hat{p}_{j,k}) } \end{aligned}


やエントロピーを最小化する場合もある。

7.3　バギング

　バギングはブートストラップと同じ考え方を決定木の生成に応用したものである。同じデータセットから同じ行数を重複を許しつつ無作為に選び、それを用いて決定木を生成する。この操作をB回繰り返して決定木\hat{f}_1,\cdots,\hat{f}_Bを得る。ここでこれらはそれぞれ回帰または分類を行う関数の形式を取り、回帰ならば実数値の出力、分類であれば事前に定めたクラスのいずれかの値を返すものとする。
　新しい入力\boldsymbol{x}\in\mathbb{R}^pについて出力\hat{f}_1(\boldsymbol{x}),\cdots,\hat{f}_B(\boldsymbol{x})が得られるが、回帰であればその算術平均、分類であればB個の出力の中で最も頻度が高い値を取る。このような処理をバギングという。

7.4　ランダムフォレスト

　生成された決定木の相関性を考慮していない点がバギングの課題であった。そこで決定木の分岐で用いる変数候補をm\lt pに限定し、そのm変数をランダムに選んでその中から最適な変数を選択する。

7.5　ブースティング

　決定木におけるブースティングは、訓練データの目的変数の数に合うように目的変数の各値を\boldsymbol{y}\in\mathbb{R}^nとして、\boldsymbol{r}_0=\boldsymbol{y}={}^{t}(r_{0,1},\cdots,r_{0,n})とおく。また適当な\lambda\gt0を用意し、分岐の個数dを制限し、以下のように木を逐次的に生成する：

  • {}^{t}(\hat{f}_1(x_1),\cdots,\hat{f}_1(x_n))が\boldsymbol{r}_0に近くなるように木\hat{f}_1を生成する。
  • 木\hat{f}_1を用いて
    \begin{aligned}r_{1,1}=r_{0,1}-\lambda\hat{f}_1(x_1),\cdots,r_{1,n}=r_{0,n}-\lambda\hat{f}_1(x_n)\end{aligned}
    として\boldsymbol{r}_1={}^{t}(r_{1,1},\cdots,r_{1,n})と\boldsymbol{r}を更新する。
  • これをB回繰り返して\boldsymbol{r}=\boldsymbol{r}_Bと更新する。
  • 以上で得られた\hat{f}_1,\cdots,\hat{f}_Bを用いて、最終的な関数
    \begin{aligned}\hat{f}(\cdots)=\displaystyle{\lambda\sum_{b=1}^{B}\hat{f}_b(\cdot)}\end{aligned}
    が得られる。

7.6　Rによるシミュレーション

7.6.1　決定木
##############
### 決定木 ###
##############

sq_loss <- function(y){
  y_bar <- mean(y,na.rm = F)
  return(sum((y - y_bar)^2))
}

branch <- function(x, y, f, S, m = ncol(x)) {  # Sはサンプル集合の添え字
  n <- length(S)
  p <- ncol(x)
  if (n == 0){
    return(NULL)
  }

  best_score <- Inf
  for (j in 1:p) {
    for (i in S) {
      left <- NULL
      right <- NULL
      
      for (k in S) {
        if (x[k, j] < x[i, j]) {
          left <- c(left, k)
        } else {
          right <- c(right, k)
        }
      }
      L <- f(y[left])
      R <- f(y[right])
      score <- L + R
      if (score < best_score) {
        best_score <- score
        info <- list(i = i, j = j, left = left, right = right,
                     score = best_score, left_score = L, right_score = R)
      }
    }
  }
  return(info)
}

###
dt <- function(x, y, f = "sq_loss", alpha = 0, n_min = 1, m = ncol(x)) {
  if (f == "sq_loss") {
    g <- sq_loss
  } else if (f == "mis_match") {
    g <- mis_match
  } else if (f == "gini") {
    g <- gini
  } else {
    g <- entropy
  }
  
  n <- length(y)
  stack <- list()
  stack[[1]] <- list(parent = 0, set = 1:n, score = g(y))
  vertex <- list()
  k <- 0
  
  while (length(stack) > 0) {
    r <- length(stack)
    node <- stack[[r]]
    stack <- stack[-r]
    k <- k + 1
    res <- branch(x, y, g, node$set, m)
    if (node$score - res$score < alpha || length(node$set) < n_min ||
        length(res$left) == 0 || length(res$right) == 0) {
      vertex[[k]] <- list(parent = node$parent, j = 0, set = node$set)
    } else {
      vertex[[k]] <- list(parent = node$parent, set = node$set,
                          th = x[res$i, res$j], j = res$j)
      stack[[r]] <- list(parent = k, set = res$right,
                         score = res$right_score)
      stack[[r + 1]] <- list(parent = k, set = res$left,
                             score=res$left_score)
    }
  }
  
  # 最頻値
  mode <- function(y) {
    return(names(sort(table(y), decreasing = TRUE))[1])
  }
  
  # 端点にその値(center)，内点にその左右の子のID (left, right)
  r <- length(vertex)
  for (h in 1:r) {
    vertex[[h]]$left <- 0
    vertex[[h]]$right <- 0
  }
  for (h in r:2) {
    pa <- vertex[[h]]$parent
    if (vertex[[pa]]$right == 0) {
      vertex[[pa]]$right <- h
    } else {
      vertex[[pa]]$left <- h
    }
  }
  if (f == "sq_loss") {
    g <- mean
  } else {
    g <- mode
  }
  for (h in 1:r)
    if (vertex[[h]]$j == 0)
      vertex[[h]]$center <- g(y[vertex[[h]]$set])
  return(vertex)
}

############################
### Bostonを用いた決定木 ###
############################
### 変数12：黒人割合が倫理面で問題視されている点に注意

library("MASS")
library("igraph")

data(Boston)

x <- as.matrix(Boston[,1:13])
y <- as.vector(Boston[,14])

vertex <- dt(x, y, n_min = 50)

# グラフの出力
r <- length(vertex)
col <- array(dim =r)
edge_list <- matrix(nrow = r, ncol = 2)

for (h in 1:r){
  col[h] <- vertex[[h]]$j
}

for (h in 1:r){
  edge_list[h, ] <- c(vertex[[h]]$parent, h)
  
}

edge_list <- edge_list[-1, ]
g <- graph_from_edgelist(edge_list)
V(g)$color <- col
plot(g, layout = layout.reingold.tilford(g, root = 1))


# 変数表の出力
VAR <- NULL
TH <- NULL

for (h in 1:r) {
  if (vertex[[h]]$j != 0) {
    j <- vertex[[h]]$j
    th <- vertex[[h]]$th
    VAR <- c(VAR, j)
    TH <- c(TH, th)
  }
}

node <- as.data.frame(cbind(VAR, TH))

node <- data.frame(nodeID = 1:18,
                   var = c("町別犯罪率","25,000平方フィートを超える区画に分類される住宅地割合",
                           "町別非小売業割合",
                           "チャールズ川に隣接しているか",
                           "NOx濃度[単位:0.1ppm]",
                           "1戸当たり平均部屋数",
                           " 1940年以前建築持家割合",
                           "ボストン雇用センター(5つ)までの加重距離",
                           "主要高速道路へのアクセス性の指標",
                           "1万ドル当たり固定資産税率",
                           "町別生徒教師比率",
                           "町ごとの黒人の割合",
                           "低所得者人口の割合",
                           "住宅価格中央値[単位:千ドル]")[VAR],
                   thr = as.numeric(TH))
# IDは上・左→右からで番号付与


7.6.2　分類の決定木
library(igraph)


###
sq_loss <- function(y){
  y_bar <- mean(y,na.rm = F)
  return(sum((y - y_bar)^2))
}

branch <- function(x, y, f, S, m = ncol(x)) {  # Sはサンプル集合の添え字
  n <- length(S)
  p <- ncol(x)
  if (n == 0){
    return(NULL)
  }
  
  best_score <- Inf
  for (j in 1:p) {
    for (i in S) {
      left <- NULL
      right <- NULL
      
      for (k in S) {
        if (x[k, j] < x[i, j]) {
          left <- c(left, k)
        } else {
          right <- c(right, k)
        }
      }
      L <- f(y[left])
      R <- f(y[right])
      score <- L + R
      if (score < best_score) {
        best_score <- score
        info <- list(i = i, j = j, left = left, right = right,
                     score = best_score, left_score = L, right_score = R)
      }
    }
  }
  return(info)
}

dt <- function(x, y, f = "sq_loss", alpha = 0, n_min = 1, m = ncol(x)) {
  if (f == "sq_loss") {
    g <- sq_loss
  } else if (f == "mis_match") {
    g <- mis_match
  } else if (f == "gini") {
    g <- gini
  } else {
    g <- entropy
  }
  
  n <- length(y)
  stack <- list()
  stack[[1]] <- list(parent = 0, set = 1:n, score = g(y))
  vertex <- list()
  k <- 0
  
  while (length(stack) > 0) {
    r <- length(stack)
    node <- stack[[r]]
    stack <- stack[-r]
    k <- k + 1
    res <- branch(x, y, g, node$set, m)
    if (node$score - res$score < alpha || length(node$set) < n_min ||
        length(res$left) == 0 || length(res$right) == 0) {
      vertex[[k]] <- list(parent = node$parent, j = 0, set = node$set)
    } else {
      vertex[[k]] <- list(parent = node$parent, set = node$set,
                          th = x[res$i, res$j], j = res$j)
      stack[[r]] <- list(parent = k, set = res$right,
                         score = res$right_score)
      stack[[r + 1]] <- list(parent = k, set = res$left,
                             score=res$left_score)
    }
  }
  
  # 最頻値
  mode <- function(y) {
    return(names(sort(table(y), decreasing = TRUE))[1])
  }
  
  # 端点にその値(center)，内点にその左右の子のID (left, right)
  r <- length(vertex)
  for (h in 1:r) {
    vertex[[h]]$left <- 0
    vertex[[h]]$right <- 0
  }
  for (h in r:2) {
    pa <- vertex[[h]]$parent
    if (vertex[[pa]]$right == 0) {
      vertex[[pa]]$right <- h
    } else {
      vertex[[pa]]$left <- h
    }
  }
  if (f == "sq_loss") {
    g <- mean
  } else {
    g <- mode
  }
  for (h in 1:r)
    if (vertex[[h]]$j == 0)
      vertex[[h]]$center <- g(y[vertex[[h]]$set])
  return(vertex)
}

# 最頻値
mode <- function(y) {
  return(names(sort(table(y), decreasing = T))[1])
}

# 誤り率
mis_match <- function(y) {
  y_hat <- mode(y)
  return(sum(y != y_hat))
}

# Gini
gini <- function(y) {
  n <- length(y)
  if (n == 0)
    return(0)
  z <- as.vector(table(y))
  m <- length(z)
  T <- 0
  for (j in 1:m)
    T <- T + z[j] * (n - z[j])
  return(T/n)
}

# エントロピー
entropy <- function(y) {
  n <- length(y)
  if (n == 0)
    return(0)
  z <- as.vector(table(y))
  m <- length(z)
  T <- 0
  for (j in 1:m)
    if (z[j] != 0)
      T <- T + z[j] * log(n / z[j])
  return(T)
}

###

df <- iris
x <- as.matrix(df[, 1:4])
y <- as.matrix(df[, 5])
vertex <- dt(x, y, "mis_match", n_min = 4)
m <- length(vertex)
u <- NULL
v <- NULL
for (h in 1:m) {
  if (vertex[[h]]$j == 0) {
    w <- y[vertex[[h]]$set]
    u <- c(u, rep(mode(w), length(w)))
    v <- c(v, w)
  }
}
table(u, v)


col <- array(dim = m)
edge_list <- matrix(nrow = m, ncol = 2)
for (h in 1:m){
  col[h] <- vertex[[h]]$j
}
for (h in 1:m){
  edge_list[h, ] <- c(vertex[[h]]$parent, h)
}
edge_list <- edge_list[-1, ]


g <- graph_from_edgelist(edge_list)
V(g)$color <- col
V(g)$name <- ""
plot(g, pin = c(4, 4), layout = layout.reingold.tilford(g, root = 1))
title("誤り率")
7.6.3　バギング
par(mfrow = c(2, 4))

# データ生成
n <- 200
p <- 5
x <- matrix(rnorm(n * p), nrow = n, ncol = p)
beta <- rnorm(p)
y <- abs(round(x %*% beta + rnorm(n))) + 1

for (h in 1:8) {
  index <- sample(1:n, n, replace = TRUE)
  x <- x[index, ]
  y <- y[index]
  vertex <- dt(x, y, "mis_match", n_min = 6)
  r <- length(vertex)
  col <- array(dim = r)
  edge_list <- matrix(nrow = r, ncol = 2)
  for (h in 1:r){
    col[h] <- vertex[[h]]$j
  }
  
  for (h in 1:r){
    edge_list[h, ] <- c(vertex[[h]]$parent, h)
  }
  edge_list <- edge_list[-1, ]
  g <- graph_from_edgelist(edge_list)
  V(g)$color <- col
  V(g)$name <- ""
  plot(g, edge.arrow.size = 0.1, layout = layout.reingold.tilford(g, root = 1))
}
7.6.4　ランダムフォレスト
sq.loss <- function(y) {
  y.bar <- mean(y)
  return(sum((y - y.bar) ^ 2))
}

branch <- function(x, y, f, S, m = ncol(x)) {  # Sはサンプル集合の添え字
  n <- length(S)
  p <- ncol(x)
  if (n == 0)
    return(NULL)
  best.score <- Inf
  for (j in 1:p) {
    for (i in S) {
      left <- NULL
      right <- NULL
      for (k in S) {
        if (x[k, j] < x[i, j]) {
          left <- c(left, k)
        } else {
          right <- c(right, k)
        }
      }
      L <- f(y[left])
      R <- f(y[right])
      score <- L + R
      if (score < best.score) {
        best.score <- score
        info <- list(i = i, j = j, left = left, right = right,
                     score = best.score, left.score = L, right.score = R)
      }
    }
  }
  return(info)
}
dt <- function(x, y, f = "sq.loss", alpha = 0, n.min = 1, m = ncol(x)) {
  if (f == "sq.loss") {
    g <- sq.loss
  } else if (f == "mis.match") {
    g <- mis.match
  } else if (f == "gini") {
    g <- gini
  } else {
    g <- entropy
  }
  n <- length(y)
  stack <- list()
  stack[[1]] <- list(parent = 0, set = 1:n, score = g(y))
  vertex <- list()
  k <- 0
  while (length(stack) > 0) {
    r <- length(stack)
    node <- stack[[r]]
    stack <- stack[-r]  # POP
    k <- k + 1  # PUSHされた頂点のid = kは，POPされてはじめて付与
    res <- branch(x, y, g, node$set, m)  # 最適な分割に関する情報
    if (node$score - res$score < alpha || length(node$set) < n.min ||
        length(res$left) == 0 || length(res$right) == 0) {
      vertex[[k]] <- list(parent = node$parent, j = 0, set = node$set)
    } else {
      vertex[[k]] <- list(parent = node$parent, set = node$set,
                          th = x[res$i, res$j], j = res$j)
      stack[[r]] <- list(parent = k, set = res$right,
                         score = res$right.score)  # PUSH
      stack[[r + 1]] <- list(parent = k, set = res$left,
                             score=res$left.score)  # PUSH
    }
  }
  
  # 最頻値
  mode <- function(y) {
    return(names(sort(table(y), decreasing = TRUE))[1])
  }
  
  # 端点にその値(center)，内点にその左右の子のID (left, right)
  r <- length(vertex)
  for (h in 1:r) {
    vertex[[h]]$left <- 0
    vertex[[h]]$right <- 0
  }
  for (h in r:2) {
    pa <- vertex[[h]]$parent
    if (vertex[[pa]]$right == 0) {
      vertex[[pa]]$right <- h
    } else {
      vertex[[pa]]$left <- h
    }
  }
  if (f == "sq.loss") {
    g <- mean
  } else {
    g <- mode
  }
  for (h in 1:r)
    if (vertex[[h]]$j == 0)
      vertex[[h]]$center <- g(y[vertex[[h]]$set])
  return(vertex)
}



value <- function(u, vertex) {
  r <- 1
  while (vertex[[r]]$j != 0) {
    if (u[vertex[[r]]$j] < vertex[[r]]$th) {
      r <- vertex[[r]]$left
    } else {
      r <- vertex[[r]]$right
    }
  }
  return(vertex[[r]]$center)
}



# 最頻値
mode <- function(y) {
  return(names(sort(table(y), decreasing = TRUE))[1])
}

# 誤り率
mis.match <- function(y) {
  y.hat <- mode(y)
  return(sum(y != y.hat))
}

# Gini
gini <- function(y) {
  n <- length(y)
  if (n == 0)
    return(0)
  z <- as.vector(table(y))
  m <- length(z)
  T <- 0
  for (j in 1:m)
    T <- T + z[j] * (n - z[j])
  return(T/n)
}

# エントロピー
entropy <- function(y) {
  n <- length(y)
  if (n == 0)
    return(0)
  z <- as.vector(table(y))
  m <- length(z)
  T <- 0
  for (j in 1:m)
    if (z[j] != 0)
      T <- T + z[j] * log(n / z[j])
  return(T)
}



library(igraph)

# Giniとエントロピーも同様に実行した（mis.matchをgini, entropyに）

branch <- function(x, y, f, S, m = ncol(x)){  # mの値の設定, デフォルトはp
  n <- length(S)
  p <- ncol(x)
  if (n == 0)
    return(NULL)
  best.score <- Inf
  if (m < p) {                               # ここが違う
    T <- sample(1:p, m, replace = FALSE)     # 〃
  } else {                                   # 〃
    T <- 1:p                                 # 〃
  }                                          # 〃
  for (j in T) {  # Tの中で最適な変数を選ぶ
    for (i in S) {
      left <- NULL
      right <- NULL
      for (k in S) {
        if (x[k, j] < x[i, j]) {
          left <- c(left, k)
        } else {
          right <- c(right, k)
        }
      }
      L <- f(y[left])
      R <- f(y[right])
      score <- L + R
      if (score < best.score) {
        best.score <- score
        info <- list(i = i, j = j, left = left, right = right,
                     score = best.score, left.score = L, right.score = R)
      }
    }
  }
  return(info)
}


rf <- function(z) {
  zz <- array(dim = c(B, 50))
  zzz <- NULL
  for (b in 1:B) {
    for (i in 1:50)
      zz[b, i] <- (mode(z[1:b, i]) == y[i + 100])
    zzz <- c(zzz, sum(zz[b, ]))
  }
  return(zzz)
}
set.seed(101)
df <- iris
n <- nrow(df)
p <- ncol(df)
index <- sample(1:n, n, replace = FALSE)
x <- as.matrix(df[index, 1:4])
y <- as.vector(df[index, 5])
train <- 1:100
test <- 101:150
B <- 100
z <- array(dim = c(B, 50))

m <- 4
for (b in 1:B) {
  index <- sample(train, 100, replace = TRUE)
  vertex <- dt(x[index, ], y[index], "mis.match", n.min = 2, m = m)
  for (i in test)
    z[b, i - 100] <- value(x[i, ], vertex)
}
z4 <- z

# m = 4をm = 3, m = 2, m = 1に変えて，それぞれz3, z2, z1に格納
# この部分は書籍に不掲載
m <- 3
for (b in 1:B) {
  index <- sample(train, 100, replace = TRUE)
  vertex <- dt(x[index, ], y[index], "mis.match", n.min = 2, m = m)
  for (i in test)
    z[b, i - 100] <- value(x[i, ], vertex)
}
z3 <- z

m <- 2
for (b in 1:B) {
  index <- sample(train, 100, replace = TRUE)
  vertex <- dt(x[index, ], y[index], "mis.match", n.min = 2, m = m)
  for (i in test)
    z[b, i - 100] <- value(x[i, ], vertex)
}
z2 <- z

m <- 1
for (b in 1:B) {
  index <- sample(train, 100, replace = TRUE)
  vertex <- dt(x[index, ], y[index], "mis.match", n.min = 2, m = m)
  for (i in test)
    z[b, i - 100] <- value(x[i, ], vertex)
}
z1 <- z
# ここまで

plot(1:B, rf(z4) - 0.1, type = "l", ylim = c(0, 50), col = 2,
     xlab = "木の生成回数", ylab = "正答回数/50回",
     main = "ランダムフォレスト")
lines(1:B, rf(z3), col = 3)
lines(1:B, rf(z2) + 0.1, col = 4)
lines(1:B, rf(z1) - 0.1, col = 5)
legend("bottomright", legend = c("m = 4", "m = 3", "m = 2", "m = 1"),
       col = c(2, 3, 4, 5), lty = 1)



7.6.5　ブースティング
######################
### ブースティング ###
######################

library("ggplot2")

sq_loss <- function(y) {
  y_bar <- mean(y)
  return(sum((y - y_bar) ^ 2))
}

value <- function(u, vertex) {
  r <- 1
  while (vertex[[r]]$j != 0) {
    if (u[vertex[[r]]$j] < vertex[[r]]$th) {
      r <- vertex[[r]]$left
    } else {
      r <- vertex[[r]]$right
    }
  }
  return(vertex[[r]]$center)
}

# ランダムフォレスト
branch <- function(x, y, f, S, m = ncol(x)){  # mの値の設定, デフォルトはp
  n <- length(S)
  p <- ncol(x)
  if (n == 0)
    return(NULL)
  best.score <- Inf
  if (m < p) {                             ここが違う
    T <- sample(1:p, m, replace = FALSE)   〃
  } else {                                 〃
    T <- 1:p                               〃
  }                                        〃
  for (j in T) {  # Tの中で最適な変数を選ぶ
    for (i in S) {
      left <- NULL
      right <- NULL
      for (k in S) {
        if (x[k, j] < x[i, j]) {
          left <- c(left, k)
        } else {
          right <- c(right, k)
        }
      }
      L <- f(y[left])
      R <- f(y[right])
      score <- L + R
      if (score < best.score) {
        best.score <- score
        info <- list(i = i, j = j, left = left, right = right,
                     score = best.score, left.score = L, right.score = R)
      }
    }
  }
  return(info)
}

# 最頻値
mode <- function(y) {
  return(names(sort(table(y), decreasing = TRUE))[1])
}

# 誤り率
mis.match <- function(y) {
  y.hat <- mode(y)
  return(sum(y != y.hat))
}

# Gini
gini <- function(y) {
  n <- length(y)
  if (n == 0)
    return(0)
  z <- as.vector(table(y))
  m <- length(z)
  T <- 0
  for (j in 1:m)
    T <- T + z[j] * (n - z[j])
  return(T/n)
}


## ブースティング
b.dt <- function(x, y, d, f = "sq_loss") {
  n <- nrow(x)
  if (f == "sq_loss") {
    g <- sq_loss
  } else if (f == "mis.match") {
    g <- mis.match
  } else if (f == "gini") {
    g <- gini
  } else {
    g <- entropy
  }
  vertex <- list()
  vertex[[1]] <- list(parent = 0, set = 1:n, score = g(y), j = 0)
  while (length(vertex) <= 2 * d - 1) {                        
    r <- length(vertex)                                        
    gain.max <- -Inf                                           
    for (h in 1:r) {                                           
      if (vertex[[h]]$j == 0) {                                
        res <- branch(x, y, g, vertex[[h]]$set)                
        gain <- vertex[[h]]$score - res$score                  
        if (gain > gain.max) {                                 
          gain.max <- gain                                     
          h.max <- h                                           
          res.max <- res                                       
        }                                                      
      }                                                        
    }                                                          
    vertex[[h.max]]$th <- x[res.max$i, res.max$j]
    vertex[[h.max]]$j <- res.max$j
    vertex[[r + 1]] <- list(parent = h.max, set = res.max$left,
                            score = res.max$left.score, j = 0)
    vertex[[r + 2]] <- list(parent = h.max, set = res.max$right,
                            score = res.max$right.score, j = 0)
  }
  r <- 2 * d + 1                                               
  for (h in 1:r) {
    vertex[[h]]$left <- 0
    vertex[[h]]$right <- 0
  }
  for (h in r:2) {
    pa <- vertex[[h]]$parent
    if (vertex[[pa]]$right == 0) {
      vertex[[pa]]$right <- h
    } else {
      vertex[[pa]]$left <- h
    }
    if (vertex[[h]]$right == 0 & vertex[[h]]$left == 0)
      vertex[[h]]$j <- 0
  }
  
  mode <- function(y) {
    return(names(sort(table(y), decreasing = TRUE))[1])
  }
  
  if (f == "sq_loss") {
    g <- mean
  } else {
    g <- mode
  }
  for (h in 1:r)
    if (vertex[[h]]$j == 0)
      vertex[[h]]$center <- g(y[vertex[[h]]$set])
  return(vertex)
}


library(MASS)

df <- Boston
x <- as.matrix(df[, 1:13])
y <- as.matrix(df[, 14])
train <- 1:200
test <- 201:300
B <- 200
lambda <- 0.1

d <- 1
trees <- list()
r <- y[train]
for (b in 1:B) {
  trees[[b]] <- b.dt(x[train, ], r, d)
  for (i in train)
    r[i] <- r[i] - lambda * value(x[i, ], trees[[b]])
}
z <- array(0, dim = c(B, 600))
for (i in test)
  z[1, i] <- lambda * value(x[i, ], trees[[1]])
for (b in 2:B)
  for (i in test)
    z[b, i] <- z[b - 1, i] + lambda * value(x[i, ], trees[[b]])
out <- NULL
for (b in 1:B)
  out <- c(out, sum((y[test] - z[b, test]) ^ 2) / length(test))
out1 <- out

# d = 2, d = 3でも実行する
d <- 2
trees <- list()
r <- y[train]
for (b in 1:B) {
  trees[[b]] <- b.dt(x[train, ], r, d)
  for (i in train)
    r[i] <- r[i] - lambda * value(x[i, ], trees[[b]])
}
z <- array(0, dim = c(B, 600))
for (i in test)
  z[1, i] <- lambda * value(x[i, ], trees[[1]])
for (b in 2:B)
  for (i in test)
    z[b, i] <- z[b - 1, i] + lambda * value(x[i, ], trees[[b]])
out <- NULL
for (b in 1:B)
  out <- c(out, sum((y[test] - z[b, test]) ^ 2) / length(test))
out2 <- out

# d = 3
d <- 3
trees <- list()
r <- y[train]
for (b in 1:B) {
  trees[[b]] <- b.dt(x[train, ], r, d)
  for (i in train)
    r[i] <- r[i] - lambda * value(x[i, ], trees[[b]])
}
z <- array(0, dim = c(B, 600))
for (i in test)
  z[1, i] <- lambda * value(x[i, ], trees[[1]])
for (b in 2:B)
  for (i in test)
    z[b, i] <- z[b - 1, i] + lambda * value(x[i, ], trees[[b]])
out <- NULL
for (b in 1:B)
  out <- c(out, sum((y[test] - z[b, test]) ^ 2) / length(test))
out3 <- out


df_out <- data.frame(x = 21:100,
                     var = "d = 1",
                     y = out1[21:100])
df_out <- rbind(df_out,data.frame(x = 21:100,
                                  var = "d = 2",
                                  y = out2[21:100]))
df_out <- rbind(df_out,data.frame(x = 21:100,
                                  var = "d = 3",
                                  y = out3[21:100]))

### 結果を図示
g <- ggplot(df_out,aes(x=x,y=y,color=var)) + geom_line()
g <- g + labs(x = "生成した木の個数", y = "テストデータでの二乗誤差",title = "ブースティング",color = "")
g <- g + theme_classic()
g <- g + theme(plot.title = element_text(hjust = 0.5),axis.text = element_text(size = 10),
               axis.title.x = element_text(size = 12),
               axis.title.y = element_text(size = 12),
               legend.position = "bottom",
               legend.title = element_text(size = 10),
               legend.text = element_text(size = 10))
g <- g + geom_hline(yintercept = 0, linetype = 3)
plot(g)

7.7　Pythonのシミュレーション

7.7.1　回帰の決定木
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import japanize_matplotlib
import scipy
from numpy.random import randn #正規乱数


### 
def sq_loss(y):
    if len(y)==0:
        return 0
    else:
        y_bar = np.mean(y)
        return np.linalg.norm(y - y_bar)**2

def branch(x, y, S, rf=0):
    if rf==0:
        m=x.shape[1]
    if x.shape[0]==0:
        return([0,0,0,0,0,0,0])
    best_score=np.inf
    for j in range(x.shape[1]):
        for i in S:
            left=[]; right=[]
            for k in S:
                if x[k,j]<x[i,j]:
                    left.append(k)
                else:
                    right.append(k)
            left_score=f(y[left]); right_score=f(y[right])
            score=left_score+right_score
            if score < best_score:
                best_score=score
                i_1=i; j_1=j
                left_1=left; right_1=right 
                left_score_1=left_score; right_score_1=right_score
    return [i_1, j_1, left_1, right_1, best_score, left_score_1, right_score_1]

class Stack:
    def __init__(self, parent, set, score):
        self.parent = parent
        self.set = set
        self.score = score

class Node:
    def __init__(self, parent, j, th, set):
        self.parent = parent
        self.j = j
        self.th=th
        self.set=set
        

def dt(x, y, alpha=0, n_min=1, rf=0):
    if rf==0:
        m=x.shape[1]
    # 1個からなるstackを構成。決定木を初期化
    stack=[Stack(0, list(range(x.shape[0])), f(y))]  # 関数 fは、大域
    node=[]
    k=-1
    # stackの最後の要素を取り出して、決定木を更新する。
    while len(stack)>0:
        popped=stack.pop()
        k=k+1
        i, j, left, right, score, left_score, right_score=branch(x, y, popped.set,rf)
        if popped.score-score<alpha or len(popped.set)<n_min or len(left)==0 or len(right)==0:
            node.append(Node(popped.parent, -1, 0, popped.set))
        else:
            node.append(Node(popped.parent, j, x[i,j], popped.set))
            stack.append(Stack(k, right, right_score))  
            stack.append(Stack(k, left, left_score))
    # これより下でnode.left, node.rightの値を設定する
    for h in range(k,-1,-1):
        node[h].left=0; node[h].right=0; 
    for h in range(k,0,-1):
        pa=node[h].parent
        if node[pa].right==0:
            node[pa].right=h
        else:
            node[pa].left=h
    # これより下で、node.centerの値を計算する
    if f==sq_loss:
        g=np.mean
    else:
        g=mode_max
    for h in range(k+1):
        if node[h].j==-1:
            node[h].center=g(y[node[h].set])
        else:
            node[h].center=0
    return(node)


from sklearn.datasets import load_boston
boston = load_boston()
X=boston.data
y=boston.target
f=sq_loss
node=dt(X,y,n_min=50)

len(node) 

##############
### テスト ###
##############

from igraph import *
r=len(node)
edge=[]
for h in range(1,r):
    edge.append([node[h].parent,h])
TAB=[];
for h in range(r):
    if not node[h].j==0:
        TAB.append([h, node[h].j, node[h].th])
TAB

def draw_graph(node):
    r=len(node)
    col=[]
    for h in range(r):
        col.append(node[h].j)
    colorlist = ['#ffffff', '#fff8ff', '#fcf9ce', '#d6fada', '#d7ffff', '#d9f2f8', '#fac8be', '#ffebff','#ffffe0','#fdf5e6','#fac8be', '#f8ecd5', '#ee82ee']
    color=[colorlist[col[i]] for i in range(r)]
    edge=[]
    for h in range(1,r):
        edge.append([node[h].parent,h])
        g = Graph(edges=edge,directed=True)
        layout=g.layout_reingold_tilford(root=[0])
    out=plot(g,vertex_size=15,layout=layout,bbox=(300,300),vertex_label=list(range(r)), vertex_color=color)
    return(out)
draw_graph(node)
7.7.2　分類の決定木
def freq(y):
    y=list(y)
    return([y.count(i) for i in set(y)])
# モード(最頻度)
def mode(y):
    n=len(y)
    if n==0:
        return(0)
    return(max(freq(y)))
# 誤り率
def mis_match(y):
    return(len(y)-mode(y))
# Gini
def gini(y):
    n=len(y)
    if n==0:
        return(0)
    fr=freq(y)
    return(sum([fr[i]*(n-fr[i])/n for i in range(len(fr))]))
# Entropy
def entropy(y):
    n=len(y)
    if n==0:
        return(0)
    freq=[y.count(i) for i in set(y)]
    return(np.sum([-freq[i]*np.log (freq[i]/n) for i in range(len(freq))]))
def table_count(m,u,v):    # 再掲
    n=u.shape[0]
    count=np.zeros([m,m])
    for i in range(n):
        count[int(u[i]),int(v[i])]+=1
    return(count)
def mode_max(y):
    if len(y)==0:
        return(-np.inf)
    count = np.bincount(y)
    return(np.argmax(count))
from sklearn.datasets import load_iris
iris = load_iris()
iris.target_names
f=mis_match
n=iris.data.shape[0]
x=iris.data
y=iris.target
n=len(x)
node=dt(x,y,n_min=4)
m=len(node)
u=[]; v=[]
for h in range(m):
    if node[h].j==-1:
        w=y[node[h].set]
        u.extend([node[h].center]*len(w))
        v.extend(w)
table_count(3,np.array(u),np.array(v))
#sum([u[i]==v[i] for i in range(150)])

draw_graph(node)

from sklearn.datasets import load_iris
iris = load_iris()
iris.target_names
f=mis_match
index=np.random.choice(n, n, replace=False) # 並び替える
X=iris.data[index,:]
y=iris.target[index]
n_min_seq=np.arange(5,51,5)
s=15
for n_min in n_min_seq:
    SS=0
    for h in range(10):
        test=list(range(h*s,h*s+s))
        train=list(set(range(n))-set(test))
        node=dt(X[train,:],y[train], n_min=n_min)
        for t in test:
            SS=SS+np.sum(y[t]!=value(X[t,:],node))
    print(SS/n)
7.7.3　バギング
n=200
p=5
X=np.random.randn(n,p)
beta=randn(p)
Y=np.array(np.abs(np.dot(X,beta)+randn(n)),dtype=np.int64)
f=mis_match
node_seq=[]
for h in range(8):
    index=np.random.choice(n, n, replace=True) # 並び替える
    x=X[index,:]
    y=Y[index]
    node_seq.append(dt(x,y,n_min=6))
7.7.4　ランダムフォレスト
def branch(x, y, S, rf=0):                                  ##
    if rf==0:                                               ##
        T=np.arange(x.shape[1])                                        ##
    else:                                                   ##
        T=np.random.choice(x.shape[1], rf, replace=False)   ##
    if x.shape[0]==0:
        return([0,0,0,0,0,0,0])
    best_score=np.inf
    for j in T:                                             ##
        for i in S:
            left=[]; right=[]
            for k in S:
                if x[k,j]<x[i,j]:
                    left.append(k)
                else:
                    right.append(k)
            left_score=f(y[left]); right_score=f(y[right])
            score=left_score+right_score
            if score < best_score:
                best_score=score
                i_1=i; j_1=j
                left_1=left; right_1=right 
                left_score_1=left_score; right_score_1=right_score
    return [i_1, j_1, left_1, right_1, best_score, left_score_1, right_score_1]

def rf(z):
    z=np.array(z,dtype=np.int64)
    zz=[]
    for b in range(B):
        u=sum([mode_max(z[range(b+1),i])==y[i+100] for i in range(50)])
        zz.append(u)
    return(zz)
iris = load_iris()
iris.target_names
f=mis_match
n=iris.data.shape[0]
order=np.random.choice(n, n, replace=False) # 並び替える
X=iris.data[order,:]
y=iris.target[order]
train=list(range(100))
test=list(range(100,150))
B=100
plt.ylim([35, 55])
m_seq=[1,2,3,4]
c_seq=["r","b","g","y"]
label_seq=['m=1','m=2','m=3','m=4']
plt.xlabel('繰り返し数 b')
plt.ylabel('テスト50データでの正答数')
plt.title('ランダムフォレスト')
for m in m_seq:
    z=np.zeros((B,50))
    for b in range(B):
        index=np.random.choice(train, 100, replace=True)
        node=dt(X[index,:],y[index],n_min=2,rf=m)
        for i in test:
            z[b,i-100]=value(X[i,],node)
    plt.plot(list(range(B)),np.array(rf(z))-0.2*(m-2), label=label_seq[m-1], linewidth=0.8, c=c_seq[m-1])
plt.legend(loc='lower right')
plt.axhline(y=50,c="b",linewidth=0.5,linestyle = "dashed")
7.7.5　ブーステイング

def b_dt(x, y, d):
n=x.shape[0]
node=
first=Node(0, -1, 0, np.arange(n))
first.score=f(y[first.set])
node.append(first)
while len(node)<=2*d-1:
r=len(node)
gain_max=-np.inf
for h in range(r):
if node[h].j==-1:
i, j, left, right, score, left_score, right_score=branch(x, y, node[h].set)
gain=node[h].score-score
if gain >gain_max:
gain_max=gain
h_max=h
i_0=i; j_0=j
left_0=left; right_0=right
left_score_0=left_score; right_score_0=right_score
node[h_max].th=x[i_0,j_0]; node[h_max].j=j_0
next=Node(h_max, -1, 0, left_0)
next.score=f(y[next.set]); node.append(next)
next=Node(h_max, -1, 0, right_0)
next.score=f(y[next.set]); node.append(next)
r=2*d+1
for h in range(r):
node[h].left=0; node[h].right=0
for h in range(r-1,1,-1):
pa=node[h].parent
if node[pa].right==0:
node[pa].right=h
else:
node[pa].left=h
if node[h].right==0 and node[h].left==0:
node[h].j=-1
if f==sq_loss:
g=np.mean
else:
g=mode_max
for h in range(r):
if node[h].j==-1:
node[h].center=g(node[h].set)
# これより下でnode.left, node.rightの値を設定する
for h in range(r-1,-1,-1):
node[h].left=0; node[h].right=0;
for h in range(r-1,0,-1):
pa=node[h].parent
if node[pa].right==0:
node[pa].right=h
else:
node[pa].left=h
# これより下で、node.centerの値を計算する
if f==sq_loss:
g=np.mean
else:
g=mode_max
for h in range(r):
if node[h].j==-1:
node[h].center=g(y[node[h].set])
else:
node[h].center=0
return(node)
from sklearn.datasets import load_boston
boston = load_boston()
B=200
lam=0.1
X=boston.data
y=boston.target
f=sq_loss
train=list(range(200))
test=list(range(200,300))
# ブースティングの木をB個生成
# 各dで5分程度、合計15分程度かかります
trees_set=
for d in range(1,4):
trees=
r=y[train]
for b in range(B):
trees.append(b_dt(X[train,:],r,d))
for i in train:
r[i]=r[i]-lam*value(X[i,:],trees[b])
print(b)
trees_set.append(trees)

# テストデータで評価
out_set=
for d in range(1,4):
trees=trees_set[d-1]
z=np.zeros*1
for i in test:
z[0,i]=lam*value(X[i,],trees[0])
for b in range(1,B):
for i in test:
z[b,i]=z[b-1,i]+lam*value(X[i,:],trees[b])
out=[]
for b in range(B):
out.append(sum*2
out_set.append(out)
# グラフで表示
plt.ylim([0, 40])
c_seq=["r","b","g"]
label_seq=['d=1','d=2','d=3']
plt.xlabel('生成した木の個数')
plt.ylabel('テストデータでの二乗誤差')
plt.title('本書のプログラム (lambda=0.1)')
for d in range(1,4):
out=out_set[d-1]
u=range(20,100)
v=out[20:100];
plt.plot(u,v,label=label_seq[d-1], linewidth=0.8, c=c_seq[d-1])
plt.legend(loc='upper right')
|

次回

power-of-awareness.com

*1:B,600

*2:y[test]-z[b,test])**2)/len(test

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