# パッケージのインストール
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")
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日)
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))
データと回帰モデル式から計画行列(説明変数の行列)を作成する。
# モデル候補
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) # カラム数
ニューラルネットワークでの学習を容易にするために説明変数を標準化する。 なお,標準化試験データは,訓練データの平均,標準偏差で標準化することに注意する。
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で置換
階層構造については,中間層の数を変えたり,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")) # 評価指標:絶対値誤差平均
バッチサイズは32(\(2^5\))が標準で多くのモデルにはこのままで対応できるが, モデルが複雑になればなるほど, 過学習防止や安定収束のため小さいバッチサイズ(16,32)が必要となる。 訓練データサイズが小さいときは, すべてのデータを使えるようにバッチサイズは訓練データサイズに合わせる。 訓練データサイズが大きいときは,64,128などを試す。 GPU利用時は,256,512などを試す。 このハイパーパラメータは試行錯誤で調整する必要がある。
エポック図において,訓練(training)ではMAEが低下し続ける(赤点)が, 評価(validation)では途中から値が増加(緑点)する場合がある。 これは過学習が発生していることを示している。 エポック数を過学習が発生する前までに設定する。
エポック数などのハイパーパラメータを決めるために 検証用データの訓練データに占める割合を設定する。 検証用データは訓練には使用しない。
本機能を使用すると,検証用データでの損失値(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) # コールバック設定
過学習する(エポック図で検証データでの損失が増加する)場合は, 次のような過学習対策を深層学習モデルに追加する。
layer_dropout(rate = 0.2)
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の仕様)で結果の再現性はないので注意
# 精度評価関数
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))
)
}
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
get.accuracy(yhat = yhat.lm, y = d.te$y)
次のデータで,坪単価を予測する回帰型深層学習モデルを構築せよ。 重回帰モデルと精度を比較せよ。
新台湾市の住宅価格
出典:【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