function p = phi(z) % % normal distribution probabilities accurate to 1.e-15. % z = no. of standard deviations from the mean. % % based upon algorithm 5666 for the error function, from: % Hart, J.F. et al, 'Computer Approximations', Wiley 1968 % % programmer: Alan Miller % % latest revision - 30 march 1986 % % modified for matlab by Alan Genz - 7/98 % p0 = 220.2068679123761; p1 = 221.2135961699311; p2 = 112.0792914978709; p3 = 33.91286607838300; p4 = 6.373962203531650; p5 = .7003830644436881; p6 = .03526249659989109; q0 = 440.4137358247522; q1 = 793.8265125199484; q2 = 637.3336333788311; q3 = 296.5642487796737; q4 = 86.78073220294608; q5 = 16.06417757920695; q6 = 1.755667163182642; q7 = .08838834764831844; rootpi = 2.506628274631001; cutoff = 7.071067811865475; zabs = abs(z); if ( zabs > 37 ) % % % |z| > 37 % p = 0; else % % |z| <= 37; % expntl = exp(-zabs^2/2); if ( zabs < cutoff ) % % |z| < cutoff = 10/sqrt(2); % pt = ((((((p6*zabs+p5)*zabs+p4)*zabs+p3)*zabs+p2)*zabs+p1)*zabs+p0); qt = ((((((q7*zabs+q6)*zabs+q5)*zabs+q4)*zabs+q3)*zabs+q2)*zabs+q1); p = expntl*pt/( qt*zabs + q0 ); else % % |z| >= cutoff % p = expntl/(zabs+1/(zabs+2/(zabs+3/(zabs+4/(zabs+0.65)))))/rootpi; end end if z > 0 , p = 1 - p; end % end phi