R语言实现相空间重构的参数选取

前言

最近写论文用到了相空间重构（PSR）技术，该方法简单来说就是将一个一维的时间序列通过重构的方法映射为一个矩阵，且该矩阵保留了原始时间序列的特征。
进行PSR的关键是确定两个参数：延迟时间和嵌入维数，关于两个参数的选取对重构的影响，这儿不多做说明，可以查阅文献。简单来说两个参数的选择直接关系到重构效果的好坏。
下面给出用互信息量法（AMI）求延迟时间的R代码

R语言求解延迟时间

###确定延迟时间tau——互信息量法（基于等间距格子法）
MI<-function(data, tau_max, n){
  I_sq = rep(NA,tau_max)  # 建立一个空数组，用来存放每个τ值的互信息
  data_length = length(data) # 计算时间序列的长度
  for (tau in 1:tau_max){
    S = data[1:(data_length-tau)]
    Q = data[(tau+1):data_length]
    as1 = min(S)
    bq = min(Q)
    delts = (max(S)-as1)/n
    deltq = (max(Q)-bq)/n
    N_sq = matrix(0,ncol=n,nrow = n)
    for (i in 1:n){
      for (j in 1:n){
        for (k in 1:(data_length-tau)){
          as_k = (S[k]-as1)/delts
          bq_k = (Q[k]-bq)/deltq
          if (as_k >= i-1 && as_k < i && bq_k >=j-1 && bq_k < j){
            N_sq[i, j] = N_sq[i, j]+1
          }
        }
      }
    }
    Ntotal = sum(N_sq)
    Ps<-rep(0,n)
    Pq<-rep(0,n)
    for (i in 1:n) {
      Ps[i] = sum(N_sq[,i])/Ntotal   # 计算位于一维s格子内的概率
      Pq[i] = sum(N_sq[i,])/Ntotal  # 计算位于一维q格子内的概率
    }
    Psq = N_sq/Ntotal         # 计算位于二维格子(ii,jj)内概率
    H_s = 0  # 计算S的熵
    H_q = 0  # 计算q的熵
    for (i in 1:n){
      if (Ps[i] != 0)
      {H_s = H_s - Ps[i] * log(Ps[i])}
      if (Pq[i] != 0)
      {H_q = H_q - Pq[i] * log(Pq[i])}
    }
   
    H_sq = 0  # 计算(s,q)的联合熵
    for (i in 1:n){
      for (j in 1:n){
        if (Psq[i, j] != 0)
        {H_sq = H_sq - Psq[i, j] * log(Psq[i, j])}
      }
    }
   
    I_sq[tau] = H_s+H_q-H_sq  # 计算tau下的互信息函数
  }
  return(I_sq)
  }

选择互信息函数的第一个极小值点对应的T作为最优延迟时间。

后来通过在网上搜索发现，R里面有现成的包可以直接确定两个参数，而且运算速度比我编的代码不知道快多少（泪目，自己的代码跑了两三天）。

延迟时间的确定——R语言

R里面nonlinearTseries包里面的timeLag函数可以用于选择延迟时间tau
参数technique = "acf"时，表示用自相关函数法，此时参数selection.method选择默认
参数technique = “ami”,表示用互信息量法，此时参数selection.method = “first.minimum”

嵌入维数的确定——R语言

嵌入维数m可用函数estimateEmbeddingDim()函数，此时使用Cao方法确定维数的
estimateEmbeddingDim(data, time.lag = tau,do.plot = F,max.embedding.dim = 20)
tseriesChaos包中的false.nearest()函数使用伪近邻方法确定嵌入维数
false.nearest(data, m = max.m, d = tau, t = 150, eps = sd(data)/10)

小结

以上就是关于用R语言做PSR的相关代码，希望能帮到那些再用这个方法的人吧，少走一些弯路。毕竟自己编代码太痛苦了，还是直接调包比较香（开玩笑啦，有能力的还是尝试自己写代码吧，对自己提升很大）