mgtsdemo.m

Package [view]: fuzzy.rar
Upload User: hnchenxi
Upload Date: 2008-11-02
Package Size: 1083k
Code Size: 5k
Category: AI-NN-PR
Development Platform: Matlab
  1. %% Chaotic Time Series Prediction
  2. % Chaotic time series prediction using ANFIS.
  3. %
  4. % Copyright 1994-2002 The MathWorks, Inc. 
  5. % $Revision: 1.15 $
  6. %%
  7. % The Mackey-Glass time-delay differential equation is defined by
  8. %  dx(t)/dt = 0.2x(t-tau)/(1+x(t-tau)^10) - 0.1x(t)
  9. % When x(0) = 1.2 and tau = 17, we have a non-periodic and
  10. % non-convergent time series that is very sensitive to
  11. % initial conditions. (We assume x(t) = 0 when t < 0.)
  12. load mgdata.dat
  13. a=mgdata;
  14. time = a(:, 1);
  15. ts = a(:, 2);
  16. plot(time, ts);
  17. xlabel('Time (sec)'); ylabel('x(t)');
  18. title('Mackey-Glass Chaotic Time Series');
  19. %%
  20. % Now we want to build an ANFIS that can predict x(t+6) from
  21. % the past values of this time series, that is, x(t-18), x(t-12),
  22. % x(t-6), and x(t). Therefore the training data format is
  23. %
  24. %  [x(t-18), x(t-12), x(t-6), x(t); x(t+6]
  25. %
  26. % From t = 118 to 1117, we collect 1000 data pairs of the above
  27. % format. The first 500 are used for training while the others
  28. % are used for checking. The plot shows the segment of the time
  29. % series where data pairs were extracted from.
  30. trn_data = zeros(500, 5);
  31. chk_data = zeros(500, 5);
  32. % prepare training data
  33. start = 101;
  34. trn_data(:, 1) = ts(start:start+500-1);
  35. start = start + 6;
  36. trn_data(:, 2) = ts(start:start+500-1);
  37. start = start + 6;
  38. trn_data(:, 3) = ts(start:start+500-1);
  39. start = start + 6;
  40. trn_data(:, 4) = ts(start:start+500-1);
  41. start = start + 6;
  42. trn_data(:, 5) = ts(start:start+500-1);
  43. % prepare checking data
  44. start = 601;
  45. chk_data(:, 1) = ts(start:start+500-1);
  46. start = start + 6;
  47. chk_data(:, 2) = ts(start:start+500-1);
  48. start = start + 6;
  49. chk_data(:, 3) = ts(start:start+500-1);
  50. start = start + 6;
  51. chk_data(:, 4) = ts(start:start+500-1);
  52. start = start + 6;
  53. chk_data(:, 5) = ts(start:start+500-1);
  54. index = 118:1117+1; % ts starts with t = 0
  55. plot(time(index), ts(index));
  56. axis([min(time(index)) max(time(index)) min(ts(index)) max(ts(index))]);
  57. xlabel('Time (sec)'); ylabel('x(t)');
  58. title('Mackey-Glass Chaotic Time Series');
  59. %%
  60. % We use GENFIS1 to generate an initial FIS matrix from training
  61. % data. The command is quite simple since default values for
  62. % MF number (2) and MF type ('gbellmf') are used:
  63.  
  64. fismat = genfis1(trn_data);
  65. % The initial MFs for training are shown in the plots.
  66. for input_index=1:4,
  67.     subplot(2,2,input_index)
  68.     [x,y]=plotmf(fismat,'input',input_index);
  69.     plot(x,y)
  70.     axis([-inf inf 0 1.2]);
  71.     xlabel(['Input ' int2str(input_index)]);
  72. end
  73. %%
  74. % There are 2^4 = 16 rules in the generated FIS matrix and the
  75. % number of fitting parameters is 108, including 24 nonlinear
  76. % parameters and 80 linear parameters. This is a proper balance
  77. % between number of fitting parameters and number of training
  78. % data (500). The ANFIS command looks like this:
  79. %  >> [trn_fismat,trn_error] = anfis(trn_data, fismat,[],[],chk_data)
  80. % To save time, we will load the training results directly.
  81. %
  82. % After ten epochs of training, the final MFs are shown in the
  83. % plots. Note that these MFs after training do not change
  84. % drastically. Obviously most of the fitting is done by the linear
  85. % parameters while the nonlinear parameters are mostly for fine-
  86. % tuning for further improvement.
  87. % load training results
  88. load mganfis
  89. % plot final MF's on x, y, z, u
  90. for input_index=1:4,
  91.     subplot(2,2,input_index)
  92.     [x,y]=plotmf(trn_fismat,'input',input_index);
  93.     plot(x,y)
  94.     axis([-inf inf 0 1.2]);
  95.     xlabel(['Input ' int2str(input_index)]);
  96. end
  97. %% Error curves 
  98. % This plot displays error curves for both training and
  99. % checking data. Note that the training error is higher than the
  100. % checking error. This phenomenon is not uncommon in ANFIS
  101. % learning or nonlinear regression in general; it could indicate
  102. % that the training process is not close to finished yet.
  103. % error curves plot
  104. epoch_n = 10;
  105. tmp = [trn_error chk_error];
  106. subplot(1,1,1)
  107. plot(tmp);
  108. hold on; plot(tmp, 'o'); hold off;
  109. xlabel('Epochs');
  110. ylabel('RMSE (Root Mean Squared Error)');
  111. title('Error Curves');
  112. axis([0 epoch_n min(tmp(:)) max(tmp(:))]);
  113. %% Comparisons
  114. % This plot shows the original time series and the one predicted
  115. % by ANFIS. The difference is so tiny that it is impossible to tell
  116. % one from another by eye inspection. That is why you probably
  117. % see only the ANFIS prediction curve. The prediction errors
  118. % must be viewed on another scale.
  119. input = [trn_data(:, 1:4); chk_data(:, 1:4)];
  120. anfis_output = evalfis(input, trn_fismat);
  121. index = 125:1124;
  122. plot(time(index), [ts(index) anfis_output]);
  123. xlabel('Time (sec)');
  124. %% Prediction errors of ANFIS
  125. % Prediction error of ANFIS is shown here. Note that the scale
  126. % is about a hundredth of the scale of the previous plot.
  127. % Remember that we have only 10 epochs of training in this case;
  128. % better performance is expected if we have extensive training.
  129. diff = ts(index)-anfis_output;
  130. plot(time(index), diff);
  131. xlabel('Time (sec)');
  132. title('ANFIS Prediction Errors');