#! /usr/bin/env python # def snorm ( ): #*****************************************************************************80 # ## SNORM samples the standard normal distribution. # # Discussion: # # This procedure corresponds to algorithm FL, with M = 5, in the reference. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 02 September 2018 # # Author: # # Original FORTRAN77 version by Barry Brown, James Lovato. # Python version by John Burkardt. # # Reference: # # Joachim Ahrens, Ulrich Dieter, # Extensions of Forsythe's Method for Random # Sampling from the Normal Distribution, # Mathematics of Computation, # Volume 27, Number 124, October 1973, page 927-937. # # Parameters: # # Output, real VALUE, a random deviate from the distribution. # import numpy as np from r4_uni_01 import r4_uni_01 a = np.array ( [ \ 0.0000000, 0.3917609E-01, 0.7841241E-01, 0.1177699, \ 0.1573107, 0.1970991, 0.2372021, 0.2776904, \ 0.3186394, 0.3601299, 0.4022501, 0.4450965, \ 0.4887764, 0.5334097, 0.5791322, 0.6260990, \ 0.6744898, 0.7245144, 0.7764218, 0.8305109, \ 0.8871466, 0.9467818, 1.009990, 1.077516, \ 1.150349, 1.229859, 1.318011, 1.417797, \ 1.534121, 1.675940, 1.862732, 2.153875 ] ) d = np.array ( [ \ 0.0000000, 0.0000000, 0.0000000, 0.0000000, \ 0.0000000, 0.2636843, 0.2425085, 0.2255674, \ 0.2116342, 0.1999243, 0.1899108, 0.1812252, \ 0.1736014, 0.1668419, 0.1607967, 0.1553497, \ 0.1504094, 0.1459026, 0.1417700, 0.1379632, \ 0.1344418, 0.1311722, 0.1281260, 0.1252791, \ 0.1226109, 0.1201036, 0.1177417, 0.1155119, \ 0.1134023, 0.1114027, 0.1095039 ] ) h = np.array ( [ \ 0.3920617E-01, 0.3932705E-01, 0.3950999E-01, 0.3975703E-01, \ 0.4007093E-01, 0.4045533E-01, 0.4091481E-01, 0.4145507E-01, \ 0.4208311E-01, 0.4280748E-01, 0.4363863E-01, 0.4458932E-01, \ 0.4567523E-01, 0.4691571E-01, 0.4833487E-01, 0.4996298E-01, \ 0.5183859E-01, 0.5401138E-01, 0.5654656E-01, 0.5953130E-01, \ 0.6308489E-01, 0.6737503E-01, 0.7264544E-01, 0.7926471E-01, \ 0.8781922E-01, 0.9930398E-01, 0.1155599, 0.1404344, \ 0.1836142, 0.2790016, 0.7010474 ] ) t = np.array ( [ \ 0.7673828E-03, 0.2306870E-02, 0.3860618E-02, 0.5438454E-02, \ 0.7050699E-02, 0.8708396E-02, 0.1042357E-01, 0.1220953E-01, \ 0.1408125E-01, 0.1605579E-01, 0.1815290E-01, 0.2039573E-01, \ 0.2281177E-01, 0.2543407E-01, 0.2830296E-01, 0.3146822E-01, \ 0.3499233E-01, 0.3895483E-01, 0.4345878E-01, 0.4864035E-01, \ 0.5468334E-01, 0.6184222E-01, 0.7047983E-01, 0.8113195E-01, \ 0.9462444E-01, 0.1123001, 0.1364980, 0.1716886, \ 0.2276241, 0.3304980, 0.5847031 ] ) u = r4_uni_01 ( ) if ( u <= 0.5 ): s = 0.0 else: s = 1.0 u = 2.0 * u - s u = 32.0 * u i = int ( u ) if ( i == 32 ): i = 31 # # Center # if ( i != 0 ): ustar = u - i aa = a[i-1] while ( True ): if ( t[i-1] < ustar ): w = ( ustar - t[i-1] ) * h[i-1] y = aa + w if ( s != 1.0 ): value = y else: value = -y return value u = r4_uni_01 ( ) w = u * ( a[i] - aa ) tt = ( 0.5 * w + aa ) * w while ( True ): if ( tt < ustar ): y = aa + w if ( s != 1.0 ): value = y else: value = - y return value u = r4_uni_01 ( ) if ( ustar < u ): break tt = u ustar = r4_uni_01 ( ) ustar = r4_uni_01 ( ) # # Tail # else: i = 6 aa = a[31] while ( True ): u = u + u if ( 1.0 <= u ): break aa = aa + d[i-1] i = i + 1 u = u - 1.0 w = u * d[i-1] tt = ( 0.5 * w + aa ) * w while ( True ): ustar = r4_uni_01 ( ) if ( tt < ustar ): y = aa + w if ( s != 1.0 ): value = y else: value = -y return value u = r4_uni_01 ( ) if ( u <= ustar ): tt = u else: u = r4_uni_01 ( ) w = u * d[i-1] tt = ( 0.5 * w + aa ) * w return value def snorm_test ( ): #*****************************************************************************80 # ## SNORM_TEST tests SNORM. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 02 September 2018 # # Author: # # John Burkardt # import platform test_num = 20 print ( '' ) print ( 'SNORM_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' SNORM generates normally distributed' ) print ( ' random values.' ) print ( '' ) for test in range ( 0, test_num ): x = snorm ( ) print ( ' %f' % ( x ) ) # # Terminate. # print ( '' ) print ( 'SNORM_TEST' ) print ( ' Normal end of execution.' ) return if ( __name__ == '__main__' ): from timestamp import timestamp timestamp ( ) snorm_test ( ) timestamp ( )