Home > wafo > wstats > wnorminv.m

wnorminv

PURPOSE ^

Inverse of the Normal distribution function

SYNOPSIS ^

x = wnorminv(F,m,v)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

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