Home > wafo > markov > smc2rfm.m

# smc2rfm

## PURPOSE

Calculates the rainflow matrix/intensity for a switching Markov chain.

## SYNOPSIS

[F_rfc,mu_rfc] = smc2rfc(P,Qc,def)

## DESCRIPTION

``` SMC2RFM Calculates the rainflow matrix/intensity for a switching Markov chain.

[F_rfc,mu_rfc] = smc2rfc(P,Q,1);
[F_rfc,mu_rfc] = smc2rfc(P,Q,[2 h]);

F_rfc  = Rainflow matrix / Rainflow intensity     [NxN]
mu_rfc = Rainflow counting intensity              [NxN]

P      = Transition matrix for regime process     [rxr]
Q      = Cell array of transition matrices        {r,1}
Q{i}   = Transition matrix for Markov chain i     [nxn]
def    = Definition 1: Markov chain (default),     N=n
2: Discretized Markov chain,   N=n+1
h      = Discretization step (ONLY Def 2!)

Calculates
(1) the rainflow matrix for a switching Markov chain OR
(2) the rainflow intensity for a discretized switching Markov chain.

Example:
param = [-1 1 32]; u = levels(param);
[F1,F2] = mktestmat(param,[-0.3 0.3],0.15,1,-Inf);
Frfc = smc2rfm(P,{F1;F2});
cmatplot(u,u,Frfc)

## CROSS-REFERENCE INFORMATION

This function calls:
 cmat2nt Calculates a counting distribution from a cycle matrix. mc2stat Calculates the stationary distribution for a Markov chain. error Display message and abort function. inv Matrix inverse. kron Kronecker tensor product. warning Display warning message; disable or enable warning messages.
This function is called by:
 mc2rfm Calculates the rainflow matrix/intensity for a Markov chain. rfcdemo1 Demo for switching AR(1)-processes. test_markov Quick test of the routines in module 'markov'

## SOURCE CODE

```001 function [F_rfc,mu_rfc] = smc2rfc(P,Qc,def)
002 %SMC2RFM Calculates the rainflow matrix/intensity for a switching Markov chain.
003 %
004 % [F_rfc,mu_rfc] = smc2rfc(P,Q,1);
005 % [F_rfc,mu_rfc] = smc2rfc(P,Q,[2 h]);
006 %
007 % F_rfc  = Rainflow matrix / Rainflow intensity     [NxN]
008 % mu_rfc = Rainflow counting intensity              [NxN]
009 %
010 % P      = Transition matrix for regime process     [rxr]
011 % Q      = Cell array of transition matrices        {r,1}
012 % Q{i}   = Transition matrix for Markov chain i     [nxn]
013 % def    = Definition 1: Markov chain (default),     N=n
014 %                     2: Discretized Markov chain,   N=n+1
015 % h      = Discretization step (ONLY Def 2!)
016 %
017 % Calculates
018 %   (1) the rainflow matrix for a switching Markov chain OR
019 %   (2) the rainflow intensity for a discretized switching Markov chain.
020 %
021 % Example:
022 %   param = [-1 1 32]; u = levels(param);
023 %   [F1,F2] = mktestmat(param,[-0.3 0.3],0.15,1,-Inf);
024 %   Frfc = smc2rfm(P,{F1;F2});
025 %   cmatplot(u,u,Frfc)
026 %
028
029 % References
030 %
031 %  P. Johannesson (1999):
032 %  Rainflow Analysis of Switching Markov Loads.
033 %  PhD thesis, Mathematical Statistics, Centre for Mathematical Sciences,
034 %  Lund Institute of Technology.
035 %
036 %  P. Johannesson (1998):
037 %  Rainflow Cycles for Switching Processes with Markov Structure.
038 %  Probability in the Engineering and Informational Sciences,
039 %  Vol. 12, No. 2, pp. 143-175.
040
041 % Tested  on Matlab  5.3
042 %
043 % History:
044 % Revised by PJ 19-May-2000
045 %   Changer disp(...) to warning(...) .
046 % Revised by PJ  23-Nov-1999
047 %   updated for WAFO
048 % Created by PJ (Pär Johannesson) 1997
049 %   Copyright (c) 1997 by Pär Johannesson
050 %   Toolbox: Rainflow Cycles for Switching Processes V.1.0, 2-Oct-1997
051
052
053 % Check input arguments
054
055 ni = nargin;
056 no = nargout;
057 error(nargchk(2,3,ni));
058
059 if ni<3, def = []; end
060 if isempty(def), def = 1; end
061
062 % Define
063
064 Zstr = '123456789';
065
066 r = length(P);   % Number of regime states
067 n = length(Qc{1});  % Number of levels
068
069 % Check that the rowsums of P are equal to 1
070
071 sumP = sum(P');
072 if sum(sumP == 1) ~= length(P)
073   warning(['Rowsums of P not equal to 1. Renormalizing.']);
074   for i = 1:length(P)
075     P(i,:) = P(i,:)/sumP(i);
076   end
077 end
078
079 % Check that the rowsums of Qc{1},...,Qc{r} are equal to 1
080
081 for i = 1:r
082   sumQi = sum(Qc{i}');
083   if sum(sumQi == 1) ~= length(Qc{i})
084     warning(['Rowsums of Q{' Zstr(i) '} not equal to 1. Renormalizing.']);
085     for j = 1:length(Qc{i})
086       Qc{i}(j,:) = Qc{i}(j,:)/sumQi(j);
087     end
088   end
089 end
090
091
092 % Make the transition matrix Q for the joint MC (X_k,Z_k)
093
094 Q = zeros(n*r,n*r);
095 I = 0:r:(n-1)*r;
096 for z=1:r
097   QQ = kron(Qc{z},P);
098   Q(I+z,:) = QQ(I+z,:);
099 end
100
101
102 % Stationary distribution = ro of Q
103
104 ro = mc2stat(Q);
105
106 % Calculate rainflow matrix / rainflow intensity
107
108 if def(1) == 1  % Markov Chain
109
110   N = n;
111   mu_rfc = zeros(N,N);
112   EYE = eye(n*r,n*r);
113   Qcumsum = fliplr(cumsum(fliplr(Q)')');
114
115   for i=2:n
116     for j=i-1:n-1
117       q = Qcumsum(:,r*j+1);
118
119       A = Q(r*(i-1)+1:r*j,r*(i-1)+1:r*j); % i:j,   i:j
120       C = Q(1:r*(i-1),r*(i-1)+1:r*j);     % 1:i-1, i:j
121       Eye = EYE(1:r*(j-i+1),1:r*(j-i+1)); % eye (size of A)
122       d = q(1:r*(i-1));                   % 1:i-1
123       e = q(r*(i-1)+1:r*j);               % i:j
124       Ro = ro(1:r*(i-1));                 % 1:i-1
125
126       if j == i-1  % isempty(A)
127         mu_rfc(i,j) = Ro*d;
128       else
129         mu_rfc(i,j) = Ro*(d+C*inv(Eye-A)*e);
130       end
131     end
132   end
133
134 elseif def(1) == 2  % Discretized Markov Process
135
136   N = n+1;
137   mu_rfc = zeros(N,N);
138   EYE = eye(n*r,n*r);
139   Qcumsum = fliplr(cumsum(fliplr(Q)')');
140
141   for i=2:N-1
142     for j=i:N-1
143       q = Qcumsum(:,r*(j-1)+1);
144
145       A = Q(r*(i-1)+1:r*(j-1),r*(i-1)+1:r*(j-1)); % i:j-1, i:j-1
146       C = Q(1:r*(i-1),r*(i-1)+1:r*(j-1));         % 1:i-1, i:j-1
147       Eye = EYE(1:r*(j-i),1:r*(j-i));             % eye (size of A)
148       d = q(1:r*(i-1));                           % 1:i-1
149       e = q(r*(i-1)+1:r*(j-1));                   % i-1:j-1
150       Ro = ro(1:r*(i-1));                         % 1:i-1
151
152       if i == j  % isempty(A)
153         mu_rfc(i,j) = Ro*d;
154       else
155         mu_rfc(i,j) = Ro*(d+C*inv(Eye-A)*e);
156       end
157     end
158   end
159
160 end
161
162 % Fill the missing triangel in the mu_rfc matrix
163 if def(1) == 1
164
165   F_rfc = zeros(n);
166   for i = 1:n-1
167     for j= i+1:n
168       F_rfc(i,j) = mu_rfc(i+1,j-1)-mu_rfc(i,j-1)-mu_rfc(i+1,j)+mu_rfc(i,j);
169     end
170   end
171
172   mu_rfc = cmat2nt(F_rfc);
173
174 % Old version
175 %  mu_rfc=flipud(mu_rfc'); % Convert to WAT matrix-format
176 %    for k = 2:N-1
177 %      for m = 1:N-k
178 %        i = k+m;
179 %        j = N+1+k-i;
180 %        mu_rfc(i,j) = mu_rfc(i,j-1)+mu_rfc(i-1,j)-mu_rfc(i-1,j-1);
181 %      end
182 %    end
183 %
184 %  mu_rfc = flipud(mu_rfc)'; % Convert to PJ-def
185 %
186 %  F_rfc = nt2cmat(mu_rfc);    % Calculate Rainflow matrix
187 %
188 %%  F_rfc = fliplr(triu(fliplr(F_rfc),1)); % Zeros blow SW-NE diagonal
189
190 end
191
192 if def(1) == 2
193
194   mu_rfc=flipud(mu_rfc'); % Convert to WAT matrix-format
195   for k = 1:N-1
196     for m = 1:N-k
197       i = k+m;
198       j = N+1+k-i;
199       mu_rfc(i,j) = mu_rfc(i,j-1)+mu_rfc(i-1,j)-mu_rfc(i-1,j-1);
200     end
201   end
202
203   % Calculate rainflow intensity by numerical differentiation
204
205   h = def(2);
206   F_rfc = zeros(N,N);
207
208   for i=2:N-1
209     for j=2:N-i
210       F_rfc(i,j) = (mu_rfc(i-1,j-1)-mu_rfc(i-1,j+1) -(mu_rfc(i+1,j-1)-mu_rfc(i+1,j+1))) / (2*h)^2;
211     end
212   end
213
214   F_rfc = flipud(F_rfc)';   % Convert to PJ-def
215   mu_rfc = flipud(mu_rfc)'; % Convert to PJ-def
216
217 end
218
219
220
221```

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