利用互信息量法确定延迟时间T
###确定延迟时间tau——互信息量法
MI<-function(data, tau_max, n){
I_sq = rep(NA,30) # 建立一个空数组,用来存放每个τ值的互信息
data_length = length(data) # 计算时间序列的长度
for (tau in 31: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-30] = H_s+H_q-H_sq # 计算tau下的互信息函数
}
return(I_sq)
}
选择互信息函数的第一个极小值点对应的T作为最优延迟时间
|