奇异谱分析脚本参考

qxsw2016

利用1951-1999年浙江的降水指数，进行SSA分析，窗口长度取为M = 10,并重建成分RC1,RC3+RC4.
def ReComponents(n,rc): #重建成分
    m = rc.shape[0]
    p = rc.shape[1]
    y = zeros(n)
    for k in arange(n):
        if(k<m-1):
            y_m = 0        
            for i in arange(m):
                for j in arange(p):
                    if((i+j)==k):
                        y_m = y_m + rc[i,j]
            y[k] = y_m/(k+1)
        elif(k>=m-1 and k<=n-m):
            y_m = 0
            for i in arange(m):
                for j in arange(p):
                    if((i+j)==k):
                         y_m = y_m + rc[i,j]
            y[k] = y_m/m
        else:
            y_m = 0
            for i in arange(m):
                for j in arange(p):
                    if((i+j)==k):
                        y_m = y_m +rc[i,j]
            y[k] = y_m/(n-k)
    return y

fn= 'F:\dataset\ZJJJA.txt'
ncol = numasciicol(fn)
nrow =  numasciirow(fn)
#print ncol,nrow
fndata = asciiread(fn,shape=(nrow,ncol))
RJJA = fndata[:,0] #标准化数据
m= 10 #窗口长度或嵌入空间维数
n = 49  #1951-1999
p = n-m+1
S = zeros((m,m))#Toeplitz矩阵
x = zeros((m,p))
c = zeros(m)
#print RJJA
for i in arange(m):
    x[i,:] = RJJA[i:i+p] #时迟排列相空间矩阵
for i in arange(m):
    c = sum(RJJA[0:n-i]*RJJA[i:n])/(n-i)
for i in arange(m):
    for j in arange(i,m):
        S[i,j] = c[j-i] #右上
for i in arange(1,m):
    for j in arange(0,i):
        S[i,j] = S[j,i]
lamda,V = linalg.eig(S)
TC = np.dot(V.T, x)
V = V[:,::-1]
TC = TC[::-1,:]
lamda = lamda[::-1]
#print S

#第一个重建成分rc1
rc1 = np.dot(V[:,0].reshape(m,1),TC[0,:].reshape(1,p))
yk = ReComponents(n,rc1)
xtick = zeros([n])
for i in arange(n):
    xtick = i+1951
subplot(2,1,1)
plot(xtick,yk,'b-')
title('RC1')
#第三、四个重建成分rc3、rc4
rc3 = np.dot(V[:,2].reshape(m,1),TC[2,:].reshape(1,p))
rc4 = np.dot(V[:,3].reshape(m,1),TC[3,:].reshape(1,p))
yk = ReComponents(n,rc3)+ReComponents(n,rc4)
subplot(2,1,2)
plot(xtick,yk,'b-')
title('RC3+RC4')

movlan

学习了··谢谢分享！

MeteoInfo

谢谢分享！已加入脚本汇总贴。

qxsw2016

最大熵谱脚本参考：
以奇异谱分析的重建成分为基础，基于Burg 算法的最大熵谱估计（数组索引开始位置为1，有点晕），过程中用到预测误差准则FPE确定最佳回归阶数。rc3+rc4的最大熵谱分析如下：
fn= r'F:\dataset\rc34.txt'
ncol = numasciicol(fn)
nrow =  numasciirow(fn)
fndata = asciiread(fn,shape=(nrow,ncol))
x = fndata[:,0]
n = 49
m = 24 #最大波数
delta2 = zeros(n+1)
a = zeros((n+1,n+1)) #处理方便，递推公式索引从1开始
r = zeros(n+1)
fpe = zeros(n-1)
shl = zeros(m+1)
#基于Burg 算法的最大熵谱估计
# k = 0
qsum = sum(x*x)
delta2[0] = qsum/float(n)
r[0]= delta2[0]
# k = 1
s12 = sum(x[0:n-1]*x[1:n])
s11 = sum(x[0:n-1]*x[0:n-1])
s22 = sum(x[1:n]*x[1:n])

a[1,1] = 2*s12/(s11+s22)        
delta2[1]=(1-a[1,1]*a[1,1])*delta2[0]
r[1] = a[1,1]*delta2[0]
# k = 2
x1 = x[2:n]-a[1,1]*x[1:n-1]
x2 = x[0:n-2]-a[1,1]*x[1:n-1]
s12 = sum(x1*x2)
s11 = sum(x1*x1)
s22 = sum(x2*x2)

a[2,2] = 2*s12/(s11+s22)
a[2,1] = a[1,1]-a[2,2]*a[1,1]
delta2[2] = (1-a[2,2]*a[2,2])*delta2[1]
r[2] = a[1,1]*r[1]+a[2,2]*delta2[1]
# k >2
for k in arange(2,n-1):      
    s12= 0   
    s11 = 0
    s22 = 0         
    for i in arange(k,n):
        v1 =0
        v2 =0
        for j in arange(1,k+1):
            v1 = v1 + a[k,j]*x[i-j]
            v2 = v2 + a[k,j]*x[i-k+j-1]            
        v1 =x - v1  
        v2 =x[i-k-1] - v2              
        s12 = s12 + 2*v1*v2   
        s11 = s11 + v1*v1
        s22 = s22 + v2*v2
    a[k+1,k+1] = s12/(s11+s22)

    for j in arange(1,k+1):
        a[k+1,j] = a[k,j]-a[k+1,k+1]*a[k,k+1-j]
    delta2[k+1] = (1-a[k+1,k+1]*a[k+1,k+1])*delta2[k]

    v1 = sum(a[k,1:k+1]*r[k+1:1:-1])
    r[k+1] = v1+a[k+1,k+1]*delta2[k]
fpe[0] = 999999
#fpe 准则确定最佳回归阶数k0
for k in arange(1,n):
    fpe[k-1] = (n+k)*delta2[k]/(n-k)
k0 = fpe.index(min(fpe))
print k0
print delta2[11]
for l in arange(0,m+1):
    v1 = 0
    v2 = 0
    for k in arange(1,k0+1):
        v1 = v1 + a[k0,k]*cos(pi*l*k/m)
        v2 = v2 + a[k0,k]*sin(pi*l*k/m)
    v1 = (1-v1)*(1-v1)
    v2 = v2*v2
    shl[l] = delta2[k0]/(v1+v2)
xtick = arange(0,m+1)
xtop = zeros(m+1)
xtop[1:m+1] = 2*m*1./xtick[1:m+1]  

ax1=axes()
plot(xtick,shl,'blacks--',linewidth=2)   
xticks(xtick)   
xlabel('L',fontsize=16)
xlim(0,m+1)
ax2 = twiny(ax1)
xaxis(location='top',tickin=False)
xlabel('T(year)',bold=False)
xlim(0,m+1)
xticks(xtick, ['%.1f' % xx for xx in xtop])
title('RC3+RC4 maximum entropy spectra')