Реконструировать временные ряды из SSA

Question

Давайте рассмотрим следующий код

clear all;
  B=xlsread('data_generations1','A1','g8:g301');
  n=length(B);
  L =input('Give the size of the interval: ' );% Number of columns in the Data matrix
   m=n-L+1;%number of rows in the Data matrix
 X = zeros(m,L);
  for i=1:m
        X(i,:)=B(i:i+L-1);
  end;
   S=X*X'; 
    [U,autoval]=eig(S);
    [d,i]=sort(-diag(autoval));  
   d=-d;
   U=U(:,i);sev=sum(d); 
    plot((d./sev)*100),hold on,plot((d./sev)*100,'rx');
    title('Singular Spectrum');xlabel('Eigenvalue Number');ylabel('Eigenvalue (% Norm of trajectory matrix retained)')
   V=(X')*U; 
   rc=U*V';
% Step 3: Grouping
   I=input('Choose the agrupation of components to reconstruct the series in the form I=[i1,i2:ik,...,iL]  ')
   Vt=V';
   rca=U(:,I)*Vt(I,:);
% Step 4: Reconstruction

   y=zeros(n,1);  
   Lp=min(L,m);
   Kp=max(L,m);

   for k=0:Lp-2
     for m1=1:k+1;
      y(k+1)=y(k+1)+(1/(k+1))*rca(m1,k-m1+2);
     end
   end

   for k=Lp-1:Kp-1
     for m1=1:Lp;
      y(k+1)=y(k+1)+(1/(Lp))*rca(m1,k-m1+2);
     end
   end

   for k=Kp:n
      for m1=k-Kp+2:n-Kp+1;
       y(k+1)=y(k+1)+(1/(n-k))*rca(m1,k-m1+2);
     end
   end

   figure;subplot(2,1,1);hold on;xlabel('Data poit');ylabel('Original and reconstructed series')
   plot(x1);grid on;plot(y,'r')
   r=x1-y;
   subplot(2,1,2);plot(r,'g');xlabel('Data poit');ylabel('Residual series');grid on
   vr=(sum(d(I))/sev)*100;

Фрагмент кода был использован с этого сайта:

http://www.mathworks.com/matlabcentral/fileexchange/8115-singular-spectrum-analysis-smoother/content/ssa.m

Когда я запускаю этот код averaging

Give the size of the interval: 15
Choose the agrupation of components to reconstruct the series in the form I=[i1,i2:ik,...,iL]  4

I =

     4

Attempted to access rca(1,16); index out of bounds because size(rca)=[280,15].

Error in averaging (line 37)
      y(k+1)=y(k+1)+(1/(Lp))*rca(m1,k-m1+2);

В чем проблема в моем коде? Пожалуйста, помогите мне.

Answer 1

Как следует из сообщения об ошибке, вы пытаетесь получить доступ к 16-му столбцу rca а также rca только 15 столбцов.

Запустив ваш код в Octave, я получаю следующие значения:

>> I
I =  4
>> Kp
Kp =  280
>> Lp
Lp =  15
>> n
n =  294
>> m
m =  280

Ошибки кода при k=15 а также m1=1, который дает k-m1+2 = 16отсюда и сообщение об ошибке.

РЕДАКТИРОВАТЬ

Вместо того, чтобы изменять функцию, почему бы просто не вызвать ее со своими данными?

clear all;
B=xlsread('data_generations1','A1','g8:g301');
L =input('Give the size of the interval: ' );% Number of columns in the Data matrix
[y,r,vr]=ssa(B,L);

Answer 2

Посмотрите на эту работающую GNU Octave-совместимую базовую реализацию SSA: ssa-octave.m. Это адаптация моего кода Scilab, который отлично работает для меня. Обратите внимание на комментарии внутри, они могут содержать некоторые важные заметки.

Если вам интересно, я также могу предоставить M-SSA (многовариантную) реализацию. Он довольно близок к базовому SSA и отличается только фазами внедрения по времени и диагональному усреднению.

Я хочу сказать, что SSA выглядит для меня как (и я верю, что так) PCA, применяемый для встраивания временных рядов с задержкой.

Answer 3

@ Answeriver ... Я только начал изучать SSA и спектральный анализ в целом. У меня есть рабочая версия 1D SSA, и она отлично работает. Я перехожу к 2D- и 3D-анализу. Если бы вы могли поделиться своим M-SSA, это было бы очень полезно. Ваше здоровье...