# パッケージのインストール
if (!require(keras3))
{
  # 【事前準備】
  # Pythonパッケージtensorfow(テンサーフロー/テンソルフロー)をインストールする。
  # Terminalで「pip install tensorflow」を入力しインストール
  # tensorflowは,Python 3.9–3.12(2024-11-05現在)にのみ対応
  
  install.packages("keras3")
  library(keras3)
  install_keras() # r-kerasという仮想環境にKerasがインストールされる(初回だけ実行)。
  # install_keras(tensorflow = "gpu") # GPU演算用(NVIDIA GPUとcuDNNが必要)
}

if (!require(plotly)) install.packages("plotly")

1 データ

set.seed(5) # 乱数シード

n <- 24*7*2 # 時間数(全14日)
x <- 1:n    # 時刻
x2 <- x^2
x3 <- 100*sin(2*pi*x/24)

f <- function(x) 100 + 0.1*x + 0.02*x2 + x3

y <- f(x) + rnorm(n, sd = 5)

d <- data.frame(y, x, x2, x3)

d
ii <- 1:(24*11)  # 11日分のインデックス
d.tr <- d[ ii, ] # 訓練データを抽出(前11日)
d.te <- d[-ii, ] # 試験データを抽出(後3日)

1.1 グラフ

vline <- list(type = "line",
              line = list(color = "blue", dash = "dash"),
              x0 = ii[length(ii)], x1 = ii[length(ii)],
              y0 = 0, y1 = max(d$y))

plot_ly(type = "scatter", mode = "markers") |>
  add_trace(x = d$x, y = d$y,  mode = 'markers', name = "観測値") |>
  add_trace(x = x,   y = f(x), mode = 'lines',   name = "真値") |>
  layout(title = "青破線の左側:訓練領域,右側:予測領域", shapes = list(vline))

2 回帰型ニューラルネットワークのモデル作成

2.1 計画行列(Design Matrix)作成

データと回帰モデル式から計画行列(説明変数の行列)を作成する。

# モデル候補
FML <- vector("list", 4)
FML[[1]] <- formula("y ~ x")
FML[[2]] <- formula("y ~ x + x2")
FML[[3]] <- formula("y ~ x + x2 + x3")
FML[[4]] <- formula("y ~ poly(x, 4)")

# 選択モデル番号
ID.MODEL <- 3

# 回帰モデルから計画行列を作成
x.tr <- model.matrix(FML[[ID.MODEL]], d.tr)
x.te <- model.matrix(FML[[ID.MODEL]], d.te)
nc <- ncol(x.tr) # カラム数

2.2 データの標準化

ニューラルネットワークでの学習を容易にするために説明変数を標準化する。 なお,標準化試験データは,訓練データの平均,標準偏差で標準化することに注意する。

me.tr <- apply(x.tr, 2, mean) # 訓練データの説明変数ごとの平均
sd.tr <- apply(x.tr, 2, sd)   # 訓練データの説明変数ごとの標準偏差
x.tr <- scale(x.tr, center = me.tr, scale = sd.tr) # 標準化訓練データ
x.te <- scale(x.te, center = me.tr, scale = sd.tr) # 標準化試験データ
x.tr[, 1] <- 1 # 切片は標準化できない(NaNなる)ので1で置換
x.te[, 1] <- 1 # 切片は標準化できない(NaNなる)ので1で置換

2.3 ニューラルネットワークの階層構造作成

階層構造については,中間層の数を変えたり,unit(ニューロン)数を変えたりして, いくつか試してみることが必要。 中間層では入力層の数より小さいunit数を設定すると情報が失われるので注意。 通常2の累乗(32,64,126など)で設定する場合が多い。

# 深層学習モデル
model <- 
  keras_model_sequential(input_shape = c(nc))    |> # 入力層
  layer_dense(units = nc*2, activation = "relu") |> # 中間層(ReLU)
  layer_dense(units = nc*4, activation = "relu") |> # 中間層(ReLU)
  layer_dense(units = nc*2, activation = "relu") |> # 中間層(ReLU)
  layer_dense(units = 1,    activation = "linear")  # 出力層(線形)

# モデル概要
summary(model)
## Model: "sequential"
## ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━┓
## ┃ Layer (type)                      ┃ Output Shape             ┃       Param # ┃
## ┡━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━┩
## │ dense (Dense)                     │ (None, 8)                │            40 │
## ├───────────────────────────────────┼──────────────────────────┼───────────────┤
## │ dense_1 (Dense)                   │ (None, 16)               │           144 │
## ├───────────────────────────────────┼──────────────────────────┼───────────────┤
## │ dense_2 (Dense)                   │ (None, 8)                │           136 │
## ├───────────────────────────────────┼──────────────────────────┼───────────────┤
## │ dense_3 (Dense)                   │ (None, 1)                │             9 │
## └───────────────────────────────────┴──────────────────────────┴───────────────┘
##  Total params: 329 (1.29 KB)
##  Trainable params: 329 (1.29 KB)
##  Non-trainable params: 0 (0.00 B)
# 高速演算のためのコンパイル(PCが素早く理解できる機械語に翻訳)
compile(model,
        loss      = "mse", # mean-squared-erros(2乗誤差平均)
        optimizer = optimizer_adam(learning_rate = 0.1), # 最適化アルゴリズム 
        metrics   = c("mean_absolute_error")) # 評価指標:絶対値誤差平均

2.4 深層学習モデルのフィッティングと予測

  • バッチサイズ(batch_size)

バッチサイズは32(\(2^5\))が標準で多くのモデルにはこのままで対応できるが, モデルが複雑になればなるほど, 過学習防止や安定収束のため小さいバッチサイズ(16,32)が必要となる。 訓練データサイズが小さいときは, すべてのデータを使えるようにバッチサイズは訓練データサイズに合わせる。 訓練データサイズが大きいときは,64,128などを試す。 GPU利用時は,256,512などを試す。 このハイパーパラメータは試行錯誤で調整する必要がある。

  • エポック数

エポック図において,訓練(training)ではMAEが低下し続ける(赤点)が, 評価(validation)では途中から値が増加(緑点)する場合がある。 これは過学習が発生していることを示している。 エポック数を過学習が発生する前までに設定する。

  • 検証用データ割合(validation_split)

エポック数などのハイパーパラメータを決めるために 検証用データの訓練データに占める割合を設定する。 検証用データは訓練には使用しない。

  • 早期停止機能(callback_early_stopping)

本機能を使用すると,検証用データでの損失値(val_loss)がこれ以上低減しないと判断された場合, 計算を早期に停止させることができる。 なお,訓練データの損失値(loss)も設定できるが, エポック数について単調減少となる場合が多いので早期停止しにくい。 早期停止判定するためのエポック回数は忍耐数(patience)で指定する。 本機能を使用する場合は,エポック数(epochs)は大きめ, 少なくとも忍耐数以上に設定しておく。 忍耐数はあまり小さいと初期停滞期に早期停止してしまうことがある。

# コールバック設定
callbacks <- list(
  # 早期停止(検証データでの損失値の改善が20エポック以上なかったら停止)
  callback_early_stopping(patience = 20, monitor = "val_loss"),
  
  # 検証データでの損失が改善されない限りモデルを上書きしない設定
  # (early_stoppingとセットで使用する)
  callback_model_checkpoint(filepath = "bestmodel.keras",
                            monitor = "val_loss", save_best_only = T),
  
  # 検証データでの損失が改善せず停滞した時(判定:5エポック)
  # に局所解を抜け出すため学習率を0.1倍に下げる設定。 
  callback_reduce_lr_on_plateau(monitor = "val_loss", 
                                factor = 0.1, patience = 5)
)

# フィッティング
fit.dp <- fit(model,                  # 深層学習モデル
              x.tr,                   # 計画行列 
              d.tr$y,                 # 目的変数
              verbose    = 0,         # 1:出力表示(低速),0:出力表示抑制
              batch_size = 2^6,       # バッチサイズ(要調整)
              epochs     = 100,       # エポック数(早期停止設定時設定不要) 
              validation_split = 0.2, # 検証用データ割合(訓練には不使用)
              callbacks  = callbacks) # コールバック設定

2.4.1 エポック図

過学習する(エポック図で検証データでの損失が増加する)場合は, 次のような過学習対策を深層学習モデルに追加する。

  • ドロップアウト層を追加

layer_dropout(rate = 0.2)

  • layer_denseに学習重みの正則化のオプション「kernel_regularizer」を追加

layer_dense(…, kernel_regularizer = regularizer_l2(0.001))

plot(fit.dp)

# 予測
yhat <- predict(model, x.te)
## 3/3 - 0s - 33ms/step
# 重回帰モデルとの比較
fit <- lm(FML[[ID.MODEL]], data = d.tr)
yhat.lm <- predict(fit, newdata = d.te)

# 予測結果のグラフ
plot_ly(type = "scatter", mode = "markers") |>
  add_trace(x = d$x,    y = d$y,     mode = "markers", name = "観測値") |>
  add_trace(x = d.te$x, y = yhat,    mode = "markers", name = "予測値(DNN)") |>
  add_trace(x = d.te$x, y = yhat.lm, mode = "lines",   name = "予測値(重回帰)") |>
  layout(shapes = list(vline))

計算のたびに生成される乱数の影響(Kerasの仕様)で結果の再現性はないので注意

  • 何度か訓練し性能の良かったモデルを採用する。
  • 正確な精度評価は交差検証法(CV)を利用する。

2.4.2 精度評価

# 精度評価関数
get.accuracy <- function(yhat, y)
{
  data.frame(MBE  = mean(yhat - y),
             MAE  = mean(abs(yhat - y)),
             MAPE = mean(abs((yhat - y) / y)) * 100,
             RMSE = sqrt(mean((yhat - y)^2))
  )
}

3 精度指標(DNN)

get.accuracy(yhat = yhat, y = d.te$y)
# Kerasパッケージでの精度指標
evaluate(model, x.te, d.te$y)
## 3/3 - 0s - 8ms/step - loss: 62411.9648 - mean_absolute_error: 240.3165
## $loss
## [1] 62411.96
## 
## $mean_absolute_error
## [1] 240.3165

4 精度指標(重回帰)

get.accuracy(yhat = yhat.lm, y = d.te$y)

5 演習課題

次のデータで,坪単価を予測する回帰型深層学習モデルを構築せよ。 重回帰モデルと精度を比較せよ。

5.1 データ

新台湾市の住宅価格

出典:【UCI Machine Learning Repository】Real estate valuation data set Data Set

説明変数 内容
id 連番
yr 住宅購入年
yrs_old 築年数
m_sta 最寄り駅までの距離(m)
nstores 近隣コンビニエンスストア数
lat 緯度(度)
lon 経度(度)
price 坪単価(万新台湾$)
d <- read.csv('https://stats.dip.jp/01_ds/data/real_estate_price.csv')
n <- nrow(d)

# 要約統計量
summary(d)
##        id              yr          yrs_old           m_sta        
##  Min.   :  1.0   Min.   :2012   Min.   : 0.000   Min.   :  23.38  
##  1st Qu.:104.2   1st Qu.:2012   1st Qu.: 9.025   1st Qu.: 289.32  
##  Median :207.5   Median :2013   Median :16.100   Median : 492.23  
##  Mean   :207.5   Mean   :2013   Mean   :17.713   Mean   :1083.89  
##  3rd Qu.:310.8   3rd Qu.:2013   3rd Qu.:28.150   3rd Qu.:1454.28  
##  Max.   :414.0   Max.   :2013   Max.   :43.800   Max.   :6488.02  
##     nstores            lat             lon            price       
##  Min.   : 0.000   Min.   :24.93   Min.   :121.5   Min.   :  7.60  
##  1st Qu.: 1.000   1st Qu.:24.96   1st Qu.:121.5   1st Qu.: 27.70  
##  Median : 4.000   Median :24.97   Median :121.5   Median : 38.45  
##  Mean   : 4.094   Mean   :24.97   Mean   :121.5   Mean   : 37.98  
##  3rd Qu.: 6.000   3rd Qu.:24.98   3rd Qu.:121.5   3rd Qu.: 46.60  
##  Max.   :10.000   Max.   :25.01   Max.   :121.6   Max.   :117.50
set.seed(7)
ii <- sample(1:n, size = floor(0.8*n))
d.tr <- d[ ii, -1] # 訓練データ(idを除く)
d.te <- d[-ii, -1] # テストデータ(idを除く)
d.tr