rEDM を用いた Empirical Dynamic Modeling: (3) 近未来予測 (S-map)

近未来予測 (S-map)

1. S-map とは

S-map とは Sugihara (1994) によって提案された非線形時系列解析の方法で、近未来予測を行うことができます。Simplex projection は複数個の近傍点を特定して、それらの挙動に基づいて注目する点の未来を予測していました。一方 S-map では全ての点を利用し、その重要性を注目する点からの距離に基づいて重み付けして、かつ局所的な線形モデルを作成することで注目する点の未来を予測します。

Sugihara (1994) では、“Sequential locally weighted global linear maps: SLWGLM (or S-map)” として提案されています。EDM の枠組みにおいては SLWGLM として紹介されることはほぼなく、もっぱら S-map と呼ばれている風に思います。

2. S-map のアルゴリズム

rEDM のロードとデモデータの作成

S-map のアルゴリズムを再び rEDM パッケージ (Ye et al. 2015; Ye et al. 2018) とモデル時系列データを利用して解説します。再び自作の 2 種系のデータを作成しましょう。

# rEDM のロード
library(rEDM)
packageVersion("rEDM") # 2019.12月時点で v0.7.5

# 2種系のデータセットを作成
d0 <- as.data.frame(matrix(rep(NA,3*1000), ncol=3))
colnames(d0) <- c("time", "x", "y")
d0[,1] <- 1:1000
d0[1,2:3] <- c(0.1, 0.1)
for(i in 1:999){
        d0[i+1,2] <- d0[i,2]*(3.8 - 3.80*d0[i,2] - 0.02*d0[i,3])
        d0[i+1,3] <- d0[i,3]*(3.5 - 0.10*d0[i,2] - 3.50*d0[i,3])
}

# 最初の 100 time step のデータを捨てる
d <- d0[100:1000,]

前回と同じく $Y$ を 2 次元で埋め込み、$t = 240$ の点をターゲットとしてその未来を S-map で予測してみましょう。ターゲットの点を設定し、その点から他の点までの距離を計算するところまでは Simplex projection と同じですが、途中の局所線形モデルの計算部分でデータに NA が入っているとまずいのでそれらを除く処理を追加しています。

埋め込みとターゲットの点の指定

# 2次元で埋め込み
y_embed <- make_block(d$y, max_lag = 2)

# ターゲットの点を指定
t_target <- 240

# 予測の計算に使うデータを指定
valid_id <- (1:nrow(y_embed))[-t_target] # 自分自身の点は予測に使わないので除去
valid_id <- valid_id[-1] # 計算上問題なので NA の行を除去
valid_id <- valid_id[-(length(valid_id))] # 最後のデータはモデル作成に利用できないので除去

# t = 240 の点から各点への距離を計算
distance <- sqrt(rowSums((y_embed[valid_id,2:3] - c(y_embed[t_target, 2:3]))^2))

各点までの距離に基づいて重みを計算

次に各点までの距離に基づいて重みを計算します。ここで $\theta$ と呼ばれる重み付けの度合いを決める重要な変数を設定します。 この $\theta$ の値は Simplex projectin のときの $E$ と同様に、慎重に決める必要があるのですが、ここではひとまず $\theta = 1$ とおきます。

# 距離の平均値を計算
d_mean <- mean(distance, na.rm = TRUE)

# Nonlinear parameter と呼ばれる theta の値を指定
theta <- 1

# 距離と theta の値に基づき、重みを計算
weights <- exp(- theta * distance / d_mean)

以下のようにターゲット以外の全ての点に対して重みが割り当てられます。

> head(weights)
        2         3         4         5         6         7 
0.9118125 0.1266185 0.7191354 0.1789585 0.9603531 0.1455265 

この場合、状態空間上に重みがどのように分布しているか可視化してみましょう。今回は ggplot2 という可視化用のパッケージを使ってみます。まずは状態空間内のグリッド点を生成して、グリッド点までの距離を計算します。

# 状態空間上での重みの分布を可視化
# グリッド点を生成
x_grid <- seq(min(y_embed[valid_id,2])-0.05, max(y_embed[valid_id,2])+0.05, by=0.01)
y_grid <- seq(min(y_embed[valid_id,3])-0.05, max(y_embed[valid_id,3])+0.05, by=0.01)
x_axis <- rep(x_grid, length(y_grid))
y_axis <- NULL
for(i in 1:length(y_grid)) y_axis <- c(y_axis, rep(y_grid[i], length(x_grid)))
grid_xy <- as.matrix(cbind(x_axis,y_axis))

# グリッドごとにターゲットからの距離を計算
target_mat <- matrix(rep(as.numeric(y_embed[t_target,2:3]), length(x_axis)), nrow = length(x_axis), byrow = T)
distance2 <- sqrt(rowSums((grid_xy - target_mat)^2))

続いて、グリッド点毎の重みを計算して可視化します。白いバツ印がターゲットの点です。

# グリッド毎の重みを計算
grid_dist1 <- exp(- theta * distance2 / d_mean)
contour_weights1 <- data.frame(cbind(grid_xy, grid_dist1))
colnames(contour_weights1) <- c("col1", "col1_1", "weights")

# ggplot2 パッケージを用いて状態空間上の重みを可視化
library(ggplot2)
packageVersion("ggplot2") # v3.2.1

# 可視化
ggplot(y_embed[,2:3], aes(x = col1, y = col1_1)) +
        theme_classic() +
        geom_tile(data = contour_weights1, aes(x = col1, y = col1_1, fill = weights)) +
        geom_point(colour = "black", alpha = 0.5, size = 0.8) +
        scale_fill_gradient2(low = "gray90", mid = "gray90", high = "red3", midpoint = 0.2) +
        geom_point(aes(x = y_embed[t_target,2], y = y_embed[t_target,3]),
                   size = 2, colour = "white", shape = 4) +
        xlab("Y(t)") + ylab("Y(t-1)")

ちなみに $\theta=0$ のときの重みは全ての点で同じはずなので、図示すると以下のようになります。

# theta = 0 の場合
grid_dist0 <- exp(- 0 * distance2 / d_mean)
contour_weights0 <- data.frame(cbind(grid_xy, grid_dist0))
colnames(contour_weights0) <- c("col1", "col1_1", "weights")

ggplot(y_embed[,2:3], aes(x = col1_1, y = col1)) +
        theme_classic() +
        geom_tile(data = contour_weights0, aes(x = col1_1, y = col1, fill = weights)) +
        geom_point(colour = "black", alpha = 0.5, size = 0.8) +
        scale_fill_gradient2(low = "gray90", mid = "gray90", high = "red3", midpoint = 0.2) +
        geom_point(aes(x = y_embed[t_target,3], y = y_embed[t_target,2]),
                   size = 2, colour = "white", shape = 4) +
        xlab("Y(t-1)") + ylab("Y(t)")

つまり、S-map では $\theta$ の値が大きいほどターゲットの点に近い点の挙動を重要視して近未来予測を行うことになります。

特異値分解を用いて局所重み付け線形回帰を実行

さて、このようにして計算した各点の重みと特異値分解を利用して、最後に重み付け線形回帰を実行します。

# 重みに基づいて局所重み付け線形回帰を実行
A <- cbind(y_embed[valid_id, 2:3], 1) * weights

# 特異値分解を実行
A_svd <- svd(A) # 特異値分解
singular_value <- A_svd$d # 特異値を取り出す
singular_inv <- diag(1/singular_value, 3, 3) # 行列に変換
map <- A_svd$v %*% singular_inv %*% t(A_svd$u) %*% (weights * y_embed[valid_id + 1, 2]) # 変換するための写像を計算

# 予測値を計算
pred <- sum(map * c(as.numeric(y_embed[t_target, 2:3]), 1))

特異値分解の部分は馴染みがないととっつきにくいかもしれませんが、要するに (重み付けした) 線形モデルを線形代数のテクニックを使って解いている、という部分です。通常の lm() 関数でも weights オプションを指定すれば解くことができますが、解が不安定になってしまうこともあります。rEDM で特異値分解を利用する場合でも小さすぎる特異値が出現する場合はそれらは除かれています (ここの 656 行目辺り)。

S-map による予測は式で書くと以下のようになります。

$$\hat{Y}(t_{target}+1) = C_0 + \sum_{i=0}^{E-1} C_i Y(t_{target}-i)$$

2 次元で埋め込んでいる場合は、$\sum$ を用いずに書くと以下のとおりです。

$$\hat{Y}(t_{target}+1) = C_0 + C_1 \times Y(t_{target}) + C_2 \times Y(t_{target}-1)$$

この式の中の $C_0$、$C_1$、$C_2$ を上の局所重み付け線形回帰で決定しています。


さて、上記のようにして求めた予測値は以下のとおりです。実測値に近い値を予測できていることが分かります。

> pred # 予測値
[1] 0.4364812

> y_embed[t_target+1, 2] # 実測値
[1] 0.4675357

rEDM::s_map() による実装

さて、上記の計算をすべてのターゲット点に対して繰り返すことで、S-map の計算は完了します。rEDM を用いてそれを行うと、以下のようなコードになります。

# rEDM による実装
smap_res <- s_map(d$y, E = 2, theta = 1, stats_only = F, silent = T)

上記のコードでは埋め込み次元 $E=2$、$\theta=1$ として、予測精度などの統計情報以外に実際の予測値も出力しています。$t = 240$ のときの予測値を出力してみると rEDM を利用せずに計算したときと同じ結果になっていることが分かります。

> smap_res$model_output[[1]][240,]
    time       obs      pred   pred_var
240  241 0.4675357 0.4364812 0.02623754

rEDM::s_map() 関数のオプションの意味は以下のとおりです。

s_map(time_series,                             # 時系列データ
      lib = c(1, NROW(time_series)),           # 埋め込みに使用するデータの範囲
      pred = lib,                              # 予測値を出すデータの範囲
      norm = 2,                                # 重み付けのための距離の計算方法
      E = 1:10,                                # 試行する埋め込み次元の範囲
      tau = 1,                                 # 時間遅れの単位
      tp = 1,                                  # 何ステップ先を予測するか
      num_neighbors = 0,                       # 使用する最近傍点の数
      theta = c(0, 1e-04, 3e-04, 0.001, 0.003, # 試行する theta の値
                0.01, 0.03, 0.1, 0.3, 0.5,
                0.75, 1, 1.5, 2, 3, 4, 6, 8),
      stats_only = TRUE,                       # 予測の統計情報のみを出力するのか、それとも予測値も実際に出力するのか
      exclusion_radius = NULL,                 # 最近傍点を探す条件として、予測したい点とどの程度時間的に近い点を除くか (数値を指定する)
      epsilon = NULL,                          # 最近傍点を探す条件として、予測したい点とどの程度距離的に遠い点を除くか (数値を指定する)
      silent = FALSE,                          # Warning message を出力するかどうか
      save_smap_coefficients = FALSE)          # S-map の局所重み付け線形回帰の係数を出力に含めるかどうか

多くのオプションは simplex() 関数と同じですが、いくつか異なるものがあります。

• theta のオプションですが、これは上で説明したとおり、挙動を予測する上で近傍点の挙動をどのくらい重要とするかを決定する S-map を特徴づける変数であり、非常に重要です。
• E も重要ですが、現在のところは Simplex projection で決定した値を利用する方法が標準的です。
• num_neighbors のオプションは Simplex projection と同じく予測に使用する近傍点の数を指定します。デフォルトでは全ての点を利用するので、num_neighbors = NROW(time_series) としても良さそうですが、ここに 1 以下の数字を指定すると自動的に全ての点を利用して S-map を実行してくれます。また、例えば s_map() 関数を利用しつつ、num_neighbors = "E+1" と指定すると、E+1 個の近傍点を利用して局所重み付け線形回帰を行う、という S-map と Simplex projection の中間のようなアルゴリズムとなります。あまり利用されている場面は見ませんが。
• save_smap_coefficients のオプションは局所重み付け線形回帰のときの係数を出力に含めるかどうかを指定します。このオプションが特に用意されている理由は、この係数に特別な意味をもたせることが出来るからなのですが、これについては後述します。予測結果のみに興味があるときはこのオプションは FALSE で問題ありません。

3. 非線形性の定量化

S-map を行う上で $\theta$ の値は非常に重要です。$\theta > 0$ のときは程度の差はあれ、よりターゲットに近い点の挙動を重要視して未来を予測します。状態空間上で互いの距離が近いということはその点が置かれている状況が過去の状況も含めて似ているということです。また、より近い点の挙動を参照したときに未来の予測精度が上がるということは、その動態はより状況依存的であると考えることができます。つまり、2次元で埋め込んだときは、$\{Y(t), Y(t-1)\}$ が 2 次元空間にプロットされていますが、動態が状況依存的であるときは $Y(t)$ だけでなく$Y(t-1)$ も同じような値でないとうまく未来は予測できません。一方で、動態が状況依存的でなくどのような状況であろうとも常に同じようなルールで動くときは (=線形なルール)、$Y(t)$ や $Y(t-1)$ の値に関係なく同じモデルで未来を予測できます。このような状況のときは $\theta = 0$ で動態の予測精度が良くなります。ちなみに $\theta = 0$ のとき、S-map はいわゆる VAR モデル (Vector Auto Regressive model) と同じになります。

$\theta$ の値は予測における状況依存性 (非線形性) の強さを制御できるため、どの $\theta$ で予測精度 (相関係数、MAE、RMSEなど) が最も良くなるかを調べることで、動態の非線形性の定量が可能です。予測精度が最も良くなる $\theta$ が大きければ大きいほど、その動態の非線形性は大きい、というわけです。

rEDM の s_map() 関数では以下のようにして簡単に非線形性の定量ができます。

# theta の値を変えて予測精度をチェック
smap_res2 <- s_map(d$y, E = 2, silent = T)

ここでは theta の値を特に指定していません。デフォルトで theta = c(0, 1e-04, 3e-04, 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 0.5, 0.75, 1, 1.5, 2, 3, 4, 6, 8) となっているので、何も指定しなければ自動的にこれらの theta を試します。もちろん、自分で theta の値を指定することもできます。RMSE と $\theta$ の関係を図示すると以下のようになります。

plot(smap_res2$theta, smap_res2$rmse, type = "b", xlab = "theta", ylab = "RMSE", las = 1)

$\theta=8$ でもっとも予測精度が良くなっています。この場合、$Y(t)$ の動態は状況依存的であると判断できます。Simplex projection のときの $E$ と同じように調べる $\theta$ の値はどのくらいの範囲か、を決める明快な基準はありません。Hsieh et al. (2005) では $\theta = 1 \ or \ 2$ までチェックされています。Ushio et al. (2018) では $\theta = 10$ まで調べました。当然ですが、調べる $\theta$ の範囲が広がれば広がるほど計算時間が増加するので、計算時間の点から調べる $\theta$ の範囲が決まることもあると思います。s_map() 関数のデフォルトでは上限が $\theta = 8$ となっています。明快な根拠はないのですが、現時点ではその程度まで調べておけば十分なのではないかと思います。


次回は S-map の続きとして、｢Simplex projection との関係｣ ｢相互作用強度の推定｣ について書きます。また、EDM の最近の手法 についても書きたいと思います。

参考文献

• Hsieh CH, Glaser SM, Lucas AJ, Sugihara G. (2005) Distinguishing random environmental fluctuations from ecological catastrophes for the North Pacific Ocean. Nature 435, 336–340
• Sugihara G. (1994) Nonlinear forecasting for the classification of natural time series. Philos. Trans. R. Soc. A 348, 477–495
• Ushio M et al. (2018). Fluctuating interaction network and time-varying stability of a natural fish community. Nature 554, 360–363.
• Ye H, et al. (2015) Equation-free mechanistic ecosystem forecasting using empirical dynamic modeling. PNAS 112 (13) E1569-E1576
• Ye H, et al. (2018) rEDM: Applications of Empirical Dynamic Modeling from Time Series. doi:10.5281/zenodo.1935847

Written with StackEdit.