RIND Computes multivariate normal expectations CALL [E, err,exTime,options] = rind(S,m,Blo,Bup,indI,xc,Nt,options); E = expectation/density as explained below size 1 x Nx err = estimated absolute error, with 99% confidence level. size 1 x Nx exTime = execution time S = Covariance matrix of X=[Xt;Xd;Xc] size Ntdc x Ntdc (Ntdc=Nt+Nd+Nc) m = the expectation of X=[Xt;Xd;Xc] size N x 1 Blo,Bup = Lower and upper barriers used to compute the integration limits, Hlo and Hup, respectively. size Mb x Nb indI = vector of indices to the different barriers in the indicator function, length NI, where NI = Nb+1 (NB! restriction indI(1)=0, indI(NI)=Nt+Nd ) (default indI = 0:Nt+Nd) xc = values to condition on (default xc = zeros(0,1)), size Nc x Nx Nt = size of Xt (default Nt = Ntdc - Nc) options = rindoptions structure or named parameters with corresponding values, see rindoptset for details RIND computes multivariate normal expectations, i.e., E[Jacobian*Indicator|Condition ]*f_{Xc}(xc(:,ix)) where "Indicator" = I{ Hlo(i) < X(i) < Hup(i), i = 1:N_t+N_d } "Jacobian" = J(X(Nt+1),...,X(Nt+Nd+Nc)), special case is "Jacobian" = |X(Nt+1)*...*X(Nt+Nd)|=|Xd(1)*Xd(2)..Xd(Nd)| "condition" = Xc=xc(:,ix), ix=1,...,Nx. X = [Xt; Xd; Xc], a stochastic vector of Multivariate Gaussian variables where Xt,Xd and Xc have the length Nt, Nd and Nc, respectively. (Recommended limitations Nx,Nt<=100, Nd<=6 and Nc<=10) Multivariate probability is computed if Nd = 0. If Mb<Nc+1 then Blo(Mb+1:Nc+1,:) is assumed to be zero. The relation to the integration limits Hlo and Hup are as follows Hlo(i) = Blo(1,j)+Blo(2:Mb,j).'*xc(1:Mb-1,ix), Hup(i) = Bup(1,j)+Bup(2:Mb,j).'*xc(1:Mb-1,ix), where i=indI(j-1)+1:indI(j), j=2:NI, ix=1:Nx NOTE : RIND is only using upper triangular part of covariance matrix, S (except for options.method=0). Examples:% A) Compute Prob{Xi<-1.2} for i=1:5 where % Xi are zero mean Gaussian variables with covariance % Cov(Xi,Xj) = 0.3 for i~=j and % Cov(Xi,Xi) = 1 otherwise % B) Compute expectation E( X1^{+}*X2^{+} ) with random % correlation coefficient,Cov(X1,X2) = rho2. N = 5; Blo =-inf; Bup=-1.2; indI=[0 N]; % Barriers Hlo = repmat(Blo,1,N); Hup = repmat(Bup,1,N); % Integration limits m = zeros(N,1); rho = 0.3; Sc =(ones(N)-eye(N))*rho+eye(N); E0 = rind(Sc,m,Blo,Bup,indI) % exact prob. 0.001946 A) E1 = rind(triu(Sc),m,Hlo,Hup) % same as E0 m2 = [0 0]; rho2 = rand(1); Sc2 = [1 rho2; rho2,1]; Blo2 = 0; Bup2 = inf; indI2 = [0 2]; Nt2 = 0; opt2 = rindoptset('method',1); g2 = inline('(x*(pi/2+asin(x))+sqrt(1-x^2))/(2*pi)'); E2 = g2(rho2) % exact value B) E3 = rind(Sc2,m2,Blo2,Bup2,indI2,[],Nt2) E4 = rind(Sc2,m2,Blo2,Bup2,indI2,[],Nt2,opt2) E5 = rind(Sc2,m2,Blo2,Bup2,indI2,[],Nt2,'abseps', 1e-4) See also mvnormprb, mvnormpcprb, rindoptset
Initializes RIND options according to speed. | |
Create or alter RIND OPTIONS structure. | |
Returns the path to executables for the WAFO Toolbox | |
Inverse of the Normal distribution function | |
Convert cell array to structure array. | |
Functions on cell array contents. | |
Create object or return object class. | |
Current date and time as date vector. | |
Delete file or graphics object. | |
Execute DOS command and return result. | |
Display message and abort function. | |
Elapsed time. | |
Check if variables or functions are defined. | |
Close file. | |
Get structure field names. | |
Open file. | |
True if arrays are numerically equal. | |
Convert string to lowercase. | |
Not-a-Number. | |
Create or convert to structure array. | |
Convert structure array to cell array. |
Calculates joint density of Maximum, minimum and period. | |
Joint density of amplitude and period/wave-length characteristics | |
Evaluates densities for various wave periods or wave lengths |
0001 function [fxind, err,exTime,options] = rind(BIG,Ex,Blo,Bup,indI,xc,Nt,varargin) 0002 %RIND Computes multivariate normal expectations 0003 % 0004 % CALL [E, err,exTime,options] = rind(S,m,Blo,Bup,indI,xc,Nt,options); 0005 % 0006 % E = expectation/density as explained below size 1 x Nx 0007 % err = estimated absolute error, with 99% confidence level. size 1 x Nx 0008 % exTime = execution time 0009 % S = Covariance matrix of X=[Xt;Xd;Xc] size Ntdc x Ntdc (Ntdc=Nt+Nd+Nc) 0010 % m = the expectation of X=[Xt;Xd;Xc] size N x 1 0011 % Blo,Bup = Lower and upper barriers used to compute the integration 0012 % limits, Hlo and Hup, respectively. size Mb x Nb 0013 % indI = vector of indices to the different barriers in the 0014 % indicator function, length NI, where NI = Nb+1 0015 % (NB! restriction indI(1)=0, indI(NI)=Nt+Nd ) 0016 % (default indI = 0:Nt+Nd) 0017 % xc = values to condition on (default xc = zeros(0,1)), size Nc x Nx 0018 % Nt = size of Xt (default Nt = Ntdc - Nc) 0019 % options = rindoptions structure or named parameters with corresponding 0020 % values, see rindoptset for details 0021 % 0022 % RIND computes multivariate normal expectations, i.e., 0023 % E[Jacobian*Indicator|Condition ]*f_{Xc}(xc(:,ix)) 0024 % where 0025 % "Indicator" = I{ Hlo(i) < X(i) < Hup(i), i = 1:N_t+N_d } 0026 % "Jacobian" = J(X(Nt+1),...,X(Nt+Nd+Nc)), special case is 0027 % "Jacobian" = |X(Nt+1)*...*X(Nt+Nd)|=|Xd(1)*Xd(2)..Xd(Nd)| 0028 % "condition" = Xc=xc(:,ix), ix=1,...,Nx. 0029 % X = [Xt; Xd; Xc], a stochastic vector of Multivariate Gaussian 0030 % variables where Xt,Xd and Xc have the length Nt, Nd and Nc, 0031 % respectively. 0032 % (Recommended limitations Nx,Nt<=100, Nd<=6 and Nc<=10) 0033 % 0034 % Multivariate probability is computed if Nd = 0. 0035 % 0036 % If Mb<Nc+1 then Blo(Mb+1:Nc+1,:) is assumed to be zero. 0037 % The relation to the integration limits Hlo and Hup are as follows 0038 % 0039 % Hlo(i) = Blo(1,j)+Blo(2:Mb,j).'*xc(1:Mb-1,ix), 0040 % Hup(i) = Bup(1,j)+Bup(2:Mb,j).'*xc(1:Mb-1,ix), 0041 % 0042 % where i=indI(j-1)+1:indI(j), j=2:NI, ix=1:Nx 0043 % 0044 % NOTE : RIND is only using upper triangular part of covariance matrix, S 0045 % (except for options.method=0). 0046 % 0047 % Examples:% A) Compute Prob{Xi<-1.2} for i=1:5 where 0048 % % Xi are zero mean Gaussian variables with covariance 0049 % % Cov(Xi,Xj) = 0.3 for i~=j and 0050 % % Cov(Xi,Xi) = 1 otherwise 0051 % % B) Compute expectation E( X1^{+}*X2^{+} ) with random 0052 % % correlation coefficient,Cov(X1,X2) = rho2. 0053 % 0054 % N = 5; 0055 % Blo =-inf; Bup=-1.2; indI=[0 N]; % Barriers 0056 % Hlo = repmat(Blo,1,N); Hup = repmat(Bup,1,N); % Integration limits 0057 % m = zeros(N,1); rho = 0.3; 0058 % Sc =(ones(N)-eye(N))*rho+eye(N); 0059 % 0060 % E0 = rind(Sc,m,Blo,Bup,indI) % exact prob. 0.001946 A) 0061 % E1 = rind(triu(Sc),m,Hlo,Hup) % same as E0 0062 % 0063 % m2 = [0 0]; rho2 = rand(1); 0064 % Sc2 = [1 rho2; rho2,1]; 0065 % Blo2 = 0; Bup2 = inf; indI2 = [0 2]; 0066 % Nt2 = 0; 0067 % opt2 = rindoptset('method',1); 0068 % g2 = inline('(x*(pi/2+asin(x))+sqrt(1-x^2))/(2*pi)'); 0069 % 0070 % E2 = g2(rho2) % exact value B) 0071 % E3 = rind(Sc2,m2,Blo2,Bup2,indI2,[],Nt2) 0072 % E4 = rind(Sc2,m2,Blo2,Bup2,indI2,[],Nt2,opt2) 0073 % E5 = rind(Sc2,m2,Blo2,Bup2,indI2,[],Nt2,'abseps', 1e-4) 0074 % 0075 % See also mvnormprb, mvnormpcprb, rindoptset 0076 0077 % References 0078 % Podgorski et al. (2000) 0079 % "Exact distributions for apparent waves in irregular seas" 0080 % Ocean Engineering, Vol 27, no 1, pp979-1016. 0081 % 0082 % P. A. Brodtkorb (2004), 0083 % Numerical evaluation of multinormal expectations 0084 % In Lund university report series 0085 % and in the Dr.Ing thesis: 0086 % The probability of Occurrence of dangerous Wave Situations at Sea. 0087 % Dr.Ing thesis, Norwegian University of Science and Technolgy, NTNU, 0088 % Trondheim, Norway. 0089 0090 0091 % 0092 % 0093 % Tested on: Matlab 5.3,DIGITAL UNIX Fortran90 compiler 0094 % History: 0095 % -revised pab Nov 2004 0096 % - removed unused code + Now Nc1c2 is passed to rindalan24 0097 % -revised pab May 2004 0098 % completely rewritten by pab 18.02.2003 0099 % - initdata is now obsolete and replaced with rindoptset and initoptions 0100 % Note: the order of input arguments changed once again. (hopefully for the 0101 % last time) 0102 % revised pab 18.02.2003 0103 % - replaced old call to rind.exe with a meximplementation of it. 0104 % revised ir 10.2 2001 File closed after call to Fortran (last fclose) 0105 % revised jr 16.2 2001 (i) speed added in the call of this routine 0106 % (ii) The first four lines after the declaration of 0107 % global variables added (call of initdata, speed) 0108 % revised ir 6.2.2001 Adapted to MATLAB 6 0109 % revised jr 1.2.2001 The definition of Nb, Mb changed again. 0110 % revised ir 15.11.2000 A bug fixed in the definition of Nb Mb. 0111 % revised ir 10.11.2000 Compilation with rind60.f 0112 % revised jr 09.11.2000 The call to init_data replaced with initdata 0113 %% revised pab 23.05.2000 replaced the matlab call with a call to the F90 0114 % program RINDD.exe 0115 % revised by Per A. Brodtkorb 14.05.1999 0116 % - No limitations on size of the inputs 0117 % - enabled recursive calls 0118 % - fixed some bugs 0119 % - added some additonal checks 0120 % - updated to fortran90 0121 % by Igor Rychlik 29.10.1998 (PROGRAM RIND11 --- Version 1.0) 0122 0123 %default options 0124 options = struct('method',3,... 0125 'xcscale',0,... 0126 'abseps',0,... 0127 'releps',1e-3,... 0128 'coveps',1e-10,... 0129 'maxpts',40000,... 0130 'minpts',0,... 0131 'seed',floor(rand*1e9),... 0132 'nit',1000,... 0133 'xcutoff',[],... 0134 'xsplit',1.5,... 0135 'quadno',[] ,... 0136 'speed',[],... 0137 'nc1c2',2); 0138 % If just 'defaults' passed in, return the default options in g 0139 if ((nargin==1) & (nargout <= 1) & isequal(BIG,'defaults')), 0140 fxind = options; 0141 return 0142 end 0143 error(nargchk(4,inf,nargin)); 0144 Ntdc = size(BIG,1); 0145 0146 if nargin<6 | isempty(xc), 0147 xc = zeros(0,1); 0148 end 0149 Nc = size(xc,1); 0150 if nargin<7 | isempty(Nt), 0151 Nt = Ntdc-Nc; 0152 end 0153 [Mb, Nb] = size(Blo); 0154 0155 Nd = Ntdc-Nt-Nc; 0156 Ntd = Nt+Nd; 0157 0158 if nargin<5 | isempty(indI), 0159 if Nb~=Ntd 0160 error('Inconsistent size of Blo and Bup') 0161 end 0162 indI = 0:Ntd; 0163 end 0164 NI = length(indI); 0165 %method,XcScale,ABSEPS,RELEPS,COVEPS,MAXPTS,MINPTS,seed,NIT,xCutOff 0166 Np = nargin-7; 0167 if (Np>0) % handle various formats for options input 0168 switch lower(class(varargin{1})) 0169 case {'char','struct'}, 0170 options = rindoptset(options,varargin{:}); 0171 case {'cell'} 0172 if Np==1 0173 options = rindoptset(options,varargin{1}{:}); 0174 else 0175 error('Invalid options') 0176 end 0177 case {'double'} 0178 % Make sure it is compatible with old calls 0179 %varargin is a cell vector with double values, i.e., 0180 %{method,XcScale,ABSEPS,RELEPS,COVEPS,MAXPTS,MINPTS,... 0181 % seed,NIT,xCutOff,XSPLIT,QUADNO,SPEED}; 0182 ind = findNonEmptyCells(varargin); 0183 if any(ind) 0184 opt0 = struct2cell(options); 0185 opt0(ind) = varargin(ind); 0186 options = cell2struct(opt0,fieldnames(options)); 0187 end 0188 otherwise 0189 error('Invalid options') 0190 end 0191 end 0192 0193 if isempty(options.xcutoff) 0194 truncError = 0.1* max(options.abseps); 0195 Nc1c2 = max(1,options.nc1c2); 0196 options.xcutoff = max(min(abs(wnorminv(truncError/(Nc1c2*2))),8.5),1.2); 0197 %options.abseps = max(options.abseps- truncError,0); 0198 %options.releps = max(options.releps- truncError,0); 0199 end 0200 0201 0202 0203 % Old call 0204 %if nargin<9 | isempty(SCIS), SCIS = 1; end 0205 %if nargin<10 | isempty(XcScale), XcScale = 0; end 0206 %if nargin<11 | isempty(ABSEPS), ABSEPS = 0; end 0207 %if nargin<12 | isempty(RELEPS), RELEPS = 1e-3; end 0208 %if nargin<13 | isempty(COVEPS), COVEPS = 1e-10; end 0209 %if nargin<14 | isempty(MAXPTS), MAXPTS = 40000; end 0210 %if nargin<15 | isempty(MINPTS), MINPTS = 0; end 0211 %if nargin<16 | isempty(seed), seed = floor(rand*1e10); end 0212 %if nargin<17 | isempty(NIT), NIT = 1000; end 0213 %if nargin<18 | isempty(xCutOff), 0214 % xCutOff = max(min(abs(wnorminv(max(RELEPS,ABSEPS)*10^(-1+0*(Nt>10)))),7),1.2); 0215 %end 0216 0217 0218 0219 0220 %funcmod1 0221 0222 if options.method>0 0223 %opt0 = {SCIS,XcScale,ABSEPS,RELEPS,COVEPS,MAXPTS,MINPTS,seed,NIT,xCutOff}; 0224 opt0 = struct2cell(options); 0225 opt0{1} = mod(opt0{1},10); 0226 0227 % INFIN = INTEGER, array of integration limits flags: size 1 x Nb 0228 % if INFIN(I) < 0, Ith limits are (-infinity, infinity); 0229 % if INFIN(I) = 0, Ith limits are (-infinity, Hup(I)]; 0230 % if INFIN(I) = 1, Ith limits are [Hlo(I), infinity); 0231 % if INFIN(I) = 2, Ith limits are [Hlo(I), Hup(I)]. 0232 infinity = 37; 0233 dev = sqrt(diag(BIG).'); % std 0234 ind = find(indI(2:end)); 0235 INFIN = repmat(2,1,length(indI)-1); 0236 INFIN(ind) = 2 - (Bup(1,ind) > infinity*dev(indI(ind+1)))-... 0237 2*(Blo(1,ind) <-infinity*dev(indI(ind+1))); 0238 0239 t0 = clock; 0240 0241 try 0242 if options.method>10 0243 [fxind,err] = mexrindalan22(BIG,Ex,indI,Blo,Bup,INFIN,xc,Nt,... 0244 opt0{1:10}); 0245 else 0246 % rindalan24 with removal of infis and stricter c1c2 0247 [fxind,err] = mexrindalan24(BIG,Ex,indI,Blo,Bup,INFIN,xc,Nt,opt0{[1:10,14]}); 0248 end 0249 catch 0250 error('Compile mexrindalan24 again') 0251 end 0252 exTime = etime(clock,t0); 0253 else 0254 0255 if isempty(options.speed) 0256 options.speed = 4; 0257 options = initoptions(options.speed,options); 0258 end 0259 0260 0261 0262 t0 = clock; 0263 try % mexcompiled function 0264 if options.method==0; 0265 NIT1 = options.nit; 0266 else 0267 NIT1 = options.method; 0268 end 0269 speed = options.speed; 0270 seed = options.seed; 0271 xcscale = options.xcscale; 0272 0273 fxind = mexrind71(BIG,Ex,xc,Nt,NIT1,speed,indI,Blo,Bup,seed,xcscale); 0274 err = repmat(NaN,size(fxind)); 0275 catch 0276 error('Compile mexrind71 again') 0277 fxind = callRindExe(BIG,Ex,indI,Blo,Bup,xc,Nt,options); 0278 if (Nc>0 & options.xcscale~=0), % scale the result 0279 CC = exp(options.xcscale); 0280 fxind = fxind*CC; 0281 end 0282 end 0283 exTime = etime(clock,t0); 0284 0285 end 0286 0287 return % main 0288 0289 function ind = findNonEmptyCells(cellArray) 0290 %FINDNONEMPTYCELLS Return index to non-empty cells 0291 try, % matlab 5.3 or higher 0292 ind = find(~cellfun('isempty',cellArray)).'; 0293 catch 0294 % Slow 0295 n = prod(size(cellArray)); 0296 ind = zeros(1,n); 0297 for ix = 1:n 0298 ind(ix) = isempty(cellArray{ix}); 0299 end 0300 ind = find(~ind); 0301 end 0302 return % findNonEmptyCells 0303 0304 function [fxind,tid] = callRindExe(BIG,Ex,indI,B_lo,B_up,xc,Nt,options) 0305 %CALLRINDEXE Call rind.exe from wafoexepath (slow) 0306 % 0307 % This is kept just in case 0308 if length(Nt)>=2 0309 Nj = Nt(2); 0310 Nt = Nt(1); 0311 else 0312 Nj=0; 0313 end 0314 Nj = min(Nj,max(Nt,0)); % make sure Nj<Nt 0315 0316 Ntdc = size(BIG,1); 0317 [Nc, Nx]=size(xc); 0318 [Mb, Nb]=size(B_lo); 0319 NI = length(indI); 0320 Nd = Ntdc-Nt-Nc; 0321 Ntd = Nt+Nd; 0322 0323 disp(' Writing data.') 0324 filename=['BIG.in']; 0325 if exist(filename) 0326 delete(filename) 0327 end 0328 fid=fopen(filename,'wt'); 0329 for ix=1:Ntdc 0330 fprintf(fid,'%12.10f \n',BIG(ix,:)); 0331 end 0332 fclose(fid); 0333 0334 filename=['Ex.in']; 0335 if exist(filename) 0336 delete(filename) 0337 end 0338 fid=fopen(filename,'wt'); 0339 fprintf(fid,'%12.10f \n',Ex); 0340 fclose(fid); 0341 0342 filename=['xc.in']; 0343 if exist(filename) 0344 delete(filename) 0345 end 0346 fid=fopen(filename,'wt'); 0347 fprintf(fid,'%12.10f \n',xc); 0348 fclose(fid); 0349 0350 filename=['indI.in']; 0351 if exist(filename) 0352 delete(filename) 0353 end 0354 fid=fopen(filename,'wt'); 0355 fprintf(fid,'%2.0f \n',indI); 0356 fclose(fid); 0357 0358 filename=['B_lo.in']; 0359 if exist(filename) 0360 delete(filename) 0361 end 0362 fid=fopen(filename,'wt'); 0363 fprintf(fid,'%12.10f \n',B_lo'); 0364 fclose(fid); 0365 0366 filename=['B_up.in']; 0367 fid=fopen(filename,'wt'); 0368 if exist(filename) 0369 delete(filename) 0370 end 0371 fid=fopen(filename,'wt'); 0372 fprintf(fid,'%12.10f \n',B_up'); 0373 fclose(fid); 0374 0375 speed = options.speed; 0376 XSPLT = options.xsplit; 0377 RELEPS = options.releps; 0378 EPSS = options.abseps; 0379 EPS2 = options.coveps; 0380 xCutOff = options.xcutoff; 0381 NIT = options.nit; 0382 SCIS = abs(options.method); 0383 seed = options.seed; 0384 N_int = options.quadno; 0385 rateLHD = 20; 0386 0387 if exist('sizeinfo.in'), 0388 delete sizeinfo.in 0389 end 0390 fid=fopen('sizeinfo.in','wt'); 0391 fprintf(fid,'%2.0f \n', speed); 0392 fprintf(fid,'%2.0f \n', Nt); 0393 fprintf(fid,'%2.0f \n', Nd); 0394 fprintf(fid,'%2.0f \n', Nc); 0395 fprintf(fid,'%2.0f \n', NI); 0396 fprintf(fid,'%2.0f \n', Mb); 0397 fprintf(fid,'%2.0f \n', Nx); 0398 fprintf(fid,'%2.0f \n', NIT); 0399 fprintf(fid,'%2.0f \n', Nj); 0400 fprintf(fid,'%2.0f \n', seed); 0401 fprintf(fid,'%2.0f \n', SCIS); 0402 fprintf(fid,'%2.0f \n', rateLHD); 0403 fprintf(fid,'%12.10f \n', XSPLT); 0404 fprintf(fid,'%12.10f \n', EPSS); % absolute error 0405 fprintf(fid,'%12.10f \n', EPS2); % 0406 fprintf(fid,'%12.10f \n', xCutOff); % truncation value 0407 fprintf(fid,'%12.10f \n', RELEPS); % relativ error 0408 fprintf(fid,'%2.0f \n', N_int(1)); 0409 fprintf(fid,'%2.0f \n', 1); % minQnr 0410 fprintf(fid,'%2.0f \n', N_int(1)); % Le2Qnr 0411 fclose(fid); 0412 0413 disp(' Starting Fortran executable.') 0414 t0 = clock; 0415 0416 dos([wafoexepath, 'rindd70.exe']); 0417 if nargout>1 0418 tid=etime(clock,t0); 0419 end 0420 fxind=load('rind.out'); 0421 return
Comments or corrections to the WAFO group