Home > wafo > trgauss > mvnortpcprb.m

# mvnortpcprb

## PURPOSE

Multivariate normal or student T probability with product correlation structure.

## SYNOPSIS

[value,bound,inform] = mvnortpcprb(RHO,A,B,D,NDF,abseps,IERC,HNC);

## DESCRIPTION

``` MVNORTPCPRB Multivariate normal or student T probability with product correlation structure.

CALL [value,bound,inform] = mvnortpcprb(RHO,A,B,D,NDF,abseps,IERC,HINC)

RHO    = REAL, array of coefficients defining the correlation
coefficient by:
correlation(I,J) =  RHO(I)*RHO(J) for J/=I
where
1 < RHO(I) < 1
A         = vector of lower integration limits.
B         = vector of upper integration limits.
NOTE: any values greater the 37 in magnitude,
are considered as infinite values.
D      = vector of means (default zeros(size(RHO)))
NDF    = Degrees of freedom, NDF<=0 gives normal probabilities (default)
ABSEPS = absolute error tolerance. (default 1e-4)
IERC   = 1 if strict error control based on fourth derivative
0 if intuitive error control based on halving the
intervals (default)
HINC   = start interval width of simpson rule (default 0.24)

OUTPUT:
VALUE  = estimated value for the integral
BOUND  = bound on the error of the approximation
INFORM = INTEGER, termination status parameter:
0, if normal completion with ERROR < EPS;
1, if N > 1000 or N < 1.
2, IF  any abs(rho)>=1
4, if  ANY(b(I)<=A(i))
5, if number of terms computed exceeds maximum number of
evaluation points
6, if fault accurs in normal subroutines
7, if subintervals are too narrow or too many
8, if bounds exceeds abseps

MVNORTPCPRB calculates multivariate normal or student T probability
with product correlation structure for rectangular regions.
The accuracy is as best around single precision, i.e., about 1e-7.

Example:
rho2 = rand(1,2);
a2   = zeros(1,2);
b2   = repmat(inf,1,2);
[val2,err2] = mvnortpcprb(rho2,a2,b2);
g2 = inline('0.25+asin(x(1)*x(2))/(2*pi)');
E2 = g2(rho2)  % exact value

rho3 = rand(1,3);
a3   = zeros(1,3);
b3   = repmat(inf,1,3);
[val3,err] = mvnortpcprb(rho3,a3,b3);
g3 = inline('0.5-sum(sort(acos([x(1)*x(2),x(1)*x(3),x(2)*x(3)])))/(4*pi)');
E3 = g3(rho3)   %  Exact value

## CROSS-REFERENCE INFORMATION

This function calls:
 error Display message and abort function. mvnprdmex mvnprdmex53 str2num Convert string matrix to numeric array. version MATLAB version number.
This function is called by:

## SOURCE CODE

```001 function [value,bound,inform] = mvnortpcprb(RHO,A,B,D,NDF,abseps,IERC,HNC);
002 %MVNORTPCPRB Multivariate normal or student T probability with product correlation structure.
003 %
004 %  CALL [value,bound,inform] = mvnortpcprb(RHO,A,B,D,NDF,abseps,IERC,HINC)
005 %
006 %     RHO    = REAL, array of coefficients defining the correlation
007 %              coefficient by:
008 %                correlation(I,J) =  RHO(I)*RHO(J) for J/=I
009 %              where
010 %                1 < RHO(I) < 1
011 %     A         = vector of lower integration limits.
012 %     B         = vector of upper integration limits.
013 %           NOTE: any values greater the 37 in magnitude,
014 %                    are considered as infinite values.
015 %     D      = vector of means (default zeros(size(RHO)))
016 %     NDF    = Degrees of freedom, NDF<=0 gives normal probabilities (default)
017 %     ABSEPS = absolute error tolerance. (default 1e-4)
018 %     IERC   = 1 if strict error control based on fourth derivative
019 %              0 if intuitive error control based on halving the
020 %                      intervals (default)
021 %     HINC   = start interval width of simpson rule (default 0.24)
022 %
023 % OUTPUT:
024 %     VALUE  = estimated value for the integral
025 %     BOUND  = bound on the error of the approximation
026 %     INFORM = INTEGER, termination status parameter:
027 %            0, if normal completion with ERROR < EPS;
028 %            1, if N > 1000 or N < 1.
029 %            2, IF  any abs(rho)>=1
030 %            4, if  ANY(b(I)<=A(i))
031 %            5, if number of terms computed exceeds maximum number of
032 %                  evaluation points
033 %            6, if fault accurs in normal subroutines
034 %            7, if subintervals are too narrow or too many
035 %            8, if bounds exceeds abseps
036 %
037 % MVNORTPCPRB calculates multivariate normal or student T probability
038 % with product correlation structure for rectangular regions.
039 % The accuracy is as best around single precision, i.e., about 1e-7.
040 %
041 % Example:
042 %  rho2 = rand(1,2);
043 %  a2   = zeros(1,2);
044 %  b2   = repmat(inf,1,2);
045 %  [val2,err2] = mvnortpcprb(rho2,a2,b2);
046 %  g2 = inline('0.25+asin(x(1)*x(2))/(2*pi)');
047 %  E2 = g2(rho2)  % exact value
048 %
049 %  rho3 = rand(1,3);
050 %  a3   = zeros(1,3);
051 %  b3   = repmat(inf,1,3);
052 %  [val3,err] = mvnortpcprb(rho3,a3,b3);
053 %  g3 = inline('0.5-sum(sort(acos([x(1)*x(2),x(1)*x(3),x(2)*x(3)])))/(4*pi)');
054 %  E3 = g3(rho3)   %  Exact value
055 %
057
058 % Reference:
059 %   Charles Dunnett (1989)
060 %  "Multivariate normal probability integrals with product correlation
061 %  structure", Applied statistics, Vol 38,No3, (Algorithm AS 251)
062
063 % The mex-interface and m-file was written by
064 %     Per Andreas Brodtkorb
065 %     Norwegian Defence Research Establishment
066 %     P.O. Box 115
067 %     N-3191 Horten
068 %     Norway
069 %     Email: Per.Brodtkorb@ffi.no
070 %
071   error(nargchk(3,8,nargin))
072   if nargin<4|isempty(D),
073     D = zeros(size(RHO));
074   end
075   if nargin<5|isempty(NDF),
076      NDF = 0;
077   end
078   if nargin<6|isempty(abseps),
079      abseps = 1.0e-4;
080   end
081   if nargin<7|isempty(IERC),
082      IERC = 0;
083   end
084   if nargin<8|isempty(HNC),
085      HNC = 0.24;
086   end
087   % Make sure integration limits are finite
088   A = min(max(A,-100),100);
089   B = max(min(B,100),-100);
090   verstr = version;
091   ind = find(verstr=='.');
092   versionNumber = str2num(verstr(1:ind(2)-1));
093   if versionNumber<=5.3
094     [value,bound,inform] = mvnprdmex53(RHO,A,B,D,NDF,abseps,IERC,HNC);
095   else
096     [value,bound,inform] = mvnprdmex(RHO,A,B,D,NDF,abseps,IERC,HNC);
097   end```

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