rEDM的用法（1）——Simplex projection and S-map analysis

Load package and time series

All the EDM analyses are carried out by the rEDM package. The red noise time series is generated by a difference equation specifying the value at next time step as proportional to the current value plus a random number drawn from the standard normal distribution. The logistic map time series is generated by a self-regulatory difference equation following Hsieh et al. (2005). Both models are simulated for 10000 time steps, but only the last 1000 data are used for further analyses in order to exclude the transient dynamics.

## loading R package: rEDM
library(rEDM)

Simplex projection (Sugihara & May 1990)

By simplex projection, we make one-step forward forecast (predict t+1 step) for the red noise or logistic map dynamics. Prior to simplex projection, the time series are normalized to zero mean and unit variance. The simplex projection is carried out by the function rEDM::simplex(). The object Red in the dataset is the normalized time series of red noise; the object Logi is the normalized time series of the logistic map. We divide the time series into two halves. The first half is used as the library set for manifold reconstruction. The second half is used as the target for out-of-sample prediction. The argument lib is the time index indicating the start (1) and the end (500) in the library set, respectively. Similarly, the argument pred indicates the time index in the prediction set. Then, we specify a sequence of testing embedding dimensions from 2 to 8 by the argument E=c(2:8)which executes the simplex projections using different embedding dimensions. Finally, we present the results showing the relationship between the predictive skill (ρ) and the embedding dimension (E) (Fig. 2c & d in the main text).

## Data loading 
dat <- read.csv('ESM2_Data_noise.csv',header=T)

## Data normalization 
Red <- ((dat[,"R"]-mean(dat[,"R"]))/sd(dat[,"R"]))
Logi <- ((dat[,"L"]-mean(dat[,"L"]))/sd(dat[,"L"]))
#######################################################################
## Simplex projection for red noise and logistic map
sim_r <- simplex(Red,lib=c(1,500),pred=c(501,1000),E=c(2:8))
sim_l <- simplex(Logi,lib=c(1,500),pred=c(501,1000),E=c(2:8))

## Plot predictive skill (rho) vs embedding dimension (E) 
par(mfrow=c(2,1),mar=c(4,4,1,1))
plot(rho~E,data=sim_r,type="l",xlab="Embedding dimension (E)",ylab=expression(rho),ylim=c(0.3,0.4),col=2,main="Red noise")
plot(rho~E,data=sim_l,type="l",xlab="Embedding dimension (E)",ylab=expression(rho),ylim=c(0.95,1.02),col=4,main="Logistic map")

## The optimal embedding dimension determined by maximizing rho
(E_r <-sim_r[which.max(sim_r$rho),"E"][1])# The optimal E of red noise
(E_l <-sim_l[which.max(sim_l$rho),"E"][1])# The optimal E of logistic map

 

最终输出E_r和E_l对应两者的最佳embedding dimension，分别为7和2，这与上图是吻合的、

S-map analysis

Using the function rEDM::s_map(), we calculate the one-step forward forecast by S-map analysis. Again, we divide the time series into two halves, one for library, lib=c(1,500), and another for out-of-sample prediction, pred=c(501,1000). The embedding dimensions, E=E_r for red noise and E=E_l for the logistic map, have been already determined by simplex projection (Fig. 2c & d) in the previous section. Then, we specify a sequence of testing state-dependency parameters θ from 0 to 2 with an increment of 0.1 by the argument, theta=seq(0,2,0.1)which executes S-map by trial-and-error using different θ. Finally, we demonstrate the results showing the relationship between the predictive skills (ρ) and state-dependency parameters (θ) (Fig. 2e & f). Again, the criterion maximizing predictive skill is applied to determine the optimal θ.

# S map for Red Noise & logistic map
smap_r <- s_map(Red,E=E_r,lib=c(1,500),pred=c(501,1000),theta=seq(0,2,0.1))
smap_l <- s_map(Logi,E=E_l,lib=c(1,500),pred=c(501,1000),theta=seq(0,2,0.1))

## Plot predictive skill (rho) vs state-dependency parameter (theta) 
plot(rho~theta,data=smap_r,type="l",xlab=expression(theta),ylab=expression(rho),ylim=c(0.4,0.5),col=2,main="Red noise")
plot(rho~theta,data=smap_l,type="l",xlab=expression(theta),ylab=expression(rho),ylim=c(0.6,1),col=4,main="Logistic map")

## The optimal theta determined by maximizing rho
(the_r <- smap_r[which.max(smap_r$rho),"theta"][1])
(the_l <- smap_l[which.max(smap_l$rho),"theta"][1])

References

Hsieh C.H, Glaser S.M., Lucas A.J. & Sugihara G. (2005). Distinguishing random environmental fluctuations from ecological catastrophes for the North Pacific Ocean. Nature, 435, 336-340.(PDF)