rEDM を用いた Empirical Dynamic Modeling: (4) S-map の周辺

S-map の周辺

前回は S-map のアルゴリズムを解説しました。S-map は様々な拡張･応用が可能で、最近でも新しい方法が提案されています。このポストでは、S-map の周辺の話題について書いていきたいと思います。

1. Simplex projection との関係

これまでに解説してきたように Simplex projection と S-map のアルゴリズムには共通する点とそうでない点があります。２つの手法では ｢時間遅れ座標系による埋め込み｣｢状態空間上でのターゲットからその他の点への距離の計算｣ は共通でした。その後は ｢数個の最近傍点の挙動を参照するのか｣ ｢全ての点に重みを付けてそれらの挙動を参照するのか｣ という違いと、予測値の計算が ｢重み付け平均か｣ ｢局所重み付け線形回帰か｣ という違いがありました。どちらも近未来の非線形予測に使われるため、使い分けが難しい部分があるかもしれません。そこで、両者の違いについて、経験的なものも含まれますがここで挙げておきます (引用をつけられないものもあるので、自分で試してみることをおすすめします)。将来的に項目を追加するかも知れません。

• EDM の枠組みでは、Simplex projection は最適埋め込み次元 $E$ の決定に用いられ、S-map は (決定した $E$ を用いて) 非線形性 $\theta$ を定量するのに用いられるのが標準的です。もちろん初めから S-map のみを用いて $E$ と $\theta$ を同時に決定する、という手順も考えられますが、この場合は試行する組み合わせが $(Eの数) \times (\thetaの数)$ となるため計算量が増大します。
• 時系列データの特性にもよりますが、S-map の方が Simplex projection より予測精度が高くなることが多いようです。
• 自分で試行して分かったこととして、Simplex projection は予測精度がしばしば S-map に劣るものの、誤差の入り方が比較的均一です。一方で S-map にはシステマティックな誤差が入ります (= 予測値が観測値に比べて小さくなりがちな領域や大きくなりがちな領域が存在する)。

2. 相互作用強度の推定

S-map の重要な利用方法として、多変数 S-map への拡張と S-map 係数 の利用があります。ここでも再び自作の２種系のデータセットを使用して説明します。

# 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,]

多変数 S-map

これまで Takens の埋め込み定理に基づき、一変数の時間遅れ座標系を利用して動態を再構成してきました。実はこの埋め込み定理は二変数以上の場合に一般化が可能です (Deyle & Sugihara 2011)。

上の時系列を例にすると動態 $Y(t)$ は $\{Y(t), Y(t-1)\}$ で再構成できますが、 $\{X(t), Y(t)\}$ でも再構成できる、ということです。この多変数埋め込みを利用した S-map を多変数 S-map (Multivariate S-map) と呼び、このときの局所線形回帰の係数に特別な意味をもたせることができます (Deyle et al. 2016)。どういうことか見てみましょう。

上の時系列を一変数での埋め込みを用いて S-map を実行したときの局所重み付け線形回帰の式は以下の通りでした。

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

そして、多変数 S-map を用いると次のようになります。
$$\hat{Y}(t_{target}+1) = C_0 + C_1 \times Y(t_{target}) + C_2 \times X(t_{target})$$

$C_i$ は S-map 係数 (S-map coefficients) と呼ばれますが、$C_i$ が何を表しているか考えてみましょう。例えば、$C_2$ は $X(t_{target})$ が少し増加したときに $Y(t_{target}+1)$ に及ぼす影響、と考えることができます。別の表し方をすれば、$C_2$ は $\partial Y(t_{target}+1)/{\partial X(t)}$ と書けるでしょう。つまり、S-map 係数 $C_2$ は $X(t)$ が $Y(t+1)$ に及ぼす影響 と解釈できます。$X$ や $Y$ がどちらも生物であれば S-map 係数が 種間相互作用強度 を表していると解釈できますし、$X$ が環境変数、$Y$ が生物であれば、環境が生物に与える影響を表していると解釈することが可能です。さらに S-map では $C_i$ を各時間ごとに計算するため、時間に伴って変数間の相互作用強度が変わる場合にもそれを柔軟に追跡することが可能です (Deyle et al. 2016)。

rEDM による多変数 S-map と S-map 係数の出力

では、rEDM を用いて多変数 S-map を実行してみましょう。まず、埋め込みに使用する変数を指定して、その後 block_lnlp() という関数を実行します (多変数 S-map は s_map() 関数では実行できません)。

# 多変数埋め込みを作成
smap_block <- cbind(d$x, d$y)

# 多変数 S-map を実行し最適な theta を決定
m_smap_res <- block_lnlp(smap_block, method = "s-map",
                         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),
                         silent = TRUE)

上では、通常の S-map と同様にいくつかの $\theta$ を試して最も最適な $\theta$ を探しています。$\theta$ と RMSE の関係は以下のとおりです。

plot(m_smap_res$theta, m_smap_res$rmse, type = "b")

$\theta=8$ が最適なようなので、それを使って多変数 S-map を実行し、予測結果や S-map 係数を出力します。

# 多変数 S-map を実行
m_smap_res2 <- block_lnlp(smap_block, method = "s-map", theta = 8, silent = TRUE, save_smap_coefficients = TRUE)

m_smap_res2 の中身は以下のとおりです (出力の後半は省いています)。

> m_smap_res2
  embedding tp nn theta num_pred       rho         mae        rmse
1      1, 2  1  0     8      900 0.9999472 0.003973874 0.004498086

ちなみに一変数で埋め込んだときの結果は以下の通りです。

> s_map(d$y, E = 2, theta = 8, silent = TRUE)
  E tau tp nn theta num_pred       rho         mae        rmse
1 2   1  1  0     8      899 0.9989587 0.006919698 0.008887697

比べると多変数 S-map では予測精度が向上していることが分かります (一変数での RMSE = 0.0088…, 二変数での RMSE = 0.0044…)。埋め込みに使用する変数 (この例では $X$) が予測する変数 ($Y$) に因果があるとき、予測精度は多変数 S-map によって向上する傾向にあります。

S-map 係数は m_smap_res2$smap_coefficients 内に格納されています。

> head(m_smap_res2$smap_coefficients[[1]])
         c_1          c_2       c_0
1  0.5747497  0.024614030 0.6633858
2 -3.0510028 -0.009597196 3.0771957
3  1.7387918 -0.029947627 0.2782008
4 -1.7066457 -0.088475705 2.0510882
5 -1.9924075  0.042627392 2.1762527
6 -1.4671749  0.005263682 1.8050193

また、S-map 係数の分布を箱ひげ図で図示すると以下のとおりです。

boxplot(m_smap_res2$smap_coefficients[[1]][,1:2], las = 1, ylab = "S-map coefficients")
abline(h = 0, lty = 2)

$C_1$ は $Y(t)$ から $Y(t+1)$ への効果、$C_2$ は $X(t)$ から $Y(t+1)$ への効果です。$C_2$ は少しわかりにくいですが、平均的には負のようです。元の方程式から考えると、$C_1$ も $C_2$ も負で $C_1 \ll C_2$ ですから、S-map 係数の結果と合っているようです (ただし、元の方程式の相互作用強度は個体あたりの相互作用強度で、S-map 係数は個体群あたりの影響なので、単純に両者が対応するわけではないので注意が必要です)。

rEDM::block_lnlp() のオプション

ここで、block_lnlp() のオプションについて解説しておきます。

block_lnlp(block,                             # 時系列データ
           lib = c(1, NROW(time_series)),     # 埋め込みに使用するデータの範囲
           pred = lib,                        # 予測値を出すデータの範囲
           norm = 2,                          # 重み付けのための距離の計算方法
           method = c("simplex", "s-map"),    # Simplex projection を利用するか、S-map を利用するか
           tp = 1,                            # 何ステップ先を予測するか
           num_neighbors = switch(match.arg(method), simplex = "e+1", `s-map` = 0),
                                              # 予測に利用する近傍点の個数 (デフォルトは方法によって自動設定)
           columns = NULL,                    # 状態空間の再構成に使用する列
           target_column = 1,                 # 予測のターゲットになる変数の列
           stats_only = TRUE,                 # 予測の統計情報のみを出力するのか、それとも予測値も実際に出力するのか
           first_column_time = FALSE,         # block の第一列目が時間かどうか
           exclusion_radius = NULL,           # 最近傍点を探す条件として、予測したい点とどの程度時間的に近い点を除くか (数値を指定する)
           epsilon = NULL,                    # 最近傍点を探す条件として、予測したい点とどの程度距離的に遠い点を除くか (数値を指定する)
           theta = NULL,                      # 試行する theta の値
           silent = FALSE,                    # Warning message を出力するかどうか
           save_smap_coefficients = FALSE)    # S-map の局所重み付け線形回帰の係数を出力に含めるかどうか

ほとんど simplex() や s_map() と同じですが、いくつか違う点があります。まず、埋め込み次元を指定する E がありません。block_lnlp() では block をこれまでの一変数のベクトルではなく、データフレームもしくは行列として与えることで E を同時に指定していることになるためです。つまり、block が２列のデータフレームであれば自動的に E = 2 と指定したことになります (columns のオプションを指定することで与えたデータフレームのどの列を埋め込みに用いるかの指定も可能です)。また、method で Simplex projection を実行するか、S-map を実行するか指定できます。多変数 S-map を実行したいときは method = "s-map" と指定する必要があります。

EDM における最近の発展手法

• S-map は局所的に見ると単なる (重み付け) 重回帰であり、そのため様々な統計手法を取り入れることで容易に拡張できます。例えば、LASSO や Ridge など正則化を入れると、正則化 S-map (Regularized S-map; Cenci et al 2019) となり、こちらでいくつか試したところ予測精度の向上が見られました。Cenci et al. (2019) に R の解析コードが付属していますので、興味ある方は試してみて下さい。また、自分でも異なる実装をして Github で公開していますので、良ければこちらもお試し下さい (https://github.com/ong8181/random-scripts)。

• EDM に明示的に階層性や空間構造を取り込むために ガウス過程 (Gausiaan Process) を組み込んだ手法も提案されています。Munch et al. (2017) においてそれらが解説されていますし、rEDM にもすでに block_gp() という関数が実装されておりガウス過程 EDM が実行できます (空間構造を明示的に指定はできないようですが)。少し試してみたのですが予測精度がかなり良く、内挿に強い、という印象です。

• 多変数埋め込みをさらに拡張し予測精度を向上させた Multi-view embedding (多視点埋め込み?とでも訳したらよいでしょうか…) という手法も考案されています (Ye & Sugihara 2016)。様々な埋め込みにより少しずつ違う視点で状態空間を再構成し、予測精度の良い埋め込みを選抜した後、それらの予測を統合する、という方法により、多変数 S-map 法よりも良い予測を行います。こちらも rEDM に multiview() 関数として実装されています。アイデアがすごく好きなのですが、計算時間がかかることが多く、今の所自分はあまり利用していません。

• Multi-view embedding と似ていますが、更に大規模かつランダムにたくさんの埋め込みを作り、それらの予測を統合することで少ない時間点･多変数の時系列データを用いて正確かつ長期の予測が可能になる、という Randomly Distributed Embedding (RDE) という方法も提案されています (Ma et al. 2018) (まだ自分で試したことはありません)。

参考文献

• Cenci S, Sugihara G, Saavedra S (2019) Regularized S-map for inference and forecasting with noisy ecological time series. Methods in Ecology and Evolution 10, 650–660
• Deyle ER, May RM, Munch SB, Sugihara G. (2016) Tracking and forecasting ecosystem interactions in real time. Proc. R. Soc. Lond. B 283, 20152258
• Deyle ER, Sugihara G (2011) Generalized theorems for nonlinear state space reconstruction. PLoS ONE 6, e18295.
• Ma H, Leng S, Aihara K, Lin W, Chen L (2018) Randomly distributed embedding making short-term high-dimensional data predictable. PNAS, 115 E9994–10002
• Munch SB, Poynor V, Arriaza JL (2017) Circumventing structural uncertainty: A Bayesian perspective on nonlinear forecasting for ecology. Ecological Complexity 32, 134–143
• Sugihara G. (1994) Nonlinear forecasting for the classification of natural time series. Philos. Trans. R. Soc. A 348, 477–495
• Takens, F. (1981) Detecting strange attractors in turbulence. In D. Rand & L.-S. Young (Eds.), Lecture Notes in Mathematics (Vol. 898, pp. 366–381).
• Ye H, Sugihara G. (2016). Information leverage in interconnected ecosystems: Overcoming the curse of dimensionality. Science, 353, 922–925
• Ye H, et al. (2018) rEDM: Applications of Empirical Dynamic Modeling from Time Series. doi:10.5281/zenodo.1935847

Written with StackEdit.

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.