Home > wafo > wsim > smctpsim.m

# smctpsim

## PURPOSE

Simulates a switching Markov chain of turning points,

## SYNOPSIS

[x,z,TT] = smctpsim(P,F,T,init,whatOut)

## DESCRIPTION

``` SMCTPSIM  Simulates a switching Markov chain of turning points,
i.e. a switching process with a Markov chain of turning points
within each regime.
The switching process x has the state space {1,2,...,n} and
the regime process z has the state space {1,2,...,r}.

[x,z] = smctpsim(P,F,T);
[x,z] = smctpsim(P,F,T,init);
[x,z] = smctpsim(P,F,T,init,'x');
[RFM,RFM0,res] = smctpsim(P,F,T,init,'RFM');
[x,z,RFM] = smctpsim(P,F,T,init,'x,RFM');

x       = Simulated switching Markov turning points.
z       = Simulated regime process.
RFM     = Rainflow matrix for x.                        [n,n]
RFM0    = Rainflow matrix for x (without the residual). [n,n]
res     = Residual from rainflow count.                 [n,2]

P       = Transition matrix for regime process.      [r,r]
F       = Cell array of min-Max and Max-min matrices {r,2}
F{i,1}  = min-Max matrix, process i                  [n,n]
F{i,2}  = Max-min matrix, process i                  [n,n]
T       = Length of simulation.
init.x0 = Initial state of process x. If not given, it will start from
the stationary distribution of minima given z(1).
init.z0 = Initial state of regime process. If not given, it will start
from the stationary distribution of the Markov chain.

If a matrix F{i,2}=[], then the process will
be assumed to be time-reversible.

One can also simulate a switching process other than Markov switching.
In that case P is a struct array describing the switching
P.P        = Embedded transition matrix
P.distr    = 'exp':       Exponential distribution (Markov switching)
'phasetype': Phasetype distribution
'const':     Constant
'other':     Other distribution (simTfun)
P.simTfun  = Function for simulating occupation times.
P.Tpar     = Parameters to distributions       {r,1}
P.sequence = Specified sequence.               [qx2]

Examples: Two regime states.
[x,z] = smctpsim(P,F,T); % Two regime states
Example: A Markov chain of turning points (One regime state).
[x,z] = smctpsim(1,F,T);```

## CROSS-REFERENCE INFORMATION

This function calls:
 dtp2rfm Calculates rainflow matrix from discrete turning points. mc2stat Calculates the stationary distribution for a Markov chain. mcsim Simulates a Markov chain. mctp2stat Calculates the stationary distribution for a MCTP. smctp2joint Calculates the joint MCTP for a SMCTP. clock Current date and time as date vector. error Display message and abort function. isa True if object is a given class. iscell True for cell array. isfield True if field is in structure array. strcmp Compare strings. tril Extract lower triangular part. triu Extract upper triangular part. warning Display warning message; disable or enable warning messages.
This function is called by:
 itmkurs_lab2 Script to computer exercises 2 itmkurs_lab4 Script to computer exercises 4 mctpsim Simulates a Markov chain of turning points rfcdemo2 Rainflow matrix for Switching Markov Chains of Turning Points. test_markov Quick test of the routines in module 'markov'

## SOURCE CODE

```0001 function [x,z,TT] = smctpsim(P,F,T,init,whatOut)
0002
0003 %SMCTPSIM  Simulates a switching Markov chain of turning points,
0004 %   i.e. a switching process with a Markov chain of turning points
0005 %   within each regime.
0006 %   The switching process x has the state space {1,2,...,n} and
0007 %   the regime process z has the state space {1,2,...,r}.
0008 %
0009 % [x,z] = smctpsim(P,F,T);
0010 % [x,z] = smctpsim(P,F,T,init);
0011 % [x,z] = smctpsim(P,F,T,init,'x');
0012 % [RFM,RFM0,res] = smctpsim(P,F,T,init,'RFM');
0013 % [x,z,RFM] = smctpsim(P,F,T,init,'x,RFM');
0014 %
0015 % x       = Simulated switching Markov turning points.
0016 % z       = Simulated regime process.
0017 % RFM     = Rainflow matrix for x.                        [n,n]
0018 % RFM0    = Rainflow matrix for x (without the residual). [n,n]
0019 % res     = Residual from rainflow count.                 [n,2]
0020 %
0021 % P       = Transition matrix for regime process.      [r,r]
0022 % F       = Cell array of min-Max and Max-min matrices {r,2}
0023 % F{i,1}  = min-Max matrix, process i                  [n,n]
0024 % F{i,2}  = Max-min matrix, process i                  [n,n]
0025 % T       = Length of simulation.
0026 % init.x0 = Initial state of process x. If not given, it will start from
0027 %          the stationary distribution of minima given z(1).
0028 % init.z0 = Initial state of regime process. If not given, it will start
0029 %          from the stationary distribution of the Markov chain.
0030 %
0031 % If a matrix F{i,2}=[], then the process will
0032 % be assumed to be time-reversible.
0033 %
0034 % One can also simulate a switching process other than Markov switching.
0035 % In that case P is a struct array describing the switching
0036 % P.P        = Embedded transition matrix
0037 % P.distr    = 'exp':       Exponential distribution (Markov switching)
0038 %              'phasetype': Phasetype distribution
0039 %              'const':     Constant
0040 %              'other':     Other distribution (simTfun)
0041 % P.simTfun  = Function for simulating occupation times.
0042 % P.Tpar     = Parameters to distributions       {r,1}
0043 % P.sequence = Specified sequence.               [qx2]
0044 %
0045 % Examples: Two regime states.
0046 %   [x,z] = smctpsim(P,F,T); % Two regime states
0047 % Example: A Markov chain of turning points (One regime state).
0048 %   [x,z] = smctpsim(1,F,T);
0049
0050 % Tested  on Matlab  5.3
0051 %
0052 % History:
0053 % Correction by PJ 20-Jan-2004
0054 %   Some corrections concerning initial state of regime process
0055 %   and simulation of non-Markovian regime process.
0056 % Correction by PJ 12-Aug-2002
0057 %   'RFM' opt didn't work cause x was not initiated. Now fixed!
0058 % Revised by PJ 19-May-2000
0059 %   Corrected method for simulating starting conditions.
0060 % Revised by PJ Jan-2000
0061 %   updated for WAFO
0062 % Created by PJ (Pär Johannesson) 1997
0063 %   Copyright (c) 1997 by Pär Johannesson
0064 %   Toolbox: Rainflow Cycles for Switching Processes V.1.0, 2-Oct-1997
0065
0066 TT(1,:) = clock;
0067
0068 ni = nargin;
0069 no = nargout;
0070 error(nargchk(3,5,ni));
0071
0072 Zstr = '123456789';
0073
0074 if ni < 4,  init = []; end
0075 if ni < 5,  whatOut = []; end
0076
0077 if isempty(init)
0078     init.x0 = [];
0079     init.z0 = [];
0080 end
0081 if ~isfield(init,'x0'), init.x0=[]; end
0082 if ~isfield(init,'z0'), init.z0=[]; end
0083
0084 if isempty(whatOut)
0085     whatOut = 'x';
0086 end
0087
0088 if isa(P,'double')
0089     Ptype = 'P';
0090 elseif isa(P,'struct')
0091     Ptype = 'struct';
0092     S = P;
0093 else
0094     error('P should be matrix or struct-array.');
0095 end
0096
0097 r = size(F,1);   % Number of regime states
0098 %r = length(P);   % Number of regime states
0099
0100 n = length(F{1,1});  % Number of levels
0101
0102 % Check that the rowsums of P are equal to 1
0103
0104 if strcmp(Ptype,'struct')
0105     if isfield(S,'P')
0106         P=S.P;
0107     else
0108         P=[];
0109     end
0110 end
0111
0112 if ~isempty(P)
0113     sumP = sum(P');
0114     if sum(sumP == 1) ~= length(P)
0115         warning([': Rowsums of P not equal to 1. Renormalizing!']);
0116         for i = 1:length(P)
0117             P(i,:) = P(i,:)/sumP(i);
0118         end
0119     end
0120 end
0121
0122 TT(2,:) = clock;
0123
0124 % Normalize the rowsums of F{1,1},...,F{r,1} to 1
0125 %  ==>  Q{1,1},...,Q{r,1}
0126
0127 for i = 1:r
0128     QQ{i,1} = triu(F{i,1},1); % Set zeros below diagonal and diagonal
0129     % Set negative elements to zero
0130     [I,J] = find(QQ{i,1}<0);
0131     if length(I) ~= 0
0132         warning(['Negative elements in Q' Zstr(i) '. Setting to zero!']);
0133         for k = 1:length(I)
0134             QQ{i,1}(I(k),J(k)) = 0;
0135         end
0136     end
0137
0138     sumQQi = sum(QQ{i,1}');
0139     % Normalize rowsums
0140     if sum(sumQQi == 1) ~= length(QQ{i,1})
0141         %disp(['Warning: Rowsums of Q' Zstr(i) ' not equal to 1. Renormalizing!']);
0142         for j = 1:n-1
0143             if sumQQi(j)~=0, QQ{i,1}(j,:) = QQ{i,1}(j,:)/sumQQi(j); end
0144         end
0145     end
0146 end
0147
0148 TT(3,:) = clock;
0149
0150 % Normalize the rowsums of F{1,2},...,F{r,2} to 1
0151 %  ==>  Q{1,2},...,Q{r,2}
0152
0153 % Normalize the rowsums of Fh1,...,Fhr to 1  ==>  Q1,...,Qr
0154
0155 for i = 1:r
0156     if isempty(F{i,2})        % Time-reversible
0157         QQ{i,2} = F{i,1}';
0158     else
0159         QQ{i,2} = F{i,2};
0160     end
0161
0162     QQ{i,2} = tril(QQ{i,2},-1); % Set zeros above diagonal and diagonal
0163     % Set negative elements to zero
0164     [I,J] = find(QQ{i,2}<0);
0165     if length(I) ~= 0
0166         warning(['Negative elements in Qh' Zstr(i) '. Setting to zero!']);
0167         for k = 1:length(I)
0168             QQ{i,2}(I(k),J(k)) = 0;
0169         end
0170     end
0171
0172     sumQQi = sum(QQ{i,2}');
0173     if sum(sumQQi == 1) ~= length(QQ{i,2})
0174         %disp(['Warning: Rowsums of Qh' Zstr(i) ' not equal to 1. Renormalizing!']);
0175         for j = 2:n
0176             if sumQQi(j)~=0, QQ{i,2}(j,:) = QQ{i,2}(j,:)/sumQQi(j); end
0177         end
0178     end
0179 end
0180
0181 TT(4,:) = clock;
0182
0183
0184
0185 % Initial state of z0, for regime process
0186 % and x0, for X-process, start from a minimum
0187
0188 % Make the transition matrix Q for the joint MC (X_k,Z_k)
0189 [Q] = smctp2joint(P,QQ);
0190
0191 % Stationary distribution = ro of minima
0192 [ro_min,ro_max] = mctp2stat(Q);
0193 ro = ro_min;
0194
0195 % Start values
0196 e0 = rand(1,1);  % Generate random numbers
0197 if isempty(init.z0) & isempty(init.x0)
0198     x0z0 = min(find( e0<=cumsum(ro) ));
0199     x0 = floor((x0z0+1)/r);
0200     z0 = mod(x0z0-1,r)+1;
0201 elseif isempty(init.x0)
0202     z0 = init.z0;
0203     rox0 = ro(z0:r:end); % Pick stat. distr. for regime z0
0204     rox0 = rox0/sum(rox0);
0205     x0 = min(find( e0<=cumsum(rox0) ));
0206 elseif isempty(init.z0)
0207     x0 = init.x0;
0208     z0 = [];  % Start from stat. distr of P
0209 else % Both z0 znd x0 are given
0210     x0 = init.x0;
0211     z0 = init.z0;
0212 end
0213
0214
0215 if strcmp(Ptype,'struct')
0216     if isfield(S,'P')
0217         S.P=P;
0218     end
0219 end
0220
0221
0222 %
0223 % Initiate vectors
0224 %
0225
0226 switch whatOut
0227
0228     case {'x','x,RFM'}
0229         if strcmp(Ptype,'P')
0230             z = mcsim(P,T,z0);  % Simulate Regime
0231         else
0232             z=zeros(T,1);
0233             z(1)=init_z(S,init);
0234         end
0235         e=rand(T,1);
0236         x=zeros(T,1);
0237
0238     case {'RFM'}
0239         if strcmp(Ptype,'P')
0240             z = mcsim(P,1,z0);  % Simulate Regime
0241         else
0242             z=init_z(S,init);
0243         end
0244         e=rand(1,1);
0245
0246 end
0247
0248 x(1) = x0;
0249
0250 TT(5,:) = clock;
0251
0252 %
0253 % Simulate Switching Markov turning points
0254 %
0255
0256 % Calculate cumulative distributions
0257
0258 cumsumP = cumsum(P,2);
0259 ones_n = ones(n,1);
0260 for i = 1:r
0261     cumsumQQ{i,1} = cumsum(QQ{i,1},2);
0262     cumsumQQ{i,2} = cumsum(QQ{i,2},2);
0263     cumsumQQ{i,1}(:,end) = ones_n;
0264     cumsumQQ{i,2}(:,end) = ones_n;
0265 end
0266
0267 % Simulate
0268
0269 switch whatOut
0270
0271     case {'x','x,RFM'}
0272         switch Ptype
0273             case 'P'
0274                 for k=2:T
0275                     if rem(k,2) == 0 % min-to-Max
0276                         x(k) = sum( cumsumQQ{z(k),1}(x(k-1),:) <= e(k) ) + 1;
0277                     else             % Max-to-min
0278                         x(k) = sum( cumsumQQ{z(k),2}(x(k-1),:) <= e(k) ) + 1;
0279                     end
0280                 end
0281
0282             case 'struct'
0283                 t0 = simT(S,z(1));
0284                 for k=2:T
0285                     if t0 > 1
0286                         z(k) = z(k-1);
0287                         t0=t0-1;
0288                     else
0289                         z(k-1:k) = mcsim(S.P,2,z(k-1));
0290                         t0 = simT(S,z(k));
0291                     end
0292                     %      fprintf(1,'k=%d, z=%d\n',k,z(k));
0293                     if rem(k,2) == 0 % min-to-Max
0294                         x(k) = sum( cumsumQQ{z(k),1}(x(k-1),:) <= e(k) ) + 1;
0295                     else             % Max-to-min
0296                         x(k) = sum( cumsumQQ{z(k),2}(x(k-1),:) <= e(k) ) + 1;
0297                     end
0298                 end
0299         end
0300
0301
0302     case 'test' %'x,RFM'
0303         [RFM0,res,nres] = dtp2rfm_init(n);
0304         [RFM0,res,nres] = dtp2rfm1(x(1),RFM0,res,nres);
0305         for k=2:T
0306             e=rand(2,1);
0307
0308             % Simulate Regime
0309             z(k) = sum( cumsumP(z(k-1),:) <= e(1) ) + 1;
0310
0311             % Simulate MCTP
0312             if rem(k,2) == 0 % min-to-Max
0313                 x(k) = sum( cumsumQQ{z(k),1}(x(k-1),:) <= e(2) ) + 1;
0314             else             % Max-to-min
0315                 x(k) = sum( cumsumQQ{z(k),2}(x(k-1),:) <= e(2) ) + 1;
0316             end
0317             [RFM0,res,nres] = dtp2rfm1(x(k),RFM0,res,nres);
0318         end
0319         [RFM] = dtp2rfm_collect(RFM0,res,nres);
0320
0321     case 'RFM'
0322         [RFM0,res,nres] = dtp2rfm_init(n);
0323         [RFM0,res,nres] = dtp2rfm1(x,RFM0,res,nres);
0324         for k=2:T
0325             e=rand(2,1);
0326
0327             % Simulate Regime
0328             z = sum( cumsumP(z,:) <= e(1) ) + 1;
0329
0330             % Simulate MCTP
0331             if rem(k,2) == 0 % min-to-Max
0332                 x = sum( cumsumQQ{z,1}(x,:) <= e(2) ) + 1;
0333             else             % Max-to-min
0334                 x = sum( cumsumQQ{z,2}(x,:) <= e(2) ) + 1;
0335             end
0336             [RFM0,res,nres] = dtp2rfm1(x,RFM0,res,nres);
0337         end
0338         [RFM] = dtp2rfm_collect(RFM0,res,nres);
0339
0340 end
0341
0342 TT(6,:) = clock;
0343
0344 % Output arguments
0345
0346 switch whatOut
0347     case 'x,RFM'
0348         RFM = dtp2rfm(x,n);
0349         TT = RFM;
0350     case 'RFM'
0351         x = RFM;
0352 end
0353
0354
0355 function z=init_z(S,init)
0356
0357 % Initiate regime.
0358
0359 if isempty(init.z0)
0360     r=length(S.P);
0361     m = zeros(1,r);
0362     switch S.distr
0363         case {'exp','const'},
0364             for i=1:r, m(i)=S.Tpar{i}; end
0365         case 'phasetype',
0366             for i=1:r, m(i)=sum(S.Tpar{i}); end
0367         end
0368         statP=mc2stat(S.P);
0369         ro = statP.*m;
0370         z = sum( cumsum(ro) <= rand ) + 1;
0371     else
0372         z=init.z0;
0373     end
0374
0375
0376     function t=simT(S,z)
0377
0378     % Simulate occupation time.
0379
0380     if iscell(S.distr)
0381         distr = S.distr{z};
0382     else
0383         distr = S.distr;
0384     end
0385
0386     switch distr
0387
0388         case 'exp'
0389             t = simFS(S.Tpar{z});
0390         case 'phasetype'
0391             t = simPhaseType(S.Tpar{z});
0392         case 'const'
0393             t = S.Tpar{z};
0394         case 'other'
0395             error('other: Not yet implemented.');
0396     end
0397
0398     % Simulate First Success R.V.
0399
0400     function t=simFS(m)
0401     % Simulation method taken from stats-toolbox 'geornd'
0402
0403     u = rand(size(m));
0404     p = 1./m;
0405     t = floor(log(u) ./ log(1 - p)) + 1;
0406     %t = - m .* log(rand(size(m)));
0407
0408     % Simulate Exponential R.V.
0409     function t=simExp(m)
0410
0411     t = - m .* log(rand(size(m)));
0412
0413     % Simulate discrete Phasetype R.V.
0414
0415     function t=simPhaseType(m)
0416
0417     t = sum(simFS(m));
0418     %t = sum(simExp(m));
0419
0420     % Description of variables for  dtp2rfm
0421     %
0422     % RFM   = Rainflow Matrix (residual included).    [nxn]
0423     % RFM0  = Rainflow matrix (without resudual).     [nxn]
0424     % res   = Residual.                               [2*n,1]
0425     % nres  = Length of residual
0426     % x     = Turning points (taking values 1,...,n). [T,1]
0427     % n     = Number of levels.
0428
0429     function [RFM0,res,nres] = dtp2rfm_init(n)
0430
0431     % Initiate variables RFM0,res,nres
0432
0433     RFM0 = zeros(n);
0434     res = zeros(2*n+1,1);
0435     nres = 0;
0436
0437
0438     function [RFM0,res,nres] = dtp2rfm1(x,RFM0,res,nres)
0439
0440     % Count one TP.
0441
0442
0443     % Calculate RFM and res
0444
0445     %for i = 1:length(x)
0446     nres = nres+1;
0447     res(nres) = x; %(i);
0448     cycleFound = 1;
0449     while cycleFound==1 & nres>=4
0450         A = sort([res(nres-1) res(nres-2)]);
0451         B = sort([res(nres) res(nres-3)]);
0452         if A(1) >= B(1) & A(2) <= B(2)
0453             RFM0(res(nres-2),res(nres-1)) = RFM0(res(nres-2),res(nres-1)) + 1;
0454             res(nres-2) = res(nres);
0455             nres = nres-2;
0456         else
0457             cycleFound = 0;
0458         end
0459     end
0460     %end
0461
0462
0463
0464     function [RFM] = dtp2rfm_collect(RFM0,res,nres)
0465
0466     % Collect RFM0 and residuals. Store in RFM.
0467
0468
0469     % Calculate RFM = RFM0 + 'cycles in res'
0470
0471     RFM = RFM0;
0472     for i=1:2:nres-1
0473         RFM(res(i),res(i+1)) = RFM(res(i),res(i+1)) + 1;
0474     end
0475
0476     % Convert to symetric rainflow
0477
0478     RFM = RFM+RFM';
0479     RFM = triu(RFM);
0480
0481```

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