Home > wafo > markov > f_smctp.m

# f_smctp

## PURPOSE

Auxiliary function used by ESTSMCTP

## SYNOPSIS

[y,F,P,FF] = f_smctp(X,Fobs,whatEst,method,known,whatKnown,init)

## DESCRIPTION

F_SMCTP  Auxiliary function used by ESTSMCTP

CALL:  [y,F,P,FF] = Fsmtp(X,Fobs,whatEst,method,known,whatKnown,init)

## CROSS-REFERENCE INFORMATION

This function calls:
 chi2cmat Chi-square distance of cycle matrix. hdcmat Hellinger distance of cycle matrix. klcmat Kullback-Leibler distance of cycle matrix. loglcmat log-Likelihood of cycle matrix. scalemat Scale and translate a cycle matrix. smctp2arfm Calculates the asymmetric rainflow matrix for a SMCTP. smctp2rfm Calculates the rainflow matrix for a SMCTP. tr_x2p Transforms a vector X to a transition matrix P. cell Create cell array. error Display message and abort function. num2str Convert number to string. (Fast version) strcmp Compare strings. triu Extract upper triangular part.
This function is called by:
 estsmctp Estimate SMCTP model from an observed rainflow matrix.

## SOURCE CODE

001 function [y,F,P,FF] = f_smctp(X,Fobs,whatEst,method,known,whatKnown,init)
002 %F_SMCTP  Auxiliary function used by ESTSMCTP
003 %
004 % CALL:  [y,F,P,FF] = Fsmtp(X,Fobs,whatEst,method,known,whatKnown,init)
005
006 %fprintf(1,'X = [%f',X(1));
007 %fprintf(1,', %f',X(2:length(X)));
008 %fprintf(1,']\n');
009
010
011 switch whatEst
012
013 case 'P'      % Estimate P
014
015   P = tr_x2p(X,1);
016   r = length(P);
017
018   FF = known.F;
019
020 case 'MeanStd' % Estimate Mean and Std
021
022   P = known.P;
023
024   r = length(X)/2;
025   MeanStd = reshape(X,r,2);
026   MeanStd(:,2) = exp(MeanStd(:,2));
027
028 %  subplot(2,2,1), plotcmat(F1,3)
029 %  subplot(2,2,2), plotcmat(F2,3)
030
031   n = length(known.F{1,1});
032   param = [-1 1 n];
033   FF = cell(r,2);
034   for i = 1:r
035     FF{i,1} = scalemat(param,known.F{i,1},MeanStd(i,1),MeanStd(i,2),param);
036   end
037
038 %  subplot(2,2,3), plotcmat(F1,3)
039 %  subplot(2,2,4), plotcmat(F2,3)
040 %  drawnow
041
042 case 'P,MeanStd' % Estimate P, Mean and Std
043
044   r=(-1+sqrt(1+4*length(X)))/2;
045   X1 = X(1:r*(r-1));
046   X2 = X(r*(r-1)+1:end);
047
048   P = tr_x2p(X1,1);
049
050   MeanStd = reshape(X2,r,2);
051   MeanStd(:,2) = exp(MeanStd(:,2));
052
053   n= length(known.F{1,1});
054   param = [-1 1 n];
055   FF = cell(r,2);
056   for i = 1:r
057     FF{i,1} = scalemat(param,known.F{i,1},MeanStd(i,1),MeanStd(i,2),param);
058   end
059
060
061 case {'P,CalcF','CalcF'} % Estimate Model parameters
062
063   if whatEst(1) == 'P' % 'P,CalcF'
064     r = length(init.P);
065     X1 = X(1:r*(r-1));
066     X = X(r*(r-1)+1:end);
067     P = tr_x2p(X1,1);
068   else
069     P = known.P;
070   end
071   r = length(P);
072
073   % transform vector to model
074   FF = cell(r,2);
075   k1=1;
076   for i = 1:r
077     k2 = k1+known.nM(i)-1;
078     M = feval(known.trX2Model,X(k1:k2),known);
079     FF{i,1} = feval(known.Ffun,known.param,M);
080     k1=k2+1;
081   end
082
083 case {'SimF'} % Estimate Model parameters
084
085   P = known.P;
086   r = length(P);
087
088   % transform vector to model
089   FF = cell(r,2);
090   k1=1;
091   for i = 1:r
092     k2 = k1+known.nM(i)-1;
093     M = feval(known.trX2Model,X(k1:k2),known);
094     FF{i,1} = feval(known.simFun,known.param,M,known.T,known.T0);
095     k1=k2+1;
096   end
097
098
099 otherwise
100
101   % This should not happen
102   error(['Unexpected whatEst: ' whatEst '.'])
103
104 end
105
106 if known.NOsubzero ~= 0
107   for i = 1:r % Set to zero below and on the sub-diagonals
108     FF{i,1} = flipud(triu(flipud(FF{i,1})',1+known.NOsubzero)');
109   end
110 end
111
112 if known.SideInfo == 0
113   F = smctp2rfm(P,FF);
114 elseif known.SideInfo == 11
115   [F,Fsid] = smctp2arfm(P,FF,1,1);
116   n = length(Fsid{1,1});
117   F = zeros(n*r,n*r);
118   FobsSid = Fobs;
119   Fobs = zeros(n*r,n*r);
120   for z = 1:r
121     for w = 1:r
122       Fsid{z,w}(1,1) = Fsid{z,w}(1,1) + Fsid{z,w}(n,n);
123       Fsid{z,w}(n,n) = 0;
124       FobsSid{z,w}(1,1) = FobsSid{z,w}(1,1) + FobsSid{z,w}(n,n);
125       FobsSid{z,w}(n,n) = 0;
126       F(1+n*(z-1):n*z,1+n*(w-1):n*w) = Fsid{z,w};
127       Fobs(1+n*(z-1):n*z,1+n*(w-1):n*w) = FobsSid{z,w};
128     end
129   end
130 else
131   error(['Method for SideInfo = ' num2str(known.SideInfo) ' not implemented']);
132 end
133
134 %F = F * sum(sum(Fobs));
135
136 if strcmp(method,'ML') == 1
137   y = -loglcmat(Fobs,F);
138 elseif strcmp(method,'chi2') == 1
139   y = chi2cmat(Fobs,F);
140 elseif strcmp(method,'HD') == 1
141   y = hdcmat(Fobs,F);
142 elseif strcmp(method,'KL') == 1
143   y = klcmat(Fobs,F);
144 else
145   fprintf(1,['Method ' method ' not implemented']);
146 end
147
148 %fprintf(1,'y=%f\n',y);
149 fprintf(1,'*');
150

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