Home > wafo > markov > test > test_markov.m

test_markov

PURPOSE ^

Quick test of the routines in module 'markov'

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

 TEST_MARKOV Quick test of the routines in module 'markov'
 
  This script is used for a quick test of the routines in
  Modul 'markov':
 
  The following routines are tested:
 
  
  MC - Markov chain
  - SIMMC - Simulates a Markov chain with state space {1,2,...,r}.
  - MC2RFM - Calculates the rainflow matrix/intensity for a Markov chain.
  - MC2STAT - Calculates the stationary distribution for a Markov chain.
 
  MCTP - Markov Chains of Turning Points
  - MCTPSIM - Simulates a Markov chain of turning points
  - MCTP2RFM - Calculates the rainflow matrix for a MCTP.
  - MCTP2STAT - Calculates the stationary distribution for a MCTP.
 
  SMC - Switching Markov Chains 
  - SMC2RFM - Calculates the rainflow matrix/intensity for a switching Markov chain.
 
  SMCTP - Switching Markov Chains of Turning Points
  - SMCTPSIM - Simulates a switching Markov chain of turning points,
  - HMMPLOT - plots a Hidden Markov Model.
  - SMCTP2RFM - Calculates the rainflow matrix for a SMCTP.
  - SMCTP2STAT - Stationary distributions for a switching MCTP.
 
  Decomposition of a Mixed rainflow matrix
  - ESTSMCTP - Estimate SMCTP model from an observed rainflow matrix.
 
  Simulation from Rainflow matrix
  - RFM2DTP  Reconstructs a sequence of turning points from a rainflow matrix.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 %TEST_MARKOV Quick test of the routines in module 'markov'
002 %
003 % This script is used for a quick test of the routines in
004 % Modul 'markov':
005 %
006 % The following routines are tested:
007 %
008 % 
009 % MC - Markov chain
010 % - SIMMC - Simulates a Markov chain with state space {1,2,...,r}.
011 % - MC2RFM - Calculates the rainflow matrix/intensity for a Markov chain.
012 % - MC2STAT - Calculates the stationary distribution for a Markov chain.
013 %
014 % MCTP - Markov Chains of Turning Points
015 % - MCTPSIM - Simulates a Markov chain of turning points
016 % - MCTP2RFM - Calculates the rainflow matrix for a MCTP.
017 % - MCTP2STAT - Calculates the stationary distribution for a MCTP.
018 %
019 % SMC - Switching Markov Chains 
020 % - SMC2RFM - Calculates the rainflow matrix/intensity for a switching Markov chain.
021 %
022 % SMCTP - Switching Markov Chains of Turning Points
023 % - SMCTPSIM - Simulates a switching Markov chain of turning points,
024 % - HMMPLOT - plots a Hidden Markov Model.
025 % - SMCTP2RFM - Calculates the rainflow matrix for a SMCTP.
026 % - SMCTP2STAT - Stationary distributions for a switching MCTP.
027 %
028 % Decomposition of a Mixed rainflow matrix
029 % - ESTSMCTP - Estimate SMCTP model from an observed rainflow matrix.
030 %
031 % Simulation from Rainflow matrix
032 % - RFM2DTP  Reconstructs a sequence of turning points from a rainflow matrix.
033 
034 % History:
035 % Created by PJ (Pär Johannesson) 18-May-2000
036 % Updated by PJ (Pär Johannesson) 07-Jul-2005
037 
038 figure(1), clf
039 
040 echo on
041 
042 %MAT2TMAT - Converts a matrix to a transition matrix.
043 
044 F = magic(5)
045 mat2tmat(F)
046 mat2tmat(F,1)
047 mat2tmat(F,-1)
048 mat2tmat(F,-1,-3)
049 mat2tmat(F,-1,0)
050 
051 pause
052 
053 %%%
054 % MC - Markov chain
055 %
056 
057 % Model Definition - MC
058 n = 16; param = [-1 1 n];                % Define discretization
059 u = levels(param);                       % Discrete levels
060 [G,Gh] = mktestmat(param,[-0.2 0.2],0.2,1);  % min-max matrix
061 P = G+Gh;
062 
063 %SIMMC - Simulates a Markov chain with state space {1,2,...,r}.
064 T = 5000;           % Length of simulation (number of TP)
065 xD = mcsim(P,T);  % Simulate
066 x = u(xD)';         % Change scale to levels -1,..,1
067 
068 t=1:200; plot(t,x(t))  
069 
070 pause
071 
072 %MC2RFM - Calculates the rainflow matrix/intensity for a Markov chain.
073 Grfc=mc2rfm(P);
074 tp = rfcfilter(xD,0,1);
075 FrfcObs = dtp2rfm(tp,n);
076 
077 cmatplot(u,u,{G Grfc FrfcObs/(T/2)},13)
078 
079 pause
080 
081 %MC2STAT - Calculates the stationary distribution for a Markov chain.
082 ro = mc2stat(P);
083 
084 clf, plot(u,ro)
085 
086 pause
087 
088 %%%
089 % MCTP - Markov Chains of Turning Points
090 %
091 
092 % Model Definition - MCTP
093 n = 16; param = [-1 1 n];                % Define discretization
094 u = levels(param);                       % Discrete levels
095 [G,Gh] = mktestmat(param,[-0.2 0.2],0.2,1);  % min-max matrix
096 
097 %MCTPSIM - Simulates a Markov chain of turning points
098 T = 5000;                % Length of simulation (number of TP)
099 xD = mctpsim({G []},T);  % Simulate
100 x = u(xD)';              % Change scale to levels -1,..,1
101 
102 t=1:200; plot(t,x(t))  
103 
104 pause
105 
106 %MCTP2RFM - Calculates the rainflow matrix for a MCTP.
107 
108 Grfc=mctp2rfm({G,[]});
109 FrfcObs = dtp2rfm(xD,n);
110 
111 cmatplot(u,u,{G Grfc FrfcObs/(T/2)},13)
112 
113 pause
114 
115 %MCTP2ARFM - Calculates the rainflow matrix for a MCTP.
116 
117 Garfc=mctp2arfm({G,[]});
118 FarfcObs = dtp2arfm(xD,n);
119 
120 cmatplot(u,u,{Garfc FarfcObs/(T/2)},13)
121 
122 pause
123 
124 %MCTP2STAT - Calculates the stationary distribution for a MCTP.
125 [ro_min,ro_max] = mctp2stat({G Gh});
126 [ro_min,ro_max] = mctp2stat({G []});
127 
128 clf, plot(u,ro_min,u,ro_max)
129 
130 pause
131 
132 %%%%
133 % SMC - Switching Markov Chains 
134 % 
135 
136 P=[0.9 0.1; 0.2 0.8]; 
137 param = [-1 1 32]; u = levels(param);
138 [F1,F2] = mktestmat(param,[-0.3 0.3],0.15,1,-Inf);
139 
140 %SIMSMC - Simulates a Switching Markov chain with state space {1,2,...,r}.
141 [x,z] = smcsim(P,{F1 F2},400); 
142 hmmplot(x,z,(1:length(z))',[1 2])         % Same colour
143 
144 pause
145 
146 %SMC2RFM - Calculates the rainflow matrix/intensity for a switching Markov chain.
147 Frfc = smc2rfm(P,{F1;F2});
148 cmatplot(u,u,Frfc)
149 
150 pause
151 
152 %%%%
153 % SMCTP - Switching Markov Chains of Turning Points
154 % 
155 
156 % Model Definition
157 
158 n=16; param = [-1 1 n];                   % Define discretization
159 u=levels(param);                          % Discrete levels
160 G1 = mktestmat(param,[-0.4 -0.3],0.15,1); % regime 1
161 G2 = mktestmat(param,[0.3 0.4],0.15,1);   % regime 2
162 
163 p1=0.10; p2=0.05;
164 P=[1-p1 p1; p2 1-p2]  % Transition matrix
165 statP=mc2stat(P)      % Stationary distribution
166 
167 % SMCTPSIM - Simulates a switching Markov chain of turning points,
168 T=5000;   % Length of simulation (number of TP)
169 [xD,z] = smctpsim(P,{G1 []; G2 []},T); % Simulate
170 x=u(xD)'; % Change scale to levels -1,..,1
171 
172 %HMMPLOT - plots a Hidden Markov Model.
173 t=1:400;
174 hmmplot(x(t),z(t),t,[1 2])         % Same colour
175 
176 pause
177 
178 hmmplot(x(t),z(t),t,[1 2],'','',1) % Different colours
179 
180 pause
181 
182 %SMCTP2RFM - Calculates the rainflow matrix for a SMCTP.
183 Grfc=smctp2rfm(P,{G1,[];G2,[]});
184 Grfc1=mctp2rfm({G1,[]});
185 Grfc2=mctp2rfm({G2,[]});
186 GrfcSum=statP(1)*Grfc1+statP(2)*Grfc2;
187 
188 cmatplot(u,u,{Grfc1 Grfc2; GrfcSum Grfc})
189 
190 pause
191 
192 %SMCTP2STAT - Stationary distributions for a switching MCTP.
193 [ro,ro_min,ro_max] = smctp2stat(P,{G1,[];G2,[]});
194 [ro,ro_min,ro_max,Ro_min,Ro_max,QQ] = smctp2stat(P,{G1,[];G2,[]});
195 
196 subplot(2,1,1), 
197 plot(u,ro_min{1},u,ro_max{1}), hold on
198 plot(u,ro_min{2},u,ro_max{2}), hold off
199 uu = u(floor(1:0.5:n+0.5));
200 subplot(2,1,2), plot(uu,Ro_min,uu,Ro_max)
201 
202 pause
203 
204 %%%
205 % Decomposition of a Mixed rainflow matrix
206 % Model and Estimation
207 
208 n = 8; param = [-1 1 n];
209 M1.x0=[-0.4 -0.3]; M1.s=0.15; M1.lam=1; 
210 M2.x0=[0.3 0.4]; M2.s=0.15; M2.lam=1;
211 G1 = mktestmat(param,M1.x0,M1.s,M1.lam);
212 G2 = mktestmat(param,M2.x0,M2.s,M2.lam);
213 P=[1-0.1 0.1; 0.05 1-0.05]                % Transition matrix
214 [xD,z] = smctpsim(P,{G1 []; G2 []},5000); % Simulate
215 Fobs = dtp2rfm(xD,n);                     % observed mixed rainflow matrix
216 
217 %ESTSMCTP - Estimate SMCTP model from an observed rainflow matrix.
218 
219 known1.F = {G1 []; G2 []};   % known min-max and max-min matrices
220 init1.P = P;                 % initial guess of P-matrix
221 [Fest1,Est1] = estsmctp(Fobs,'P','ML',known1,[],init1);
222 Est1.P     % Estimated P-matrix
223 
224 known3.Ffun = 'f_funm';      % Function for calculating a submodel
225 known3.trModel2X = 'tr_m2x'; % transform from Model to X-vector
226 known3.trX2Model = 'tr_x2m'; % transform from X-vector to model
227 known3.param =param;
228 init3.P = P;       % initial guess of P-matrix
229 init3.M = {M1 M2}; % initial guess of Models for min-max mat
230 [Fest3,Est3] = estsmctp(Fobs,'P,CalcF','ML',known3,[],init3);
231 Est3.P     % Estimated P-matrix
232 Est3.M{:}  % Estimated parameters in models
233 
234 Pest_z = estmc(z,2)
235 
236 mc2stat(Est1.P)
237 mc2stat(Est3.P)
238 mc2stat(Pest_z)
239 mc2stat(P)
240 
241 pause
242 
243 % Methods for Estimation (Optional)
244 
245 known.F = {G1 []; G2 []};   % known min-max and max-min matrices
246 init.P = P;                 % initial guess of P-matrix
247 [Fest,Est] = estsmctp(Fobs,'P','ML',known,[],init); Est.P
248 [Fest,Est] = estsmctp(Fobs,'P','chi2',known,[],init); Est.P
249 [Fest,Est] = estsmctp(Fobs,'P','HD',known,[],init); Est.P
250 [Fest,Est] = estsmctp(Fobs,'P','KL',known,[],init); Est.P
251 
252 %%%%%
253 % RFM2DTP  Reconstructs a sequence of turning points from a rainflow matrix.
254 
255 x=load('sea.dat');
256 param = [-2 2 64]; n=param(3);
257 dtp0 = dat2dtp(param,x(:,2));
258 [RFM,RFM0,res0] = dtp2rfm(dtp0,n);
259 dtp = rfm2dtp(RFM0,res0);
260 plot(1:length(dtp0),dtp0,'b',1:length(dtp),dtp,'r')
261 
262 echo off
263 
264 
265 
266

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