2019年12月11日水曜日

rEDM を用いた Empirical Dynamic Modeling: (2) 近未来予測 (Simplex projection)

20191211_Blogger0005

再構成した動態の利用方法

前回のポストでは一つの変数から興味ある系の動態を再構成する方法について紹介しました。この操作は時間遅れ時系列による状態空間の再構成 (State Space Reconstruction; SSR) とも呼ばれます。SSR により得られた｢形 (= 動態)｣ は、測定していない変数の情報や系の駆動メカニズムに関する情報など、いくつもの有用な情報を含んでいます。Empirical Dynamic Modeling (EDM) と呼ばれる非線形時系列解析の枠組みを利用することで、これらの情報を取り出し、様々な解析をすることが可能になります。

現時点での EDM の主な応用としては、近未来予測 (Sugihara & May 1990; Sugihara 1994; Ye & Sugihara 2016)･動態特性の定量 (Hsieh et al. 2005)･変数間の因果関係の検出 (Sugihara et al. 2012)･変数間の相互作用強度の定量 (Deyle et al. 2016)･動態の安定性の定量 (Ushio et al. 2018)、などが挙げられます。

ここからモデル時系列を用いて EDM による解析を具体的に説明していきます。ここでは、R 言語を使用し、rEDM という EDM のためのパッケージを用いることにします (Ye et al. 2015; Ye et al. 2018)。Python で EDM を行うための pyEDM というパッケージも存在しますが、それについてはまたの機会に解説したいと思います (こちらのサイトで解析が試されています)。

また、具体的な解説を始める前に現時点で EDM についてのまとまった情報が得られる論文やサイトなどを以下に紹介しておきます。

また、EDM の理論的背景となる多様体やカオス時系列解析に関しては自分は以下の資料を参考にして勉強しました (というか未だに勉強中という感じですが)。

近未来予測 (Simplex projection)

ここからは R と rEDM を使い、解析を行いつつ EDM の手法について解説していきます。まずは、近未来予測についてです。

参考: Simplex projection は異なるいくつかの説明の仕方が可能です。このポストでの説明は比較的一般的なものと思いますが、状態空間上ではなく元の時系列 (横軸が時間･縦軸が解析対象の変数) で考えたほうが分かりやすいこともあります。そのような説明も読みたい方は、以下のポストを読んだあとにぜひ Simplex projection walkthrough も御覧ください (英語ですが)。

1. rEDM とデモ時系列データの準備

始めに rEDM のインストールとデモ用の時系列データの準備です。rEDM パッケージ内にもでもデータが用意されていますが、今回はモデルを明示的にするために自作しています。
(こちらの環境は R のバージョンは 3.6.1 で、OS は macOS Mojave です)

注意点: v0.7.3以前とそれ以降では確か計算結果の出力フォーマットが異なるので注意が必要です。v0.7.5 を使用して下さい。

# ライブラリのインストール
install.packages("remotes")
library(remotes)
remotes::install_github("ha0ye/rEDM")

# 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){ # Sugihara et al. (2012) Fig.1 のモデル
        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,]

# 一部を表示してみる
plot(d$time, d$x, type = "l", xlim = c(100,150),
     las = 1, xlab = "Time", ylab = "Value", col = "royalblue")
lines(d$time, d$y, col = "red3") # Xが青・Yが赤です

enter image description here

ここで使用するモデルデータはロトカ･ボルテラの方程式に基づいて生成されたカオス動態で、XX と YY が相互作用しています。

注意点: 通常、EDM では時系列は全て標準化 (平均 0, 分散 1 に変換) してから解析を行います。EDM では解析の途中で点間のユークリッド距離を求めたりしますが、標準化しないと、距離が変数の単位 (g, kg など) や平均値に大きく依存してしまいます。それを防ぐためです。ただ、今回は誤差なしのモデル時系列であることと、図示したときの直感的な理解しやすさを優先して、標準化せずに解析を進めます。自分のデータで解析するときは scale() などの関数で標準化をしてから解析を進めて下さい。今回の場合はd$x <- as.numeric(scale(d$x)); d$y <- as.numeric(scale(d$y)) などとすれば良いでしょう。

2. Simplex projection による近未来予測と最適埋め込み次元の推定

Simplex projection とは

Simplex projection は注目するデータ点の周りにある点 (近傍点) の挙動を元に、注目する点の近未来を予測する手法です (Sugihara & May 1990)。上のモデルデータを利用してそのアルゴリズムを解説します。

時間遅れ時系列の作成

ここでは YY に注目して解析を進めてみることにします。まずは時間遅れ埋め込みにより YY の動態を再構成します。

# Yの時間遅れ時系列を作る
y_embed <- make_block(d$y, max_lag = 2) # 暫定的に Y(t), Y(t-1)まで作成

make_block() 関数を使うと、y_embed の中身は以下のようになります。ラグを取るだけなので、マニュアルで作成することも容易ですが、ここでは rEDM の関数を使うことにします。head(y_embed)の出力は以下のとおりです。

> head(y_embed) # 中身を確認
  time      col1    col1_1
1    1 0.4485858        NA
2    2 0.8468057 0.4485858
3    3 0.3758589 0.8468057
4    4 0.8115287 0.3758589
5    5 0.4771027 0.8115287
6    6 0.8369673 0.4771027

次に、2 次元で埋め込んだ場合、どのような動態を描けるか図示してみます。

plot(y_embed$col1_1, y_embed$col1, col = rgb(0,0,0,0.5),
     xlab = "Y(t-1)", ylab = "Y(t)", las = 1)

何らかの ｢軌道｣ のようなものが見えていますね。この軌道に対して Simplex projection を適用して、1ステップ未来の予測を試みます。

最近傍点の特定

Simplex projection では、まず近未来を予測したい点を再構成した動態の中で一つ設定します。例えば、t=240t = 240 の点の一時間点先の未来に興味があるとし、ttarget=240t_{target} = 240 とします。

t_target <- 240

YY を埋め込んだものを Y(t)=[Y(t),Y(t−1)]\boldsymbol{Y}(t) =[Y(t), Y(t-1)] とすると、次の点 Y(ttaregt)\boldsymbol{Y}(t_{taregt}) がひとまずのターゲットとなります。

> y_embed[t_target,]
    time      col1    col1_1
240  240 0.8293677 0.4683887

再構成した動態の中では ttarget=240t_{target} = 240 は以下の赤色の点です。

plot(y_embed$col1_1, y_embed$col1, col = rgb(0,0,0,0.5),
     xlab = "Y(t-1)", ylab = "Y(t)", las = 1)
points(y_embed[t_target,3:2], col = "black", bg = "red3", pch = 21, cex = 2)

もう少し拡大すると以下のようになります。

plot(y_embed$col1_1, y_embed$col1, col = rgb(0,0,0,0.5),
     xlab = "Y(t-1)", ylab = "Y(t)", las = 1,
     xlim = c(y_embed[t_target,3] - 0.005, y_embed[t_target,3] + 0.005),
     ylim = c(y_embed[t_target,2] - 0.005, y_embed[t_target,2] + 0.005))
points(y_embed[t_target,3:2], col = "black", bg = "red3", pch = 21, cex = 0.9)

最近傍点はターゲットの点からのユークリッド距離を計算して求めます。Simplex projection ではデフォルトでは (埋め込み次元 +1+ 1) 個の最近傍点を選びます。今は E=2E=2 で埋め込んでいるので、3 個の最近傍点を探すことになります (この選ぶ個数はオプションで変更可能)。

# t = 240 の点から各点への距離を計算
distance <- sqrt(rowSums((y_embed[,2:3] - c(y_embed[t_target, 2:3]))^2))
nn_index <- order(distance)[2:4] # 自分自身を除いて 3 点近い点を選ぶ

> nn_index
[1] 628 564 684

628,564,684628, 564, 684 の 3 点が最近傍点のようです。それぞれを t1=628,t2=564,t3=684t_1 = 628, t_2 = 564, t_3 = 684 としましょう。Y(t1),Y(t2),Y(t3)\boldsymbol{Y}(t_1), \boldsymbol{Y}(t_2), \boldsymbol{Y}(t_3) の位置を青色で点で確認してみます。

plot(y_embed$col1_1, y_embed$col1, col = rgb(0,0,0,0.5),
     xlab = "Y(t-1)", ylab = "Y(t)", las = 1,
     xlim = c(y_embed[t_target,3] - 0.005, y_embed[t_target,3] + 0.005),
     ylim = c(y_embed[t_target,2] - 0.005, y_embed[t_target,2] + 0.005))
points(y_embed[t_target,3:2], col = "black", bg = "red3", pch = 21, cex = 0.9)
points(y_embed[nn_index,3:2], col = "black", bg = "royalblue", pch = 21, cex = 0.9)

しっかりと最近傍点が特定できているようです。

近未来予測

さて、次に Y(t1),Y(t2),Y(t3)\boldsymbol{Y}(t_1), \boldsymbol{Y}(t_2), \boldsymbol{Y}(t_3) の挙動を元に Y(ttarget)\boldsymbol{Y}(t_{target}) の近未来予測を行います。まずは青色の点の挙動を調べましょう。

> # 最近傍点の挙動を調べる
> # 最近傍点の現在地
> y_embed[nn_index,]
    time      col1    col1_1
628  628 0.8294636 0.4680146
564  564 0.8288434 0.4697511
684  684 0.8310493 0.4683390

> # 最近傍点の１時点未来
> y_embed[nn_index + 1,]
    time      col1    col1_1
629  629 0.4665281 0.8294636
565  565 0.4725765 0.8288434
685  685 0.4548780 0.8310493

y_embed[nn_index + 1,]が最近傍点の一つ先の未来の点です。これらについても図示してみましょう。プロットの右下の方に移動しているようです。

plot(y_embed$col1_1, y_embed$col1, col = rgb(0,0,0,0.5),
     xlab = "Y(t-1)", ylab = "Y(t)", las = 1)
points(y_embed[t_target,3:2], col = "black", bg = "red3", pch = 21, cex = 1.2)
points(y_embed[nn_index+1,3:2], col = "black", bg = "green3", pch = 24, cex = 1.2)

Y(t1),Y(t2),Y(t3)\boldsymbol{Y}(t_1), \boldsymbol{Y}(t_2), \boldsymbol{Y}(t_3) が 次の時間には右下 (水色の三角の点; Y(t1+1),Y(t2+1),Y(t3+1)\boldsymbol{Y}(t_1+1), \boldsymbol{Y}(t_2+1), \boldsymbol{Y}(t_3+1)) に移動しています。ここで (うまく動態が再構成できていて、しっかり軌道が捉えられていれば) Y(ttarget)\boldsymbol{Y}(t_{target}) の挙動は再構成した動態上での Y(t1),Y(t2),Y(t3)\boldsymbol{Y}(t_1), \boldsymbol{Y}(t_2), \boldsymbol{Y}(t_3) の挙動と似ているはず、ということを利用して Y(ttarget+1)\boldsymbol{Y}(t_{target}+1) を予測します。実際の計算では最近傍点のターゲットの点までの距離に応じた重み付けをすることで より近い点の挙動を重視して Y(ttarget+1)\boldsymbol{Y}(t_{target}+1) を予測します。

# 最短距離を取得
min_distance <- distance[nn_index][1]
# 最短距離で補正した重みを計算
weights <- exp(-distance[nn_index]/min_distance)
# 合計の重みを計算
total_weight <- sum(weights)
# 重み付け平均を計算
pred <- (weights %*% as.matrix(y_embed[nn_index + 1,3:2])) / total_weight

上のRコードは式で表すと、以下のようにして Y(ttarget+1)\boldsymbol{Y}(t_{target}+1) を予測しています。Y^(ttarget+1)\boldsymbol{\hat{Y}}(t_{target}+1) は Y(ttarget+1)\boldsymbol{Y}(t_{target}+1) の予測値、wtotalw_{total} は total_weight、w1w_1, w2w_2, w3w_3 は weigths の中に格納されているそれぞれの重みを表しています。

Y^(ttarget+1)={w1×Y(t1+1)+w2×Y(t2+1)+w3×Y(t3+1)}wtotal\boldsymbol{\hat{Y}}(t_{target}+1) = \frac{ \{ w_1 \times \boldsymbol{Y}(t_1+1) + w_2 \times \boldsymbol{Y}(t_2+1) + w_3 \times \boldsymbol{Y}(t_3+1)\} }{w_{total}}

予測の結果は以下の通りで非常に近い値を予測できていることがわかります。

> pred[2] # 予測値
[1] 0.4664999

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

予測を図示すると以下のような感じです。水色の三角が Y(t1+1),Y(t2+1),Y(t3+1)\boldsymbol{Y}(t_1+1), \boldsymbol{Y}(t_2+1), \boldsymbol{Y}(t_3+1)で、バツ印がそれらの重み付け平均で予測した値です。きちんと Y(ttarget+1)\boldsymbol{Y}(t_{target}+1) (赤三角の点) の実測値の近くを予測できていることが分かります。

plot(y_embed$col1_1, y_embed$col1, col = rgb(0,0,0,0.5),
     xlab = "Y(t-1)", ylab = "Y(t)", las = 1,
     xlim = c(y_embed[t_target+1,3] - 0.001, y_embed[t_target+1,3] + 0.002),
     ylim = c(y_embed[t_target+1,2] - 0.02, y_embed[t_target+1,2] + 0.01))
points(y_embed[nn_index+1,3:2], col = "black", bg = "green3", pch = 24, cex = 1.2)
points(y_embed[t_target+1,3:2], col = "black", bg = "red3", pch = 24, cex = 1.2)
points(pred, col = "red3", pch = 4, cex = 2)

以上の計算を全ての点について繰り返すことで Simplex projection による予測精度が計算できます。

rEDM::simplex() による実装

以上、一つ一つ手順を確認しながら計算していきましたが、rEDM::simplex()関数を利用することで簡単に計算を実行できます。simplex()関数にはいくつかのオプションが指定できますが、簡単には以下のコードで Simplex projection が実行できます。

# Simplex 関数による計算
simp_res <- simplex(d$y, E = 2, stats_only = F)

t=240t = 240 のときの予測値は以下の通りであり、マニュアルで計算したときと同じことが分かります。stats_only = F とすることで統計情報だけではなく実際の予測値も出力してくれます。

> simp_res$model_output[[1]][240,]
    time       obs      pred     pred_var
240  241 0.4675357 0.4664999 6.384118e-06

各オプションとそれらの意味は以下のとおりです。

simplex(time_series,                     # 解析する時系列 (1変数)
        lib = c(1, NROW(time_series)),   # 埋め込みに使用するデータの範囲
        pred = lib,                      # 予測値を出すデータの範囲
        norm = 2,                        # 重み付けのための距離の計算方法
        E = 1:10,                        # 試行する埋め込み次元の範囲
        tau = 1,                         # 時間遅れの単位
        tp = 1,                          # 何ステップ先を予測するか
        num_neighbors = "e+1",           # 使用する最近傍点の数
        stats_only = TRUE,               # 予測の統計情報のみを出力するのか、それとも予測値も実際に出力するのか
        exclusion_radius = NULL,         # 最近傍点を探す条件として、予測したい点とどの程度時間的に近い点を除くか (数値を指定する)
        epsilon = NULL,                  # 最近傍点を探す条件として、予測したい点とどの程度距離的に遠い点を除くか (数値を指定する)
        silent = FALSE)                  # Warning message を出力するかどうか

以下、いくつかオプションを変更して Simplex projection を実行してみます。最も重要なオプションは E ですが、それについては次のセクションで解説します。

  1. 統計情報のみを表示し、Warning を表示しない (いくつかの出力は長くなるので省いています)。stats_only = T がデフォルトなので関数内には何も記載していません。
> simplex(d$y, E = 2, silent = T)
  E tau tp nn num_pred       rho         mae        rmse perc p_val
1 2   1  1  3      899 0.9995513 0.002529848 0.005832568    1     0
  1. 予測する未来の時間点を tp = 1:10 と指定。つまり、1 時間点先から 10 時間点先までを予測。予測する未来が遠くなると予測精度が悪くなっていきます。つまり、予測値と実測値の相関係数 rho が下がっていったり、誤差の指標 mae (Mean Square Error) や rmse (Root Mean Square Error) が上がっていったりします。これは初期値鋭敏性をもつカオス動態の特徴の一つです。
> simplex(d$y, E = 2, tp = 1:10, stats_only = T, silent = T)
   E tau tp nn num_pred       rho         mae        rmse perc p_val
1  2   1  1  3      899 0.9995513 0.002529848 0.005832568    1     0
2  2   1  2  3      898 0.9990229 0.003728225 0.008609239    1     0
3  2   1  3  3      897 0.9974562 0.006427730 0.013887861    1     0
4  2   1  4  3      896 0.9959668 0.008843418 0.017530732    1     0
5  2   1  5  3      895 0.9936978 0.012128339 0.021840383    1     0
6  2   1  6  3      894 0.9906415 0.015322926 0.026606249    1     0
7  2   1  7  3      893 0.9887471 0.018156332 0.029168297    1     0
8  2   1  8  3      892 0.9860934 0.021092473 0.032412541    1     0
9  2   1  9  3      891 0.9834527 0.023851365 0.035337366    1     0
10 2   1 10  3      890 0.9816744 0.025640024 0.037179981    1     0
  1. データの半分をトレーニングデータとして動態の再構成に利用し、それを用いて残り半分を予測する。データが十分長い場合はこのような方法で予測性能を評価することがよくあります。データが短い場合は lib と pred を完全に一致させることで Leave-One-Out Cross Validation (LOOCV) と呼ばれる予測性能の評価を実行します (こちらがデフォルト設定)。Cross validation (交差検証) については Wikipedia にも載っているので参照して下さい。
> simplex(d$y, E = 2,
          lib = c(1, 450),     # 時間 t = 1-450 の点で動態を再構成
          pred = c(451, 901),  # 時間 t = 451-901 の点の1時間点先を予測
          tp = 1, stats_only = T, silent = T)
  E tau tp nn num_pred       rho         mae        rmse perc p_val
1 2   1  1  3      449 0.9991972 0.003711699 0.007815253    1     0

その他、tau もしばしば重要な変数となります。tau は埋め込みの際のタイムラグの値をどうするか、という変数です。つまり、tau = 1 のときは、埋め込みを {Y(t),Y(t−1),Y(t−2),...}\{Y(t), Y(t-1), Y(t-2), ...\} として作っていきますが、tau = 2 のときは埋め込みを {Y(t),Y(t−2),Y(t−4),...}\{Y(t), Y(t-2), Y(t-4), ...\} として作っていきます。実用上は tau = 1 が使用されている場合がほとんどのようですが、対象の時系列データの特性を考えて tau = 1 以外を使うべき場合も十分ありえます。適切な tau を決めるためには、いくつかの tau を試してみて予測性能を比較してみるというやり方もできますし、対象の時系列データの特性を考えて決めるというやり方もありえます。また、EDM の枠組み以外の方法で ```tau`` を決める方法も多数提案されています (例えば、自己相関関数や相互情報量を利用した方法)。

最適埋め込み次元の推定

さて、simplex() 関数を利用することで異なる E に対する予測精度を一度に求めることができます。Simplex projection で予測性能が良い、ということはどういう状況か少し考えてみましょう。

最近傍点の挙動に基づいた予測がうまくいっているということは、注目する点と最近傍点の挙動が似ており、状態空間上で近い点はほぼ同じように挙動している、ということです。従って、状態空間内で動態がなめらかに描けている、ということが言えるでしょう。逆に、予測精度が低いということは、状態空間上で複数の点がある時点で同じような場所にあったとしても次の時間点ではてんでばらばらの場所に移動しており、描いた軌道はぐちゃぐちゃだった、ということが言えそうです。

｢埋め込み｣ とは ｢(観測できないかもしれない) 真の動態｣ と ｢再構成した動態｣ が ｢1:1 かつなめらか｣ に対応している状態です。従って、時間遅れ埋め込みにより系の状態がきちんと再構成できたかどうかを判断するために Simplex projection による予測精度が指標として使えるかもしれません。なめらかできれいな軌道が描けている場合は、Simplex projection による予測精度が高いはずだからで、｢きちんと｣ 動態を描けていると言ってもいいでしょう (前回のポストで ｢きちんと｣ の定義は後で考えると書きましたが、それがここでの考察に該当します)。

では、E を様々に変えて予測精度の変化を見てみましょう。

simp_res2 <- simplex(d$y, E = 1:10, silent = T)

E に伴う予測精度の変化は以下のとおりです。今回は RMSE (Root Mean Square Error) を予測精度の指標として用いていて、低いほど予測精度が良いということを表しています。

plot(simp_res2$E, simp_res2$rmse, type = "b", xlab = "E", ylab = "RMSE", las = 1)

E=2E=2 で RMSE が最も低く、最も予測精度が良くなっています。従って、YY は E=2E=2 が最も埋め込み次元として適切なようです。このようにして推定された埋め込み次元を 最適埋め込み次元 (optimal embedding dimension) と言ったりします。

EDM の枠組みでは Simplex projection を用いて最適埋め込み次元を推定するのが一般的ですが、それ以外の方法も提案されています。他の方法としては、例えば、誤り近傍法 (False Nearest Neighbour) が有名だと思います。最適埋め込み次元の推定方法は議論が多いところで、注意深い試行錯誤が必要です。

注意点: 最適埋め込み次元の決定は EDM による解析のほぼ全てに影響を及ぼす重要なステップで、注意深く考えて行う必要があります。実際に解析を行う際に多くの人が悩む点がいくつかありますので、ここで列挙しておきます。

  • 予測精度を調べる EE はどこからどこまでか?

    • 下限は E=1E = 1 から調べれば良いと思いますが (ただし一番最後も読んで下さい)、上限は研究によってまちまちです。(おそらく何となく) E=10E = 10 となっている研究もありますが、これらについては ｢データの長さ｣ ｢対象の時系列の特性｣ ｢対象の時系列が含む誤差｣ などを元に決めると良いでしょう。
    • Munch et al. (2019) には、最大の埋め込み次元はせいぜいデータ長の平方根程度になりがち、との記述があります (Question 5 への答えの部分)。これに従うとデータ点 100 の時系列であれば E=10E = 10 かそれよりもう少し大きめまで調べれば OK ということになります (ただこれはどこまで一般化出来るのかは自分は検討していません)。
    • E=2E = 2 とすれば、過去の 2 時間点の情報が未来の挙動に影響していることになり、E=3E = 3 とすれば過去の 3 時間点までの情報が未来の挙動に影響していると考えることができます。自分が興味ある変数の挙動が過去どのくらいの時間までの影響を受けていそうか考えることが EE の上限を決める助けになるかもしれません。
    • 実データの場合は、多かれ少なかれ値に誤差を含んでいます。EE を大きくすればするほど状態空間内での点の位置にも誤差が積み重なることになるので、誤差の大きい時系列であるほど大きな EE を調べることが難しくなっていきます。
    • E=1E=1 を調べない、というやり方もあり得るようです。埋め込みは ｢なめらかに 1:1 対応｣ というのが定義でしたので、1 次元の時系列がなめらかになる、とは単調減少 or 単調増加の場合以外は想像しにくく、その意味で E=1E=1 で埋め込みになる状況はあり得ない、という考え方です。理論的には E=1E=1 が最適埋め込み次元というのはあり得るようですが。ただし、この点についてはまだコンセンサスはないようで、自分は念のため E=1E=1 も調べています。

  • どの予測精度の指標を利用すべきか?

    • rEDM が出力する予測精度は 3 種類で、相関係数･MAE･RMSE です。初期は相関係数が指標として利用されることも多かったですが (Sugihara et al. 2012 など)、相関係数は外れ値に影響されることも多いことから、MAE や RMSE の方がよい場合が多いように思います。自分は最近は RMSE に頼ることが多いです。

次回のポストでは同じく近未来予測の手法である S-map 法について解説します。

参考文献

  • Chang CW, Ushio M, Hsieh C. (2017) Empirical dynamic modeling for beginners. Ecological Research, 32, 785–796
  • Deyle ER, May RM, Munch SB, Sugihara G. (2016) Tracking and forecasting ecosystem interactions in real time. Proc. R. Soc. Lond. B 283, 20152258
  • 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
  • Munch, S. B., Brias, A., Sugihara, G., Rogers, T. L. (2019) Frequently asked questions about nonlinear dynamics and empirical dynamic modelling. ICES Journal of Marine Science, doi:10.1093/icesjms/fsz209
  • Sugihara, G. et al. (2012) Detecting causality in complex ecosystems. Science 338, 496–500
  • Sugihara G, May RM. (1990) Nonlinear forecasting as a way of distinguishing chaos from measurement error in time series. Nature 344, 734–741
  • 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, 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
  • 中山新一朗, 阿部真人, 岡村寛. (2015) Convergent cross mapping の紹介：生態学における時系列間の因果関係推定法. 日本生態学会誌, 65, 241-253

Written with StackEdit.

時刻:

0 件のコメント:

コメントを投稿

登録: コメントの投稿 (Atom)

Empirical Dynamic Modeling with rEDM: (2) Near-future forecasting (Simplex projection)

20231015_Blogger0021 *This is an English version of my previous post (first translated by ChatGPT, checked and edited ...