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.
Plots a cycle matrix, e.g. a rainflow matrix. | |
The sequence of discretized turning points from a signal. | |
Calculates asymmetric RFM from discrete turning points. | |
Calculates rainflow matrix from discrete turning points. | |
estimates transition matrix P from a time series of a Markov chain. | |
Estimate SMCTP model from an observed rainflow matrix. | |
plots a Hidden Markov Model. | |
Calculates discrete levels given the parameter matrix. | |
Converts a matrix to a transition matrix. | |
Calculates the rainflow matrix/intensity for a Markov chain. | |
Calculates the stationary distribution for a Markov chain. | |
Simulates a Markov chain. | |
Calculates asymmetric rainflow matrix for a MCTP. | |
Calculates the rainflow matrix for a MCTP. | |
Calculates the stationary distribution for a MCTP. | |
Simulates a Markov chain of turning points | |
Makes test matrices for min-max (and max-min) matrices. | |
Rainflow filter a signal. | |
Reconstructs a sequence of turning points from a rainflow matrix. | |
Calculates the rainflow matrix/intensity for a switching Markov chain. | |
Simulates a Switching Markov chain with state space. | |
Calculates the rainflow matrix for a SMCTP. | |
Stationary distributions for a switching MCTP. | |
Simulates a switching Markov chain of turning points, |
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
Comments or corrections to the WAFO group