WNORMINV Inverse of the Normal distribution function CALL: x = wnorminv(F,m,v) x = inverse cdf for the Normal distribution evaluated at F m = mean (default 0) v = variance (default 1) Example: F = linspace(0,1,100); x = wnorminv(F,2.5,0.6); plot(F,x) See also wnorminv
Check if all input arguments are either scalar or of common size. | |
Normal cumulative distribution function | |
Display message and abort function. | |
Not-a-Number. |
Parameter estimates for Beta-Rayleigh data. | |
Estimate transformation, g, from observed CDF. | |
Parameter estimates for DIST2D data. | |
Initializes RIND options according to speed. | |
Estimate transformation, g, from observed crossing intensity, version2. | |
Inverse of the conditional cdf of X2 given X1. | |
Parameter estimates for MDIST2D data. | |
Multivariate Normal probabilities with product correlation | |
Normal probability plot of effects | |
Parameter estimates and confidence intervals for Ochi data. | |
Computes multivariate normal expectations | |
Parameter estimates for Beta data. | |
Parameter estimates for Chi squared data. | |
Inverse of the conditional 2D weibull cdf of X2 given X1. | |
Parameter estimates for 2D Weibull data. | |
Parameter estimates for Gamma data. | |
Inverse of the Gamma distribution function | |
Parameter estimates for GEV data | |
Parameter estimates for Generalized Gamma data. | |
Inverse of the Generalized Gamma distribution function | |
Parameter estimates for Gumbel data. | |
Parameter estimates for Inverse Gaussian data. | |
Inverse of the Lognormal distribution function | |
Parameter estimates for Rayleigh data. | |
Parameter estimates for Student's T data. | |
Inverse of the Student's T distribution function | |
Parameter estimates for Truncated Rayleigh data. | |
Parameter estimates for truncated Weibull data. | |
Parameter estimates for Weibull data. |
001 function x = wnorminv(F,m,v) 002 %WNORMINV Inverse of the Normal distribution function 003 % 004 % CALL: x = wnorminv(F,m,v) 005 % 006 % x = inverse cdf for the Normal distribution evaluated at F 007 % m = mean (default 0) 008 % v = variance (default 1) 009 % 010 % Example: 011 % F = linspace(0,1,100); 012 % x = wnorminv(F,2.5,0.6); 013 % plot(F,x) 014 % 015 % See also wnorminv 016 017 % Tested on: Matlab 6.0, 5.3 018 % History: 019 % revised pab 23.03.2003 020 % -replace erfinv with PHIINV which perform more accurate 021 % inversion in the lower tail. This also results in faster execution. 022 % revised jr 03.04.2001 023 % - fixed a bug in the last if statement 024 % - updated information, example 025 % revised pab 24.10.2000 026 % - added comnsize, nargchk 027 % - added default values 028 % added ms 17.06.2000 029 030 error(nargchk(1,3,nargin)) 031 if nargin<2|isempty(m), m=0; end 032 if nargin<3|isempty(v), v=1; end 033 034 [errorcode, F, m, v] = comnsize (F,m, v); 035 if (errorcode > 0) 036 error ('F, m and v must be of common size or scalar'); 037 end 038 039 x=zeros(size(F)); 040 041 ok = ((v>0)&(F>=0)&(F<=1)); 042 k = find (ok); 043 if any(k) 044 % old call 045 %x(k) = sqrt(2*v(k)).*erfinv(2 * F(k) - 1) + m(k); 046 % new call 047 x(k) = sqrt(v(k)).*PHIINV(F(k))+m(k); 048 end 049 050 k1 = find (~ok); 051 if any(k1) 052 tmp=NaN; 053 x(k1) = tmp(ones(size(k1))); 054 end 055 056 return 057 058 059 function PHINV = PHIINV( P ) 060 % 061 % ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3 062 % 063 % Produces the normal deviate Z corresponding to a given lower 064 % tail area of P. 065 % 066 % Absolute error less than 1e-13 067 % Relative error less than 1e-15 for abs(VAL)>0.1 068 % 069 % The hash sums below are the sums of the mantissas of the 070 % coefficients. They are included for use in checking 071 % transcription. 072 % 073 074 SPLIT1 = 0.425; 075 SPLIT2 = 5; 076 CONST1 = 0.180625; 077 CONST2 = 1.6; 078 079 080 % Coefficients in rational approximations. 081 %! 082 %! Coefficients for P close to 0.5 083 %! 084 A0 = 3.3871328727963666080E00; 085 A1 = 1.3314166789178437745E+2; 086 A2 = 1.9715909503065514427E+3; 087 A3 = 1.3731693765509461125E+4; 088 A4 = 4.5921953931549871457E+4; 089 A5 = 6.7265770927008700853E+4; 090 A6 = 3.3430575583588128105E+4; 091 A7 = 2.5090809287301226727E+3; 092 093 B1 = 4.2313330701600911252E+1; 094 B2 = 6.8718700749205790830E+2; 095 B3 = 5.3941960214247511077E+3; 096 B4 = 2.1213794301586595867E+4; 097 B5 = 3.9307895800092710610E+4; 098 B6 = 2.8729085735721942674E+4; 099 B7 = 5.2264952788528545610E+3; 100 %* HASH SUM AB 55.88319 28806 14901 4439 101 %! 102 %! Coefficients for P not close to 0, 0.5 or 1. 103 %! 104 C0 = 1.42343711074968357734E00; 105 C1 = 4.63033784615654529590E00; 106 C2 = 5.76949722146069140550E00; 107 C3 = 3.64784832476320460504E00; 108 C4 = 1.27045825245236838258E00; 109 C5 = 2.41780725177450611770E-1; 110 C6 = 2.27238449892691845833E-2; 111 C7 = 7.74545014278341407640E-4; 112 D1 = 2.05319162663775882187E00; 113 D2 = 1.67638483018380384940E00; 114 D3 = 6.89767334985100004550E-1; 115 D4 = 1.48103976427480074590E-1; 116 D5 = 1.51986665636164571966E-2; 117 D6 = 5.47593808499534494600E-4; 118 D7 = 1.05075007164441684324E-9 ; 119 %* HASH SUM CD 49.33206 50330 16102 89036 120 %! 121 %! Coefficients for P near 0 or 1. 122 %! 123 E0 = 6.65790464350110377720E00; 124 E1 = 5.46378491116411436990E00; 125 E2 = 1.78482653991729133580E00; 126 E3 = 2.96560571828504891230E-1; 127 E4 = 2.65321895265761230930E-2; 128 E5 = 1.24266094738807843860E-3; 129 E6 = 2.71155556874348757815E-5; 130 E7 = 2.01033439929228813265E-7; 131 F1 = 5.99832206555887937690E-1; 132 F2 = 1.36929880922735805310E-1; 133 F3 = 1.48753612908506148525E-2; 134 F4 = 7.86869131145613259100E-4; 135 F5 = 1.84631831751005468180E-5; 136 F6 = 1.42151175831644588870E-7; 137 F7 = 2.04426310338993978564E-15; 138 % HASH SUM EF 47.52583 31754 92896 71629 139 140 PHINV = zeros(size(P)); 141 Q = ( 2*P - 1 )/2; 142 region1 = abs(Q) <= SPLIT1; 143 k1 = find(region1); 144 if any(k1) 145 R1 = CONST1 - Q(k1).^2; 146 PHINV(k1) = Q(k1).*( ( ( ((((A7.*R1 + A6).*R1 + A5).*R1 + A4).*R1 + A3) ... 147 .*R1 + A2 ).*R1 + A1 ).*R1 + A0 ) ... 148 ./( ( ( ((((B7.*R1 + B6).*R1 + B5).*R1 + B4).*R1 + B3) ... 149 .*R1 + B2 ).*R1 + B1 ).*R1 + 1 ); 150 end 151 152 k2 = find(~region1); 153 if any(k2) 154 R2 = min( P(k2), 1 - P(k2) ); 155 k3 = find(R2>0); 156 if any(k3) % if ( R2 > 0 ) %THEN 157 R3 = sqrt( -log(R2(k3)) ); 158 k4 = find(R3 <= SPLIT2); 159 if any(k4) %if ( R2 <= SPLIT2 ) % THEN 160 R4 = R3(k4) - CONST2; 161 PHINV(k2(k3(k4))) = ( ( ( ((((C7.*R4 + C6).*R4 + C5).*R4 + ... 162 C4).*R4 + C3).*R4 + C2 ).*R4 + ... 163 C1 ).*R4 + C0 ) ... 164 ./( ( ( ((((D7.*R4 + D6).*R4 + D5).*R4 + D4).*R4 + D3) ... 165 .*R4 + D2 ).*R4 + D1 ).*R4 + 1 ); 166 end 167 k5 = find(R3 > SPLIT2); 168 if any(k5) 169 R5 = R3(k5) - SPLIT2; 170 PHINV(k2(k3(k5))) = ( ( ( ((((E7.*R5 + E6).*R5 + E5).*R5 + ... 171 E4).*R5 + E3).*R5 + E2 ).*R5 + ... 172 E1 ).*R5 + E0 ) ... 173 ./( ( ( ((((F7.*R5 + F6).*R5 + F5).*R5 + F4).*R5 + F3) ... 174 .*R5 + F2 ).*R5 + F1 ).*R5 + 1 ); 175 end % IF 176 end 177 k6 = find(R2<=0); 178 if any(k6) 179 PHINV(k2(k6)) = inf; 180 end %IF 181 k7 = find(Q(k2)<0); 182 if any(k7) %( Q < 0 ) %THEN 183 PHINV(k2(k7)) = - PHINV(k2(k7)); 184 end %IF 185 end %IF 186 187 % The relative error of the approximation has absolute value less 188 % than 1.e-13. One iteration of Halley's rational method (third 189 % order) gives full machine precision. 190 %k = find(0 < P & P < 1); 191 if 0 %any(k) 192 e = wnormcdf(PHINV(k)) - P(k); % error 193 u = e * sqrt(2*pi) .* exp(PHINV(k).^2/2); % f(z)/df(z) 194 % Newton's method: 195 %PHINV(k) = PHINV(k) - u; 196 % Halley's method: 197 PHINV(k) = PHINV(k) - u./( 1 + PHINV(k).*u/2 ); 198 end 199 return 200 201 202 203
Comments or corrections to the WAFO group