% Exercise 2: AR models

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ndata　と　M　の値を変えて試して見よ
% 学習データ数
Ndata = 50;
% AR 次数
M = 5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
 % 学習すべきデータ x の生成
 %  x は 長さ2Nの時系列で、2N次元の縦ベクトル。
 %  randn()は、正規ノイズを成分とする行列を出力する関数。
 N = Ndata;
 t = (1:N*2)'/100;
 x = sin(t*4*pi) ...
    + sin(t*4.1*pi) ...
    + randn(N*2,1)*0.1;
 % Lorenz attractor のデータを使うならば、
 % 下の行の'%'をはずすべし。
   %load lorenz

 % パラメータベクトル推定
 %  ここでx(1:N)はxベクトルの前半分。
 a = AR_estimate( x(1:N), M );

 % パラメータベクトルと x の前半分に基づいて、
 %  x の後ろ半分をマルチステップ予測
 xp = AR_predict( x(1:N), a, N );
 
 % 結果のプロット
 plot( x,'.' );
 hold on
 plot( N+1:N*2, xp, 'r' )
 hold off
 title( ['AR sample  degree=' num2str(M)] )