#! /usr/bin/env python # def normal_01_cdf ( x ): #*****************************************************************************80 # ## NORMAL_01_CDF evaluates the Normal 01 CDF. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 27 February 2015 # # Author: # # John Burkardt # # Reference: # # A G Adams, # Areas Under the Normal Curve, # Algorithm 39, # Computer j., # Volume 12, pages 197-198, 1969. # # Parameters: # # Input, real X, the argument of the CDF. # # Output, real VALUE, the value of the CDF. # import numpy as np a1 = 0.398942280444 a2 = 0.399903438504 a3 = 5.75885480458 a4 = 29.8213557808 a5 = 2.62433121679 a6 = 48.6959930692 a7 = 5.92885724438 b0 = 0.398942280385 b1 = 3.8052E-08 b2 = 1.00000615302 b3 = 3.98064794E-04 b4 = 1.98615381364 b5 = 0.151679116635 b6 = 5.29330324926 b7 = 4.8385912808 b8 = 15.1508972451 b9 = 0.742380924027 b10 = 30.789933034 b11 = 3.99019417011 # # |X| <= 1.28. # if ( abs ( x ) <= 1.28 ): y = 0.5 * x * x q = 0.5 - abs ( x ) * ( a1 - a2 * y / ( y + a3 \ - a4 / ( y + a5 \ + a6 / ( y + a7 ) ) ) ) # # 1.28 < |X| <= 12.7 # elif ( abs ( x ) <= 12.7 ): y = 0.5 * x * x q = np.exp ( - y ) \ * b0 / ( abs ( x ) - b1 \ + b2 / ( abs ( x ) + b3 \ + b4 / ( abs ( x ) - b5 \ + b6 / ( abs ( x ) + b7 \ - b8 / ( abs ( x ) + b9 \ + b10 / ( abs ( x ) + b11 ) ) ) ) ) ) # # 12.7 < |X| # else: q = 0.0 # # Take account of negative X. # if ( x < 0.0 ): value = q else: value = 1.0 - q return value def normal_01_cdf_test ( ): #*****************************************************************************80 # ## NORMAL_01_CDF_TEST tests NORMAL_01_CDF. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 27 February 2015 # # Author: # # John Burkardt # import platform print ( '' ) print ( 'NORMAL_01_CDF_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' NORMAL_01_CDF evaluates the CDF;' ) print ( '' ) print ( ' X CDF CDF' ) print ( ' (exact) (computed)' ) print ( '' ) n_data = 0 while ( True ): n_data, x, cdf1 = normal_01_cdf_values ( n_data ) if ( n_data == 0 ): break cdf2 = normal_01_cdf ( x ) print ( ' %14.6g %24.16g %24.16g' % ( x, cdf1, cdf2 ) ) # # Terminate. # print ( '' ) print ( 'NORMAL_01_CDF_TEST:' ) print ( ' Normal end of execution.' ) return def normal_01_cdf_values ( n_data ): #*****************************************************************************80 # ## NORMAL_01_CDF_VALUES returns some values of the Normal 01 CDF. # # Discussion: # # In Mathematica, the function can be evaluated by: # # Needs["Statistics`ContinuousDistributions`"] # dist = NormalDistribution [ 0, 1 ] # CDF [ dist, x ] # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 14 February 2015 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz and Irene Stegun, # Handbook of Mathematical Functions, # US Department of Commerce, 1964. # # Stephen Wolfram, # The Mathematica Book, # Fourth Edition, # Wolfram Media / Cambridge University Press, 1999. # # Parameters: # # Input/output, integer N_DATA. The user sets N_DATA to 0 before the # first call. On each call, the routine increments N_DATA by 1, and # returns the corresponding data; when there is no more data, the # output value of N_DATA will be 0 again. # # Output, real X, the argument of the function. # # Output, real F, the value of the function. # import numpy as np n_max = 17 f_vec = np.array ( (\ 0.5000000000000000E+00, \ 0.5398278372770290E+00, \ 0.5792597094391030E+00, \ 0.6179114221889526E+00, \ 0.6554217416103242E+00, \ 0.6914624612740131E+00, \ 0.7257468822499270E+00, \ 0.7580363477769270E+00, \ 0.7881446014166033E+00, \ 0.8159398746532405E+00, \ 0.8413447460685429E+00, \ 0.9331927987311419E+00, \ 0.9772498680518208E+00, \ 0.9937903346742239E+00, \ 0.9986501019683699E+00, \ 0.9997673709209645E+00, \ 0.9999683287581669E+00 )) x_vec = np.array ((\ 0.0000000000000000E+00, \ 0.1000000000000000E+00, \ 0.2000000000000000E+00, \ 0.3000000000000000E+00, \ 0.4000000000000000E+00, \ 0.5000000000000000E+00, \ 0.6000000000000000E+00, \ 0.7000000000000000E+00, \ 0.8000000000000000E+00, \ 0.9000000000000000E+00, \ 0.1000000000000000E+01, \ 0.1500000000000000E+01, \ 0.2000000000000000E+01, \ 0.2500000000000000E+01, \ 0.3000000000000000E+01, \ 0.3500000000000000E+01, \ 0.4000000000000000E+01 )) if ( n_data < 0 ): n_data = 0 if ( n_max <= n_data ): n_data = 0 x = 0.0 f = 0.0 else: x = x_vec[n_data] f = f_vec[n_data] n_data = n_data + 1 return n_data, x, f if ( __name__ == '__main__' ): from timestamp import timestamp timestamp ( ) normal_01_cdf_test ( ) timestamp ( )