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)
PDF class constructor | |
Computes the cumulant generating function formula for noncentral Chi2-process. | |
Normal cumulative distribution function | |
Normal probability density function | |
Display message and abort function. | |
Trapezoidal numerical integration. |
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
Comments or corrections to the WAFO group