最近写论文用到了相空间重构(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的相关代码,希望能帮到那些再用这个方法的人吧,少走一些弯路。毕竟自己编代码太痛苦了,还是直接调包比较香(开玩笑啦,有能力的还是尝试自己写代码吧,对自己提升很大)

相空间重构求关联维数——GP算法、自相关法求时间延迟tau、最近邻算法求嵌入维数m GP算法: 若有一维时间序列为{x1,x2,…,xn},对其进行相空间重构得到高维相空间的一系列向量: xi(τ,m)=(xi,xi1,⋯ ,xi+(m−1)τ){x_i}(\tau ,m) = \left( {{x_i},{x_{i1}}, \cdots ,{x_{i + {{(m - 1)}_\tau }}}} \right)xi​(τ,m)=(xi​,xi1​,⋯,xi+(m−1)τ​​) 式中:τ\tauτ为时间延迟 布谷鸟搜索(Cuckoo Search, CS)算法是一种新型群体智能优化算法,不仅结合了许多鸟类及果蝇特殊的利维飞行模式进行搜索,而且增加了群体之间的信息交流,加快收敛速度,为BP神经网络机参数优化提供了一种新的研究工具。为了提高电力系统负荷和天气预测精度,提出一种CS算法优化BP神经网络参数的负荷和天气预测预测模型(CS-BP)。 首先,输入的时间序列向量被重构为一个低维度矩阵 X,其中每行代表一个点,每列代表一个特征。然后,可以对矩阵进行描点和缩放。接下来,可以使用主成分分析(PCA)对矩阵进行降维,以便更好地可视化和理解数据。最后,可以绘制主成分贡献度图和二维散点图来展示数据的主要特征。需要注意的是,相空间重构是一个复杂的过程,需要根据具体的数据和应用场景进行调整和优化。上述代码仅提供了一个基本的框架,需要根据实际需求进行修改和扩展。它可以将高维时间序列数据转换为低维表示,以便更好地理解和分析系统的行为。 1)异常模式发生时,程序把通用寄存器压入堆,SP一直指向栈顶的位置。返回时再出栈,保证程序状态的完整性。 2)有MSP 和PSP(两者只需一个,不能同时使用,默认MSP)。 MSP :主堆栈指针,系统用。PSP : 进程堆栈指针,个人堆栈指针。 R14:LR连接寄存器 功能:保... 互信息法求时延参数 matlab 对于一维时间序列,可以用时延嵌入的方式扩展到高维空间来学习其动态特性。时延嵌入的参数:一个是时间延迟 τ,一个是嵌入维数 m。时间延迟参数常用互信息法求得。 互信息(mutual information):原始时间序列x 与时延τ之后的时间序列 x(t+τ) 之间的依赖关系。公式如下: 1. MI 计算函数 function mi = MI(x, nbins, maxlag) % 时间序列x, 分箱数nbins,最大时延maxlag mi = zero