CHITWO2LC_SORM SORM-approximation of the crossing intensity for the noncentral Chi^2 process The noncentral Chi^2 process is defined in IR-(27). X(t)=sum beta_j Z_j(t) + gamma_j Z_j(t)^2 CALL: mu = chitwo2lc_sorm(x0,gamma,beta,S12,S22,ass); mu = two column vector containing: levels in the first column, the SORM-approximation of crossing intensitity for quadratic sea in the second column. u = column with levels, note that u can not take value zero. g,b = the vectors containing constants defining the process approximation of the cdf of the quadratic sea. S12 = covariance between vector Z(t) and Z'(t), where Z(t)=(Z_1(t),...,Z_n(t)). S22 = covariance between vector Z'(t) and Z'(t)., where Z(t)=(Z_1(t),...,Z_n(t)). ass = 0, SORM, 0< FORM, >0 higher order very slow, optional ass=0; n0 = parameter used to find a starting point in the optimization, default value n0=10 epsi = tolerance level in the optimization, default value epsi=1e-9 Example: S=jonswap; [gamma beta S12 S22]=dirsp2chitwo(S.S,S.w); L0=sqrt(sum(beta.^2)); u=(0.01:0.1:4*L0)'; u=[(-4*L0:0.1:-0.1)';u]; mu = chitwo2lc_sorm(u,gamma,beta,S22,S12); semilogy(mu(:,1),mu(:,2))
Finds minimal distance to the origin on the surface b'*x+x'*diag(g)*x=u | |
Gives first two terms in an asymptotic expansion of the | |
Display message and abort function. | |
Kronecker tensor product. |
001 function mu=chitwo2lc_sorm(u,g,b,S22,S12,ass,n0,epsi) 002 % CHITWO2LC_SORM SORM-approximation of the crossing intensity for the noncentral Chi^2 process 003 % 004 % The noncentral Chi^2 process is defined in IR-(27). 005 % X(t)=sum beta_j Z_j(t) + gamma_j Z_j(t)^2 006 % 007 % CALL: mu = chitwo2lc_sorm(x0,gamma,beta,S12,S22,ass); 008 % 009 % mu = two column vector containing: levels in the first column, 010 % the SORM-approximation of crossing intensitity for quadratic sea 011 % in the second column. 012 % u = column with levels, note that u can not take value zero. 013 % g,b = the vectors containing constants defining the process 014 % approximation of the cdf of the quadratic sea. 015 % S12 = covariance between vector Z(t) and Z'(t), where Z(t)=(Z_1(t),...,Z_n(t)). 016 % S22 = covariance between vector Z'(t) and Z'(t)., where Z(t)=(Z_1(t),...,Z_n(t)). 017 % ass = 0, SORM, 0< FORM, >0 higher order very slow, optional ass=0; 018 % n0 = parameter used to find a starting point in the optimization, default value n0=10 019 % epsi = tolerance level in the optimization, default value epsi=1e-9 020 % 021 % Example: S=jonswap; [gamma beta S12 S22]=dirsp2chitwo(S.S,S.w); 022 % L0=sqrt(sum(beta.^2)); u=(0.01:0.1:4*L0)'; u=[(-4*L0:0.1:-0.1)';u]; 023 % mu = chitwo2lc_sorm(u,gamma,beta,S22,S12); 024 % semilogy(mu(:,1),mu(:,2)) 025 % 026 027 028 % References: 029 % K. Breitung, (1988), "Asymptotic crossing rates for stationary {G}aussian vector 030 % processes.", Stochastic Processes and 031 % their Applications, 29,pp. 195-207. 032 % 033 % U. Machado, I. Rychlik (2002) "Wave statistics in nonlinear sea", Extremes, 6, pp. 125--146. 034 % Hagberg, O. (2005) PhD - thesis, Dept. of Math. Statistics, Univ. of Lund. 035 % Baxevani, A., Hagberg, O. and Rychlik, I. Note on the distribution of extreme waves crests (OMAE 2005). 036 037 038 % Calls: mindist 039 % Revised by I.R 14.04.05 040 % By OH. 24.10.04 041 % 042 %------------------------------------------------------------------------------------ 043 044 if nargin<7 045 n0=10; 046 end 047 if nargin<8 048 epsi=1e-9; 049 end 050 % 051 % if ass>0 one is using second order assymptotics, see PhD thesis of Hagberg 052 % 053 if nargin<6 054 ass=0; 055 end 056 057 n=length(u); 058 x0=[]; 059 mu=[]; 060 for i=1:n; 061 xx1=mindist(g,b,u(i),n0,epsi); 062 x0=[x0 xx1]; 063 if ass>0 064 [Cb0 Cb1]=rqlf_asympt(u(i),g,b,S22,S12); 065 %mu=[mu; u(i) muu]; 066 mu=[mu; u(i) 0.5*(Cb0+Cb1)]; 067 end 068 end 069 if ass<=0 070 071 g=g(:); 072 b=b(:); 073 d=length(g); 074 [d1,n]=size(x0); 075 if d1~=d 076 error('the observations of x0 should be stored as columns') 077 end 078 b0=sqrt(sum(x0.^2,1)); 079 080 G=g*(-2*b0./sqrt(sum((kron(b,ones(1,n))+2*kron(g,ones(1,n)).* ... 081 x0).^2))); 082 % cofactor matrix: 083 D=zeros(size(x0)); 084 for k=1:d 085 D(k,:)=prod(1+G([1:k-1 k+1:d],:),1); 086 end 087 088 cc=[sqrt((sum(x0.*(S22*x0),1)+sum(G.*(S12'*x0).^2,1))) sqrt(sum(D.*x0.^2,1))]; 089 090 mu=exp(-b0.^2/2).*sqrt((sum(x0.*(S22*x0),1)+sum(G.*(S12'*x0).^2,1))./(sum(D.*x0.^2,1)))/pi/2; 091 092 mu1=exp(-b0.^2/2);%.*sqrt((sum(x0.*((S22)*x0),1))./(sum(x0.^2,1)))/2/pi; 093 mu2=exp(-b0.^2/2).*sqrt((sum(x0.*((S22)*x0),1))./(sum(x0.^2,1)))/2/pi; 094 095 cc=sqrt((sum(x0.*(S22*x0),1)+sum(G.*(S12'*x0).^2,1))./(sum(D.*x0.^2,1)))/pi/2; 096 mu=[u mu']; 097 if ass<0 098 mu=[u mu2']; 099 end 100 end
Comments or corrections to the WAFO group