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
Display message and abort function. | |
Convert string matrix to numeric array. | |
MATLAB version number. |
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
Comments or corrections to the WAFO group