特異スペクトル解析による線形再帰関係式の導出がよくない気がする -

特異スペクトル解析(SSA)は、原理的に
$x_n = \sum_{k=1}^p a_k x_{n-k}$
という形の予測式を与えることができる(式の形は、ARモデル/線形予測分析と同一)

Singular spectrum analysis based on the perturbation theory
http://www.cf.ac.uk/maths/subsites/zhigljavskyaa/pdfs/SSA/SSA%20perturbation.pdf

On the choice of parameters in Singular Spectrum Analysis and related subspace-based methods
http://arxiv.org/abs/1005.4374

などを見ると、簡単に書いてある。実装は平易で

import numpy as np

def ssa(dat,Fs,Cmax=5):
   N = len(dat)
   L = N/2
   M = N+1-L
   Kmax = Cmax*2
   H = np.zeros(L*M).reshape((L,M))
   for n in xrange(L):
      for m in xrange(M):
          H[n,m] = dat[n+m]
   U,S,V = np.linalg.svd(H)
   for k,(s1,s2) in enumerate(zip(S,S[1:])):
       if s2<s1*1.0e-3 or k==len(S)-2 or k>=Kmax-1:
          K=k+1
          U_K,V_K = U[: , :K],V[: , :K]
          break
   A = np.zeros( L-1 )
   v2 = 0
   for i in xrange(K):
       pi_i = U_K[-1,i]
       v2 += pi_i**2
       for j in xrange(L-1):
          A[j] += pi_i*U_K[j,i]
   assert(v2<1.0),"SSA failed"
   A= A/(1-v2)
   return A

みたいな感じになる。係数の個数は、L-1であるが、これは本質的には任意で、ものすごく大きくなりうる。そして、正しくは、コード中のKが係数の個数(=モデルの次数-1)になるべきだと思う

import numpy as np
Fs=400.0
dat = np.sin(2*np.pi*123.4*np.arange(100)/Fs)

で与えられる、単純な正弦波を予測する、最良の線形モデルは
$y_n = -0.7186908 y_{n-1} - y_{n-2}$
で、実際に

for n in xrange(2,10):
    print dat[n] , -(0.7186908*dat[n-1]+dat[n-2])

の出力は

-0.670685576537 -0.670685583991
-0.451189084442 -0.451189079084
0.994951016981 0.994951020585
-0.263873049965 -0.263873057913
-0.805307885711 -0.805307883603
0.84264041216 0.842640418593
0.199709980514 0.199709973783
-0.986170136229 -0.986170137824

となって高い精度で一致する(まぁ当たり前だが)
$z^2 + 0.7186908 z + 1 = 0$
の根は、-0.359349+0.93320324j, -0.359349-0.93320324jの2つで、単位円上に乗っている。偏角は、(123.4/200.0)πに等しい。-0.7186908は2cos(123.4π/200.0)から来ている。これは、微分方程式を単純に差分化したものよりよい。

要するに、単純な正弦波は、2階の定数係数線形離散方程式で記述されるが、もっと高階の定数係数線形離散方程式の解にもなりうる。後者の高階方程式の解には、本来知りたかった波とは何の関係もない周波数を持った正弦波が解として含まれうる。こいつは、解析の邪魔でしかない。一般には、信号を、複数の正弦波に分解したいわけで、じゃぁ、これは、何個の正弦波から構成されていると考えるべきかというと、部分空間の次元=Kが求める答えだと思う(sinとcosがあるので、普通、これは偶数)。多分、予測精度としては、高階の定数係数線形離散方程式を使っても、問題ないのだろうけど(逆に言うと、部分空間の次元を大きくとりすぎると、本来いらない周波数成分が見えてしまうということである。振幅・位相は、最小二乗推定で計算するから、これは、現実の信号に対して、振幅・位相を大きくずらしてしまう場合がある)

よりよい実装は例えば以下のようなものと思う(多分)。Kが決まっているなら、L=K+1と置いてしまっても、(少なくとも、人工的なデータに対しては)よいようである

import numpy as np

def ssa(dat,Fs,Cmax=5):
   N = len(dat)
   L = N/2
   M = N+1-L
   Kmax = Cmax*2
   H = np.zeros(L*M).reshape((L,M))
   for n in xrange(L):
      for m in xrange(M):
          H[n,m] = dat[n+m]
   U,S,V = np.linalg.svd(H)
   for k,(s1,s2) in enumerate(zip(S,S[1:])):
       if s2<s1*1.0e-3 or k==len(S)-2 or k>=Kmax-1:
          K=k+1
          U_K,V_K = U[: , :K],V[: , :K]
          break
   U2 = np.zeros( (L-1)*2*K ).reshape( (L-1,2*K) )
   for k1 in xrange(L-1):
       for k2 in xrange(K):
            U2[k1,k2] = U_K[k1,k2]
            U2[k1,k2+K] = U_K[k1+1,k2]
   _,_,V2 = np.linalg.svd(U2)
   V12 = np.matrix(V2[K:,:K])
   V22 = np.matrix(V2[K:,K:])
   Q = -V12*np.linalg.inv(V22)
   return -np.poly(np.matrix(Q))[1:]   #-- 最初は1であるはず