Home > wafo > trgauss > rind.m

# rind

## PURPOSE

Computes multivariate normal expectations

## SYNOPSIS

[fxind, err,exTime,options] = rind(BIG,Ex,Blo,Bup,indI,xc,Nt,varargin)

## DESCRIPTION

``` 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)

## CROSS-REFERENCE INFORMATION

This function calls:
 initoptions Initializes RIND options according to speed. rindoptset Create or alter RIND OPTIONS structure. wafoexepath Returns the path to executables for the WAFO Toolbox wnorminv Inverse of the Normal distribution function cell2struct Convert cell array to structure array. cellfun Functions on cell array contents. class Create object or return object class. clock Current date and time as date vector. delete Delete file or graphics object. dos Execute DOS command and return result. error Display message and abort function. etime Elapsed time. exist Check if variables or functions are defined. fclose Close file. fieldnames Get structure field names. fopen Open file. isequal True if arrays are numerically equal. lower Convert string to lowercase. mexrind71 mexrindalan22 mexrindalan24 nan Not-a-Number. struct Create or convert to structure array. struct2cell Convert structure array to cell array.
This function is called by:
 spec2mmtpdf Calculates joint density of Maximum, minimum and period. spec2thpdf Joint density of amplitude and period/wave-length characteristics spec2tpdf2 Evaluates densities for various wave periods or wave lengths

## SOURCE CODE

```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 %
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
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,...
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,...
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;
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