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   
    
  See also  mvnormpcprb, mvnormprb, rind

CROSS-REFERENCE INFORMATION ^

This function calls: 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 %   
056 % See also  mvnormpcprb, mvnormprb, rind 
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