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: This function is called by:

SUBFUNCTIONS ^

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