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,
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
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:
 createpdf PDF class constructor kchitwo Computes the cumulant generating function formula for noncentral Chi2-process. wnormcdf Normal cumulative distribution function wnormpdf Normal probability density function error Display message and abort function. trapz Trapezoidal numerical integration.
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,
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