Home > wafo > markov > estsmctp.m

estsmctp

PURPOSE ^

Estimate SMCTP model from an observed rainflow matrix.

SYNOPSIS ^

[Fest,Est,OPTIONS] = estsmctp(Fobs,whatEst,method,known,whatKnown,init,OPTIONS)

DESCRIPTION ^

 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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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');

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