ESTSMCTP Estimate SMCTP model from an observed rainflow matrix. Estimates parameters in a Switching Markov Chain of Turning Points from an observation of the rainflow matrix. CALL: [Fest,Est] = estsmtp(Fobs,whatEst,method,known,whatKnown,init,OPTIONS) Fest = Estimated SMCTP model. [SA] Est = Estimated parameters. [SA] Fobs = Observation of rainflow matrix. [nxn] whatEst = See below. method = 'ML' / 'chi2' / 'HD' / 'KL' (See below.) known = Values of known parameters of the model. [SA] whatKnown = Which parameters are known? (Not used!) [SA] init = Initial guess. (for optimization) [SA] OPTIONS = Options to optimization routine. (Optional) (see 'help foptions') [SA]=[structure array] method: 'ML' : Approximate Maximum Likelihood (Multinomial) 'chi2' : Chi-square distance 'HD' : Hellinger distance 'KL' : Kullback-Leibler distance whatEst: 'P' : Estimate P-matrix, min-max matrices for the subloads are known [known.F]. 'MeanStd' : Estimate mean and std of subloads. The shape of min-max matrices for the subloads are known [known.F]. P-matrix known [known.P]. 'P,MeanStd' : Also estimate P-matrix otherwise as above. 'CalcF' : 'P,CalcF' : 'SimF' : 'P,SimF' : Side Information: known.SideInfo = 0: No side information 11: Mark min&max, y = 'regime process' 12: Mark min&max, y = 'scrambled regime process' 21: Mark when counted, y = 'regime process' 22: Mark when counted, y = 'scrambled regime process' (Optional, Default = 0, No side information) known.NOsubzero = Number of subdiagonals that are zero (Optional, Default = 0, only the diagonal is zero) Example: M1.x0=[-0.4 -0.3]; M1.s=0.15; M1.lam=1; M2.x0=[0.3 0.4]; M2.s=0.15; M2.lam=1; F1 = mktestmat([-1 1 8],M1.x0,M1.s,M1.lam); F2 = mktestmat([-1 1 8],M2.x0,M2.s,M2.lam); P=[1-0.1 0.1; 0.05 1-0.05] % Transition matrix [xD,z] = smctpsim(P,{F1 []; F2 []},5000); % Simulate Fobs = dtp2rfm(xD,8); known.F = {F1 []; F2 []}; % known min-max and max-min matrices init.P = P; % initial guess of P-matrix [Fest,Est] = estsmctp(Fobs,'P','ML',known,[],init); known.Ffun = 'f_funm'; % Function for calculating a submodel known.trModel2X = 'tr_m2x'; % transform from Model to X-vector known.trX2Model = 'tr_x2m'; % transform from X-vector to model known.param = [-1 1 8]; init.P = P; % initial guess of P-matrix init.M = {M1 M2}; % initial guess of Models for min-max mat [Fest,Est] = estsmctp(Fobs,'P,CalcF','ML',known,[],init); Further examples of using ESTSMCP can be found in WDEMOS/ITMKURS, especially in the script ITMKURS_LAB2.
Auxiliary function used by ESTSMCTP | |
Transform Model-structure to X-vector. | |
Transform P-matrix to X-vector | |
Transform X-vector to Model-structure. | |
Transforms a vector X to a transition matrix P. | |
Create cell array. | |
Display message and abort function. | |
Multidimensional unconstrained nonlinear minimization (Nelder-Mead). | |
True if field is in structure array. |
Script to computer exercises 2 | |
Script to computer exercises 4 | |
Quick test of the routines in module 'markov' |
001 function [Fest,Est,OPTIONS] = estsmctp(Fobs,whatEst,method,known,whatKnown,init,OPTIONS) 002 %ESTSMCTP Estimate SMCTP model from an observed rainflow matrix. 003 % 004 % Estimates parameters in a Switching Markov Chain 005 % of Turning Points from an observation of the rainflow matrix. 006 % 007 % CALL: [Fest,Est] = estsmtp(Fobs,whatEst,method,known,whatKnown,init,OPTIONS) 008 % 009 % Fest = Estimated SMCTP model. [SA] 010 % Est = Estimated parameters. [SA] 011 % 012 % Fobs = Observation of rainflow matrix. [nxn] 013 % whatEst = See below. 014 % method = 'ML' / 'chi2' / 'HD' / 'KL' (See below.) 015 % known = Values of known parameters of the model. [SA] 016 % whatKnown = Which parameters are known? (Not used!) [SA] 017 % init = Initial guess. (for optimization) [SA] 018 % OPTIONS = Options to optimization routine. (Optional) 019 % (see 'help foptions') 020 % 021 % [SA]=[structure array] 022 % 023 % method: 024 % 'ML' : Approximate Maximum Likelihood (Multinomial) 025 % 'chi2' : Chi-square distance 026 % 'HD' : Hellinger distance 027 % 'KL' : Kullback-Leibler distance 028 % 029 % whatEst: 030 % 'P' : Estimate P-matrix, min-max matrices for the 031 % subloads are known [known.F]. 032 % 'MeanStd' : Estimate mean and std of subloads. 033 % The shape of min-max matrices for the subloads are known 034 % [known.F]. P-matrix known [known.P]. 035 % 'P,MeanStd' : Also estimate P-matrix otherwise as above. 036 % 'CalcF' : 037 % 'P,CalcF' : 038 % 'SimF' : 039 % 'P,SimF' : 040 % 041 % Side Information: known.SideInfo = 042 % 0: No side information 043 % 11: Mark min&max, y = 'regime process' 044 % 12: Mark min&max, y = 'scrambled regime process' 045 % 21: Mark when counted, y = 'regime process' 046 % 22: Mark when counted, y = 'scrambled regime process' 047 % (Optional, Default = 0, No side information) 048 % 049 % known.NOsubzero = Number of subdiagonals that are zero 050 % (Optional, Default = 0, only the diagonal is zero) 051 % 052 % Example: 053 % M1.x0=[-0.4 -0.3]; M1.s=0.15; M1.lam=1; 054 % M2.x0=[0.3 0.4]; M2.s=0.15; M2.lam=1; 055 % F1 = mktestmat([-1 1 8],M1.x0,M1.s,M1.lam); 056 % F2 = mktestmat([-1 1 8],M2.x0,M2.s,M2.lam); 057 % P=[1-0.1 0.1; 0.05 1-0.05] % Transition matrix 058 % [xD,z] = smctpsim(P,{F1 []; F2 []},5000); % Simulate 059 % Fobs = dtp2rfm(xD,8); 060 % 061 % known.F = {F1 []; F2 []}; % known min-max and max-min matrices 062 % init.P = P; % initial guess of P-matrix 063 % [Fest,Est] = estsmctp(Fobs,'P','ML',known,[],init); 064 % 065 % known.Ffun = 'f_funm'; % Function for calculating a submodel 066 % known.trModel2X = 'tr_m2x'; % transform from Model to X-vector 067 % known.trX2Model = 'tr_x2m'; % transform from X-vector to model 068 % known.param = [-1 1 8]; 069 % init.P = P; % initial guess of P-matrix 070 % init.M = {M1 M2}; % initial guess of Models for min-max mat 071 % [Fest,Est] = estsmctp(Fobs,'P,CalcF','ML',known,[],init); 072 % 073 % Further examples of using ESTSMCP can be found in WDEMOS/ITMKURS, 074 % especially in the script ITMKURS_LAB2. 075 076 % Updated by PJ 07-Apr-2005 077 % Adaptation for Matlab 7. 078 % Changed 'fmins' to 'fminsearch'. 079 080 % Check input aruments 081 082 ni = nargin; 083 no = nargout; 084 error(nargchk(6,7,ni)); 085 086 if ni < 7 087 OPTIONS = []; 088 end 089 090 if ~isfield(known,'NOsubzero') 091 known.NOsubzero = 0; 092 end 093 094 if ~isfield(known,'SideInfo') 095 known.SideInfo = 0; 096 end 097 098 % Options to fmins 099 100 % OPTIONS(1) = 0; % 1 = intermediate steps in the solution are displayed 101 % OPTIONS(2) = 1e-2; % the termination tolerance for x; 102 % OPTIONS(3) = 1e-2; % the termination tolerance for F(x); 103 % OPTIONS(2) = 1e-1; % the termination tolerance for x; 104 % OPTIONS(3) = 1e-1; % the termination tolerance for F(x); 105 if 0 106 % The HTML documentation software is unable to track dependencies to 107 % functions evaluated with feval. Here is a 108 % trick to get the html documentation right (pab 2005) 109 M.s = 1; 110 X = tr_m2x(M); 111 M = tr_x2m(X,known); 112 end 113 114 switch whatEst 115 116 case 'P' % Estimate P 117 118 P=init.P; 119 X0 = tr_p2x(P,1); 120 121 try 122 % For Matlab 5.3 and higher ??? 123 [X,OPTIONS] = fminsearch('f_smctp',X0,OPTIONS,Fobs,whatEst,method,known,whatKnown,init); 124 catch 125 % For Matlab 5.2 and lower ??? 126 [X,OPTIONS] = fmins('f_smctp',X0,OPTIONS,[],Fobs,whatEst,method,known,whatKnown,init); 127 end 128 129 Pest = tr_x2p(X,1); 130 Est.P=Pest; 131 132 case 'MeanStd' % Estimate Mean and Std 133 134 init.MeanStd(:,2) = log(init.MeanStd(:,2)); 135 X0 = init.MeanStd(:); 136 137 try 138 % For Matlab 5.3 and higher ??? 139 [X,OPTIONS] = fminsearch('f_smctp',X0,OPTIONS,Fobs,whatEst,method,known,whatKnown,init); 140 catch 141 % For Matlab 5.2 and lower ??? 142 [X,OPTIONS] = fmins('f_smctp',X0,OPTIONS,[],Fobs,whatEst,method,known,whatKnown,init); 143 end 144 145 r = length(X)/2; 146 MeanStd = reshape(X,r,2); 147 MeanStd(:,2) = exp(MeanStd(:,2)); 148 Est.MeanStd = MeanStd; 149 150 case 'P,MeanStd' % Estimate P, Mean and Std 151 152 P=init.P; 153 X1 = tr_p2x(P,1); 154 init.MeanStd(:,2) = log(init.MeanStd(:,2)); 155 X2 = init.MeanStd(:); 156 X0 = [X1; X2]; 157 158 try 159 % For Matlab 5.3 and higher ??? 160 [X,OPTIONS] = fminsearch('f_smctp',X0,OPTIONS,Fobs,whatEst,method,known,whatKnown,init); 161 catch 162 % For Matlab 5.2 and lower ??? 163 [X,OPTIONS] = fmins('f_smctp',X0,OPTIONS,[],Fobs,whatEst,method,known,whatKnown,init); 164 end 165 166 r=(-1+sqrt(1+4*length(X)))/2; 167 X1 = X(1:r*(r-1)); 168 X2 = X(r*(r-1)+1:end); 169 170 Pest = tr_x2p(X1,1); 171 Est.P=Pest; 172 173 MeanStd = reshape(X2,r,2); 174 MeanStd(:,2) = exp(MeanStd(:,2)); 175 Est.MeanStd = MeanStd; 176 177 % Fest = smctp(Pest,known.F); 178 179 180 case {'CalcF','P,CalcF'} % Estimate P, Model parameters 181 182 % transform model to vector 183 if whatEst(1) == 'P' % 'P,CalcF' 184 P=init.P; 185 X0 = tr_p2x(P,1); 186 r = length(init.P); 187 else 188 X0 = []; 189 r = length(known.P); 190 end 191 192 for i = 1:r 193 X = feval(known.trModel2X,init.M{i}); 194 nM(i) = length(X); 195 X0 = [X0; X]; 196 end 197 198 known.nM = nM; 199 200 try 201 % For Matlab 5.3 and higher ??? 202 [X,OPTIONS] = fminsearch('f_smctp',X0,OPTIONS,Fobs,whatEst,method,known,whatKnown,init); 203 catch 204 % For Matlab 5.2 and lower ??? 205 [X,OPTIONS] = fmins('f_smctp',X0,OPTIONS,[],Fobs,whatEst,method,known,whatKnown,init); 206 end 207 208 % transform vector to model 209 if whatEst(1) == 'P' % 'P,CalcF' 210 r = length(init.P); 211 X1 = X(1:r*(r-1)); 212 X2 = X(r*(r-1)+1:end); 213 P = tr_x2p(X1,1); 214 Est.P = P; 215 else 216 P = known.P; 217 X2 = X; 218 end 219 r = length(P); 220 221 % transform vector to model 222 k1=1; 223 for i = 1:r 224 k2 = k1+nM(i)-1; 225 M{i} = feval(known.trX2Model,X2(k1:k2),known); 226 k1=k2+1; 227 end 228 229 Est.M = M; 230 231 232 case {'SimF','P,SimF'} % Estimate P, Model parameters 233 234 r = length(known.P); 235 236 % 1. Initial estimate 237 238 M = init.M; 239 240 % 2. Simulate each subload 241 242 F = cell(2,1); 243 for i = 1:r 244 F{i} = feval(known.simFun,known.param,M{i},known.T,known.T0)/(known.T/2); 245 end 246 247 % while ~slut 248 249 % 3. Uppdatera skattning 250 251 % transform model to vector 252 X0 = []; 253 for i = 1:r 254 X = feval(known.trModel2X,init.M{i}); 255 nM(i) = length(X); 256 X0 = [X0; X]; 257 end 258 259 known.nM = nM; 260 known.F = F; 261 262 % Estimate 263 try 264 % For Matlab 5.3 and higher ??? 265 [X,OPTIONS] = fminsearch('f_smctp',X0,OPTIONS,Fobs,whatEst,method,known,whatKnown,init); 266 catch 267 % For Matlab 5.2 and lower ??? 268 [X,OPTIONS] = fmins('f_smctp',X0,OPTIONS,[],Fobs,whatEst,method,known,whatKnown,init); 269 end 270 271 % transform vector to model 272 k1=1; 273 for i = 1:r 274 k2 = k1+nM(i)-1; 275 M{i} = feval(known.trX2Model,X(k1:k2),known); 276 k1=k2+1; 277 end 278 279 % 4. Simulate each subload and update 280 281 for i = 1:r 282 F1 = feval(known.simFun,M{i},known.T,known.T0)/(known.T/2); 283 F{i} = known.theta*F{i} + (1-known.theta)*F1; 284 end 285 286 Est.M = M; 287 % end 288 289 otherwise 290 291 % This should not happen 292 error(['Unexpected whatEst: ' whatEst '.']) 293 294 end % switch 295 296 % Calculate Estimated SMCTP-model 297 298 299 [y,F,P,FF] = f_smctp(X,Fobs,whatEst,method,known,whatKnown,init); 300 301 %Fest = smctp(P,FF); 302 Fest.P = P; 303 Fest.F = FF; 304 305 fprintf(1,'\n');
Comments or corrections to the WAFO group