Home > wafo > trgauss > chitwo2lc_sp.m

chitwo2lc_sp

PURPOSE ^

Saddlepoint approximation of the crossing intensity for the noncentral Chi^2 process

SYNOPSIS ^

[mu,mu2,mu1,muG,beta,cc,fu,Fu]= chitwo2lc_sp(gamma,beta,S12,S22,ds0)

DESCRIPTION ^

  CHITWO2LC_SP  Saddlepoint 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,fu,Fu]= chitwo2lc_sp(gamma,beta,S12,S22,ds);
   
        mu   = three column vector containing: levels in the first column, 
               the saddlepoint approximation of crossing intensitity for quadratic sea  
               in the second column and for Gaussian sea in the third column.
               Note that the levels $u$ are not equally spaced. If the spacing 
               in the grid is not fine enough, choose smaller value of parameter ds0. 
        fu   = the pdf structure: f.x contains levels u, f.f is a two column matrix containing
               the Danniel's-saddlepoint approximation of the density (pdf) of the quadratic 
               sea in the first column and pdf of the Gaussian sea in the second.  
  gamma,beta = the vector 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)).
      ds0    = parameter defining the levels: default value is ds0=0.1.
 
    Example: S=jonswap; [gamma beta S12 S22]=dirsp2chitwo(S.S,S.w);
             [mu, fu, Fu]=chitwo2lc_sp(gamma,beta,S12,S22); 
             semilogy(mu(:,1),mu(:,2)); figure; pdfplot(fu)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [mu,mu2,mu1,muG,beta,cc,fu,Fu]= chitwo2lc_sp(gamma,beta,S12,S22,ds0)
002 % CHITWO2LC_SP  Saddlepoint 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,fu,Fu]= chitwo2lc_sp(gamma,beta,S12,S22,ds);
008 %  
009 %       mu   = three column vector containing: levels in the first column, 
010 %              the saddlepoint approximation of crossing intensitity for quadratic sea  
011 %              in the second column and for Gaussian sea in the third column.
012 %              Note that the levels $u$ are not equally spaced. If the spacing 
013 %              in the grid is not fine enough, choose smaller value of parameter ds0. 
014 %       fu   = the pdf structure: f.x contains levels u, f.f is a two column matrix containing
015 %              the Danniel's-saddlepoint approximation of the density (pdf) of the quadratic 
016 %              sea in the first column and pdf of the Gaussian sea in the second.  
017 % gamma,beta = the vector containing constants defining the process
018 %              approximation of the cdf of the quadratic sea.  
019 %     S12    = covariance between vector Z(t) and Z'(t), where Z(t)=(Z_1(t),...,Z_n(t)).
020 %     S22    = covariance between vector Z'(t) and Z'(t)., where Z(t)=(Z_1(t),...,Z_n(t)).
021 %     ds0    = parameter defining the levels: default value is ds0=0.1.
022 %
023 %   Example: S=jonswap; [gamma beta S12 S22]=dirsp2chitwo(S.S,S.w);
024 %            [mu, fu, Fu]=chitwo2lc_sp(gamma,beta,S12,S22); 
025 %            semilogy(mu(:,1),mu(:,2)); figure; pdfplot(fu)
026 %
027 
028 
029 % References: U. Machado, I. Rychlik (2003) "Wave statistics in nonlinear
030 %                                                 sea", Extremes, 6, pp. 125--146.           
031 %             Butler, R., Machado, U. Rychlik, I. (2003) "Distribution of wave crests in non-
032 %             linear random sea - application of saddlepoint methods", presented at ISOPE 2003.
033 %             Hagberg, O. (2005) PhD - thesis, Dept. of Math. Statistics, Univ. of Lund. 
034 %   Calls: kchitwo
035 % Revised by I.R 14.04.05
036 % Revised by I.R 24.11.03
037 %
038 %------------------------------------------------------------------------------------
039 
040 
041 
042 if nargin<1 
043   error('problem is not defined');
044 end
045 
046 n=length(gamma);
047 if nargin<2 
048   beta=zeros(size(gamma));
049 end
050 variance=2*(sum(gamma.^2));
051 if nargin<3 
052   S12=zeros(n);
053 end
054 if nargin<4 
055   S22=eye(n);
056 end
057 if nargin<5 | ds0<=0
058   ds0=0.1;
059 end
060 
061   ucrt=5*sqrt(sum(beta.^2)+2*variance);
062  
063 
064 % ====================================================================================
065 
066 
067 K=[];
068 DK=[];
069 DDK=[];
070 DK111=[];
071 
072 DK1122=[];
073 DK122=[];
074 
075 DDKt=[];
076 
077 DK2222=[];
078 R0=[];
079 S=[];
080 
081 s=0;
082 j=0;
083 flag=0;
084 
085 ds1=ds0;
086 
087 % In these formulas dt must be equal to ds
088 
089 %==>positive side
090 K10old=0;
091 while flag~=1
092  
093   j=j+1;
094   
095   if (j>400)
096       flag=1;
097       disp('step dt is too small ')
098 %      error('step dt is too small - stop')
099   end 
100 % ds=0.2*ds1;
101  ds=0.05*ds1;
102  dt=ds;
103      
104  [K0]=kchitwo(s,0,gamma,beta,S12,S22);
105  [K1]=kchitwo(s-ds,0,gamma,beta,S12,S22);
106  [K2]=kchitwo(s+ds,0,gamma,beta,S12,S22);   
107  [K3]=kchitwo(s-2*ds,0,gamma,beta,S12,S22);
108  [K4]=kchitwo(s+2*ds,0,gamma,beta,S12,S22);
109  
110  [K1t]=kchitwo(s,-dt,gamma,beta,S12,S22);
111  [K2t]=kchitwo(s,dt,gamma,beta,S12,S22);
112  
113  [K11]=kchitwo(s-ds,-dt,gamma,beta,S12,S22);
114  [K22]=kchitwo(s+ds,dt,gamma,beta,S12,S22);
115  
116  [K12t]=kchitwo(s-ds,dt,gamma,beta,S12,S22);
117  [K21t]=kchitwo(s+ds,-dt,gamma,beta,S12,S22);
118  
119  [K2tt]=kchitwo(s,2*dt,gamma,beta,S12,S22);
120  [K1tt]=kchitwo(s,-2*dt,gamma,beta,S12,S22);
121 
122  K10=(K2-K1)/(2*ds);
123  K20=(K1+K2-2*K0)/(ds*ds);
124  
125  K30=(K4-2*K2+2*K1-K3)/(2*ds*ds*ds);
126  K40=(K4-4*K2+6*K0-4*K1+K3)/(ds*ds*ds*ds);
127 
128  K12=(K22-2*K2+K21t-K12t+2*K1-K11)/(2*ds*ds*ds);
129  
130  K02=(K1t+K2t-2*K0)/(dt*dt);
131  K04=(K2tt-4*K2t+6*K0-4*K1t+K1tt)/(dt*dt*dt*dt);
132  K22=(K22-2*K2t+K12t-2*K2+4*K0-2*K1+K21t-2*K1t+K11)/(ds*ds*ds*ds);  
133 [j K10];
134 dds=abs(K10-K10old);
135 K10old=K10;
136 if abs(dds)>0.075
137     ds1=ds1/2;
138 end
139  if (j~=1) & ((s*K10-K0)>25 | 0>(s*K10-K0) | abs(K10)>ucrt | K20<0)
140 
141 %if (j~=1) & (abs(K0-s*K10)>20)
142 %    if (j~=1) & (abs(K0-s*K10)>10) 
143   flag=1;
144 else
145 
146  S=[S; s]; 
147  s=s+ds1/2; 
148  
149  K=[K; K0];  
150  R0=[R0; (K40/K20^2/8-5*K30*K30/K20^3/24)];
151  DK=[DK; K10];
152  DDK=[DDK; K20];
153  
154  DK111=[DK111; K30];
155  DDKt=[DDKt; K02];
156  DK1122=[DK1122; K22]; 
157  DK122=[DK122; K12]; 
158  DK2222=[DK2222; K04];
159  end    
160 end     
161 
162 %=======================================================
163 
164 s=-ds0/2;
165 j=0;
166 flag=0;
167 ds1=ds0;
168 K10old=DK(1);
169 %==>negative side
170 ds=0.25*ds1;
171 dt=ds;
172 
173 while flag~=1
174  j=j+1;
175 
176  
177   if (j>400)
178       flag=1;
179       disp('step dt is too small ')
180 %      error('step dt is too small - stop')
181   end 
182  [K0]=kchitwo(s,0,gamma,beta,S12,S22);
183  [K1]=kchitwo(s-ds,0,gamma,beta,S12,S22);
184  [K2]=kchitwo(s+ds,0,gamma,beta,S12,S22);
185  [K3]=kchitwo(s-2*ds,0,gamma,beta,S12,S22);
186  [K4]=kchitwo(s+2*ds,0,gamma,beta,S12,S22);
187   
188  [K1t]=kchitwo(s,-dt,gamma,beta,S12,S22);
189  [K2t]=kchitwo(s,dt,gamma,beta,S12,S22); 
190  
191  [K11]=kchitwo(s-ds,-dt,gamma,beta,S12,S22);
192  [K22]=kchitwo(s+ds,dt,gamma,beta,S12,S22);
193 
194  [K12t]=kchitwo(s-ds,dt,gamma,beta,S12,S22);
195  [K21t]=kchitwo(s+ds,-dt,gamma,beta,S12,S22);
196  
197  [K2tt]=kchitwo(s,2*dt,gamma,beta,S12,S22);
198  [K1tt]=kchitwo(s,-2*dt,gamma,beta,S12,S22);
199  
200  K10=(K2-K1)/(2*ds);
201  K20=(K1+K2-2*K0)/(ds*ds);
202  
203  K30=(K4-2*K2+2*K1-K3)/(2*ds*ds*ds); 
204  K40=(K4-4*K2+6*K0-4*K1+K3)/(ds*ds*ds*ds);
205 
206  K12=(K22-2*K2+K21t-K12t+2*K1-K11)/(2*ds*ds*ds);
207  
208  K02=(K1t+K2t-2*K0)/(dt*dt);
209  K04=(K2tt-4*K2t+6*K0-4*K1t+K1tt)/(dt*dt*dt*dt);
210  K22=(K22-2*K2t+K12t-2*K2+4*K0-2*K1+K21t-2*K1t+K11)/(ds*ds*ds*ds); 
211  [j K10];
212  if (j~=1) & ((s*K10-K0)>25 | 0>(s*K10-K0) | abs(K10)>ucrt | K20<0)
213 %     ((s*K10-K0)>25 | 0>(s*K10-K0) )  & abs(K10)>ucrt  & K20<0
214 %         if (j~=1) & (abs(K0-s*K10)>10) 
215     flag=1;
216 else
217 dds=abs(K10-K10old);
218 K10old=K10;
219 if dds>0.1
220     ds1=ds1/2;
221 end
222  
223  S=[s;S];
224  s=s-ds1/2;
225  ds=0.25*ds1;
226  dt=ds;
227  K=[K0;K];
228  R0=[(K40/K20^2/8-5*K30*K30/K20^3/24); R0];
229 
230  DK=[K10;DK];
231  DDK=[K20;DDK];
232  
233  DK111=[K30;DK111];
234  DDKt=[K02; DDKt];
235  DK1122=[K22; DK1122]; 
236  DK122=[K12; DK122]; 
237  DK2222=[K04; DK2222];
238  
239 end
240 end     
241 
242 % ====================================================================================
243 % ====================================================================================
244 % marginal density
245 
246 f=1/sqrt(2*pi)*exp(K-S.*DK)./sqrt(DDK);
247 wh=sign(S).*sqrt(2*(S.*DK-K));
248 f1=wnormpdf(wh)./sqrt(DDK);
249 
250 [area]=trapz(DK,f);
251 beta=sqrt(2*abs(S.*DK-K)-2*log(area));
252 %area=1;
253 fu=[DK f/area]; 
254 fu=createpdf;
255   Htxt = ['Density of X(0) non-central Chi^2 process'];
256   xtxt = ['X(0)'];
257 fu.title=Htxt;
258 fu.labx{1}=xtxt;
259 fu.x{1}=DK;
260 fu.f=[f/area];   
261 
262 %fu1=[DK f1/area];
263 F1=wnormcdf(wh)-wnormpdf(wh).*((1./S)./sqrt(DDK)-1./wh);
264 Fu=[DK F1];
265 
266 % ====================================================================================
267 % 1-TERM SADDLEPONT APRROXIMATION (GAUSSIAN APPR)
268 % ====================================================================================
269 
270 muG=(f/area).*sqrt(DDKt)/sqrt(2*pi);
271 cc=sqrt(DDKt)/sqrt(2*pi);
272 %figure(3)
273 %plot(DK,muG,'k:o')
274 
275 % ====================================================================================
276 % 2-TERMS SADDLEPONT APRROXIMATION (GAUSSIAN APPR)
277 % ====================================================================================
278 
279 % equation 22
280 eq22=-1/4*(DK1122.*DDK-DK111.*DK122)./((DDK.^2).*DDKt);
281 
282 %mu_2=muG.*(1+eq22);
283 
284 %figure(3)
285 %hold on
286 %plot(DK,mu_2,'k:*')
287 
288 % ====================================================================================
289 % 3-TERMS SADDLEPONT APRROXIMATION (GAUSSIAN APPR)
290 % ====================================================================================
291 
292 % equation 23
293 eq23=(DK2222.*DDK-3*DK122.^2)./(DDK.*(DDKt.^2));
294 
295 mu1=muG.*(1+eq22);
296 mu=muG.*(1+R0).*(1+eq22-(1/24)*eq23);
297 mu2=muG.*(1+eq22-(1/24)*eq23.*(1-eq22));
298 %figure(3)
299 %hold on
300 %plot(DK,mu,'k:d')
301 
302 %mu1=DDKt;
303 %mu2=(DK2222.*DDK-3*DK122.^2)./DDK;
304 
305 mu=[DK mu];
306 muG=[DK muG];
307 mu1=[DK mu1];
308 mu2=[DK mu2];
309 
310 
311 
312 
313 
314 % ====================================================================================
315

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