Home > wafo > wdemos > rfcdemo1.m

rfcdemo1

PURPOSE ^

Demo for switching AR(1)-processes.

SYNOPSIS ^

[F_RFC] = refdemo1(demoNr,P,A,m,s2,param)

DESCRIPTION ^

  RFCDEMO1  Demo for switching AR(1)-processes.  
  
  [F_RFC] = rfcdemo1; 
  [F_RFC] = rfcdemo1(demoNr); 
  [F_RFC] = rfcdemo1(demoNr,P); 
  [F_RFC] = rfcdemo1(demoNr,P,A,m,s2,param); 
  
  Input: 
  demoNr = 1 / 2 / 3 / 0 
           Default: 1. 0 = Define your own process. 
  P      = Transition matrix for regime process.  [rxr] 
  A      = AR coefficients.                       [rx2] 
  m      = Mean koefficients.                     [rx1] 
  s2     = Innovation variances.                  [rx1] 
  param  = Parameter vector, [a b n], defines discretization. 
  
  Output: 
  F_RFC  = Calculated rainflow intensity          [nxn] 
  
    Demo for switching AR(1)-processes, by using  
    switching Markov chains.  
    See Examples 3.1-3.3 in PhD-thesis: 
      P. Johannesson (1999): Rainflow Analysis of 
      Switching Markov Loads. 
      [http://www.maths.lth.se/matstat/staff/pj/PhD/] 
  
  Example: 
    rfcdemo1; 
    rfcdemo1(2); 
    P=[0.98 0.02; 0.05 0.95]; rfcdemo1(1,P); 
  
  See also  rfcdemo2, fatigue

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [F_RFC] = refdemo1(demoNr,P,A,m,s2,param) 
002  
003 % RFCDEMO1  Demo for switching AR(1)-processes.  
004 % 
005 % [F_RFC] = rfcdemo1; 
006 % [F_RFC] = rfcdemo1(demoNr); 
007 % [F_RFC] = rfcdemo1(demoNr,P); 
008 % [F_RFC] = rfcdemo1(demoNr,P,A,m,s2,param); 
009 % 
010 % Input: 
011 % demoNr = 1 / 2 / 3 / 0 
012 %          Default: 1. 0 = Define your own process. 
013 % P      = Transition matrix for regime process.  [rxr] 
014 % A      = AR coefficients.                       [rx2] 
015 % m      = Mean koefficients.                     [rx1] 
016 % s2     = Innovation variances.                  [rx1] 
017 % param  = Parameter vector, [a b n], defines discretization. 
018 % 
019 % Output: 
020 % F_RFC  = Calculated rainflow intensity          [nxn] 
021 % 
022 %   Demo for switching AR(1)-processes, by using  
023 %   switching Markov chains.  
024 %   See Examples 3.1-3.3 in PhD-thesis: 
025 %     P. Johannesson (1999): Rainflow Analysis of 
026 %     Switching Markov Loads. 
027 %     [http://www.maths.lth.se/matstat/staff/pj/PhD/] 
028 % 
029 % Example: 
030 %   rfcdemo1; 
031 %   rfcdemo1(2); 
032 %   P=[0.98 0.02; 0.05 0.95]; rfcdemo1(1,P); 
033 % 
034 % See also  rfcdemo2, fatigue 
035  
036 % Tested  on Matlab  5.3 
037 % 
038 % History: 
039 % Created by PJ (Pär Johannesson) 1997 
040 % Updated by PJ 19-May-2000 
041 %   updated for WAFO 
042 % Updated by PJ 07-Jul-2005 
043 %   warning off 
044  
045 warning_state=warning;  % Get current warning-settings 
046 warning off 
047  
048 fprintf(1,'Welcome to demo 1 for "Rainflow Cycles for Switching Processes"!\n'); 
049 fprintf(1,'It demonstrates calculation of the rainflow intensity\n'); 
050 fprintf(1,'for a switching AR(1)-process.\n'); 
051 fprintf(1,'Type "help rfcdemo1" for further information\n\n'); 
052  
053 % Check input arguments 
054  
055 ni = nargin; 
056 no = nargout; 
057 error(nargchk(0,6,ni)); 
058  
059 Zstr='123456789'; 
060  
061 % Add path to rfcdemo1 
062 demoDir = [waforoot filesep 'wdemos' filesep 'rfcdemo1']; 
063 addpath(demoDir); 
064  
065 % Delete all figure windows 
066  
067 slut=0; 
068 h=gcf; 
069 while ~slut 
070   h_prev=h; 
071   delete(h); 
072   h=gcf; 
073   if h==h_prev, slut=1; end 
074 end 
075    
076  
077 if ni == 0, demoNr=1; end 
078  
079 % Define switching AR(1)-process 
080  
081 if demoNr == 1 
082   p1=0.02; p2=0.01; 
083   P0 = [1-p1 p1; p2 1-p2]; 
084   A = [1 -0.5; 1 0.5]; 
085   m = 2*[-0.5 1.5]'; 
086   s2 = [1 1]'; 
087   param = [-8 8 64]; 
088 elseif demoNr == 2 
089   p1=0.01; p2=0.02; 
090   P0 = [1-p1 p1; p2 1-p2]; 
091   A = [1 -0.5; 1 0.5]; 
092   m = [0 0]'; 
093   s2 = [4 1]'; 
094   param = [-10 10 64]; 
095 elseif demoNr == 3 
096   p12=0.02; p13=0.005; 
097   p21=0.01; p23=0.01; 
098   p31=0.005; p32=0.02; 
099   P0 = [1-p12-p13 p12 p13; p21 1-p21-p23 p23; p31 p32 1-p31-p32]; 
100   A = [1 -0.5; 1 -0.3; 1 0.5]; 
101   m = 2*[-0.5 0 1.5]'; 
102   s2 = [1 1 1.44]'; 
103   param = [-8 8 64]; 
104 elseif demoNr == 0 
105   if ni < 6 
106     error('demoNr=0: Too few input arguments.'); 
107   end 
108 end 
109  
110 if ~exist('P') 
111   P=P0; 
112 end 
113  
114 % Initiation 
115  
116 u=levels(param);    % Discretization levels 
117 delta = u(2)-u(1);  % Discretization step 
118 r=length(P);        % Number of regime states 
119 n=param(3);         % Number of discretization levels 
120 C=ones(r,1);        %  
121  
122 % Calculate statistics: mean , std, proportion, mean time 
123  
124 mX = m./sum(A')'; 
125 sX = sqrt(s2./(1-A(:,2).^2)); 
126 statP = mc2stat(P); 
127 tau = 1./(1-diag(P))'; 
128  
129 % Parameters 
130  
131 fprintf('Parameters\n'); 
132 fprintf(' z      a(z)      m(z)      s(z)\n'); 
133 fprintf('%2d  %8g  %8g  %8g\n',[(1:r)' A(:,2) m sqrt(s2)]') 
134 P 
135 fprintf('\n') 
136  
137 % Statistics 
138 fprintf('Statistics\n'); 
139 fprintf('\n z    m_X(z)    s_X(z)     ro(z)    tau(z)\n'); 
140 fprintf('%2d  %8g  %8g  %8g  %8g\n',[(1:r)' mX sX statP' tau']') 
141 fprintf('\n') 
142  
143 % Calculate spectral densities 
144  
145 fprintf(1,' -- Calculate spectral densities.\n'); 
146  
147 set(gcf,'Name','Power spectra') 
148  
149 for i=1:r 
150   Si = armaspec(1,A(i,:),s2(i,:),100); 
151   subplot(r,1,i) 
152   plot(Si(:,1),Si(:,2)), grid 
153   ylabel(['S(f), regime ' Zstr(i)]) 
154 end 
155 xlabel('frequency, f') 
156  
157 drawnow, disp('Hit any key to continue.'); pause; 
158  
159 % Simulate MARX(2,1)-process 
160  
161 fprintf(1,' -- Simulating.\n'); 
162  
163 T=10000;   % Length of simulation 
164 Tinit=500; % Length of initiation (to remove transients) 
165  
166 [x,z,e] = sarmasim(C,A,m,s2,P,T,Tinit); 
167  
168 figure 
169 set(gcf,'Name','Switching AR(1)-process') 
170  
171 I = 1:501; 
172 hmmplot(x(I),z(I),I-1,[1 r],'','Switching Process, X(k)',1); 
173 ylabel('regime, Z(k)'),xlabel('time, k') 
174  
175 drawnow, disp('Hit any key to continue.'); pause; 
176  
177 % 
178 % Discretize AR(1)-processes 
179 % Approximate with Markov chain 
180 % 
181  
182 ud = u; ud(1)=u(1)-0.25*(u(n)-u(1)); ud(n)=u(n)+0.25*(u(n)-u(1)); 
183  
184 figure 
185 set(gcf,'Name','Markov chain for Discretized AR(1)-processes') 
186  
187 if demoNr == 1 | demoNr ==2 | demoNr ==3 
188   fprintf(1,' -- Loading discretized AR-processes.\n'); 
189  
190   eval(['load rfcdem1' Zstr(demoNr)]); % load uu,Q1,Q2,Q3 
191   for i=1:r 
192     subplot(1,r,i) 
193     eval(['Q{i}=Q' num2str(i) ';']); 
194     cmatplot(uu,uu,Q{i},3), axis square 
195     title(['Regime ' Zstr(i)]) 
196   end; 
197   drawnow 
198 else 
199  
200   fprintf(1,' -- Discretizing AR-processes. (This can take some time!)\n'); 
201   for i=1:r 
202     fprintf(1,' -- Discretizing process %d.\n',i); 
203     a1i = A(i,2); 
204     si = sqrt(s2(i)); 
205     mi = mX(i); 
206     [Qi,uu] = discar(ud,a1i,mi,si); 
207     Q{i}=Qi; 
208     %eval(['Q{i}' num2str(i) '=Qi;']); 
209     subplot(r,1,i) 
210     cmatplot(uu,uu,Qi,3), axis square 
211     title(['Regime ' Zstr(i)]) 
212     drawnow 
213   end 
214 end 
215  
216 drawnow, disp('Hit any key to continue.'); pause; 
217  
218 % 
219 % Compute RFC counting distribution for SMC 
220 % 
221  
222 fprintf(1,' -- Calculating rainflow intensity.\n'); 
223  
224 Qstr=''; 
225 for i=1:r 
226   Qstr = [Qstr ',Q' Zstr(i)]; 
227 end 
228  
229 [F_RFC,mu_RFC] = smc2rfm(P,Q,[2 delta]);    % definition 2 
230 %eval( ['[F_RFC,mu_RFC] = smc2rfc(P,[2 delta]' Qstr ');'] );    % definition 2 
231  
232 for i=1:r 
233   %eval( ['[F_RFC' Zstr(i) ',mu_RFC' Zstr(i) '] = mc2rfc(Q' Zstr(i) ',[2 delta]);'] );     % Regime i 
234   [F_RFC0{i},mu_RFC0{i}] = mc2rfm(Q{i},[2 delta]);     % Regime i 
235 end 
236  
237 figure 
238  
239 cmatplot(u,u,F_RFC,2); 
240 set(gcf,'Name','Rainflow intensity for switching AR(1)-process') 
241 title('3D-plot') 
242  
243 drawnow, disp('Hit any key to continue.'); pause; 
244  
245 figure 
246  
247 % Calculate Cycles and crossings 
248  
249 fprintf(1,' -- Cycles and crossings for simulated load.\n'); 
250  
251 TP = dat2tp([(1:T)' x]); 
252 RFC = tp2rfc(TP); 
253 cross = cc2lc(RFC); 
254  
255 cocc(param,RFC,F_RFC) 
256 set(gcf,'Name','Rainflow intensity for switching AR(1)-process') 
257 title('Iso-lines compared with observed cycles') 
258 xlabel('min'), ylabel('Max') 
259  
260 drawnow, disp('Hit any key to continue.'); pause; 
261  
262 % Plot Crossing intensity 
263  
264 figure 
265 set(gcf,'Name','Crossing intensity') 
266  
267 subplot(1,2,1) 
268 plot(u,diag(mu_RFC),'g--'), hold on 
269 stairs(cross(:,1),cross(:,2)/T), hold off 
270 title('Calculated vs. observed') 
271 xlabel('level, u'); ylabel('Crossing intensity, mu(u)'); 
272  
273 subplot(1,2,2) 
274 plot(u,diag(mu_RFC),'g'), hold on 
275 cross_sum=zeros(size(diag(mu_RFC))); 
276 for i=1:r 
277   crossi=diag(mu_RFC0{i}); 
278   plot(u,statP(i)*crossi,'r-.') 
279   cross_sum=cross_sum+statP(i)*crossi; 
280 end 
281 plot(u,cross_sum,'y--'), hold off 
282 title('Calculated vs. individual') 
283 xlabel('level, u'); ylabel('Crossing intensity, mu(u)'); 
284  
285 fprintf(1,'\nThank you for running this demo!\n'); 
286  
287 % Remove path to rfcdemo1 
288 demoDir = [waforoot filesep 'wdemos' filesep 'rfcdemo1']; 
289 rmpath(demoDir); 
290  
291 warning(warning_state); 
292  
293

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003