How to create a subroutine based on random_number that fills an array with random numbers following a normal distribution with mean mu and standard deviation sigma?
Maybe, Box-Muller transformation
Zc = sqrt(-2.0 * log(U1)) * cos(2.0*pi*U2))
Zs = sqrt(-2.0 * log(U1))* sin(2.0*pi*U2))
x = mu + sigma * Zc,s
example:
x = mu + sigma*sqrt(-2.0*log(U1))*cos(2*pi*U2)
U1,2 ~Uniform[0,1)
x ~ N[mu, sigma^2]
Reference here
Here is a code adapted from Alan Miller’s rnorm.f90 to use double precision:
module rnorm_mod
implicit none
private
public :: dp, rnorm, rnorm_vec
integer, parameter :: dp = kind(1.0d0)
contains
function rnorm_vec(n,mu,sigma) result(variates)
integer , intent(in) :: n
real(kind=dp), intent(in), optional :: mu, sigma
real(kind=dp) :: variates(n)
integer :: i
do i=1,n
variates(i) = rnorm()
end do
if (present(sigma)) variates = sigma*variates
if (present(mu)) variates = variates + mu
end function rnorm_vec
!
FUNCTION rnorm() RESULT(fn_val)
! Generate a random normal deviate using the polar method.
! Reference: Marsaglia,G. & Bray,T.A. 'A convenient method for generating
! normal variables', Siam Rev., vol.6, 260-264, 1964.
IMPLICIT NONE
REAL(kind=dp) :: fn_val
! Local variables
REAL(kind=dp) :: u, sum
REAL(kind=dp), SAVE :: v, sln
LOGICAL, SAVE :: second = .FALSE.
REAL(kind=dp), PARAMETER :: one = 1.0_dp, vsmall = TINY( one )
IF (second) THEN
! If second, use the second random number generated on last call
second = .false.
fn_val = v*sln
ELSE
! First call; generate a pair of random normals
second = .true.
DO
CALL RANDOM_NUMBER( u )
CALL RANDOM_NUMBER( v )
u = SCALE( u, 1 ) - one
v = SCALE( v, 1 ) - one
sum = u*u + v*v + vsmall ! vsmall added to prevent LOG(zero) / zero
IF(sum < one) EXIT
END DO
sln = SQRT(-SCALE(LOG(sum),1)/sum)
fn_val = u*sln
END IF
END FUNCTION rnorm
end module rnorm_mod
!
program main
use rnorm_mod, only: dp, rnorm_vec
implicit none
real(kind=dp), allocatable :: x(:)
integer, parameter :: n = 10**7
integer :: ipow, iter
real(kind=dp), parameter :: mu = 2.0_dp, sigma = 3.0_dp
call random_seed()
allocate (x(n))
print "(*(a8))", "n", "mu", "sigma"
print "(i8,2f8.4)", n , mu , sigma
print "(/,'central moments',/,*(i10))",(ipow,ipow=1,4)
do iter=1,5
x = rnorm_vec(n,mu,sigma)
x = x - sum(x)/n
print "(*(f10.4))",(sum(x**ipow)/n,ipow=1,4)
end do
print*,"theoretical"
print "(*(f10.4))",0.0_dp,sigma**2,0.0_dp,3*sigma**4
end program main
sample gfortran output:
n mu sigma
10000000 2.0000 3.0000
central moments
1 2 3 4
-0.0000 8.9980 0.0211 242.9815
0.0000 9.0006 -0.0265 243.0072
0.0000 8.9938 -0.0025 242.6678
0.0000 8.9995 -0.0236 243.0038
-0.0000 8.9914 0.0200 242.3521
theoretical
0.0000 9.0000 0.0000 243.0000
What if U1 = 0?
U ~ Uniform(0,1)
Ann. Math. Statist. 29(2): 610-611 (June, 1958). DOI: 10.1214/aoms/1177706645
1 Like
Still working on updating the NIST DATAPAC library, but it might be of interest
in particular
https://urbanjost.github.io/M_datapac/norran.3m_datapac.html
it still lacks a unit test, so caveat emptor. There are several dozen routines for various random distributions in the package.
1 Like
I updated the previous code to add the Box-Muller method and the ratio of uniforms method in addition to the polar method and to add timings. For Box-Muller I guard against zero uniform variates, which does not reduce speed significantly. For the code
module rnorm_mod
implicit none
private
public :: dp, rnorm_polar, rnorm_ratio_uni, rnorm_vec, rnorm_polar_vec, rnorm_ratio_uni_vec, &
rnorm_box_muller, rnorm_box_muller_vec
integer, parameter :: dp = kind(1.0d0)
real(kind=dp), parameter :: half = 0.5_dp, two_pi = 6.28318530718_dp
contains
function rnorm_vec(n,method,mu,sigma) result(variates)
integer , intent(in) :: n ! # of normal variates
character (len=*), intent(in) :: method
real(kind=dp), intent(in), optional :: mu ! target mean
real(kind=dp), intent(in), optional :: sigma ! target standard deviation
real(kind=dp) :: variates(n) ! normal variates
select case (method)
case ("polar") ; variates = rnorm_polar_vec(n,mu,sigma)
case ("ratio_uni") ; variates = rnorm_ratio_uni_vec(n,mu,sigma)
case ("box_muller"); variates = rnorm_box_muller_vec(n,mu,sigma)
case default ; error stop trim(method) // " is invalid"
end select
end function rnorm_vec
!
function rnorm_polar_vec(n,mu,sigma) result(variates)
integer , intent(in) :: n ! # of normal variates
real(kind=dp), intent(in), optional :: mu ! target mean
real(kind=dp), intent(in), optional :: sigma ! target standard deviation
real(kind=dp) :: variates(n) ! normal variates
integer :: i
do i=1,n
variates(i) = rnorm_polar()
end do
if (present(sigma)) variates = sigma*variates
if (present(mu)) variates = variates + mu
end function rnorm_polar_vec
!
function rnorm_ratio_uni_vec(n,mu,sigma) result(variates)
integer , intent(in) :: n ! # of normal variates
real(kind=dp), intent(in), optional :: mu ! target mean
real(kind=dp), intent(in), optional :: sigma ! target standard deviation
real(kind=dp) :: variates(n) ! normal variates
integer :: i
do i=1,n
variates(i) = rnorm_ratio_uni()
end do
if (present(sigma)) variates = sigma*variates
if (present(mu)) variates = variates + mu
end function rnorm_ratio_uni_vec
!
function rnorm_box_muller_vec(n,mu,sigma) result(variates)
integer , intent(in) :: n ! # of normal variates
real(kind=dp), intent(in), optional :: mu ! target mean
real(kind=dp), intent(in), optional :: sigma ! target standard deviation
real(kind=dp) :: variates(n) ! normal variates
integer :: i,j
logical :: n_odd
n_odd = mod(n,2) /= 0
do i=1,n/2
j = 2*i - 1
variates(j:j+1) = rnorm_box_muller()
end do
if (n_odd) variates(n) = rnorm_box_muller_single_variate()
if (present(sigma)) variates = sigma*variates
if (present(mu)) variates = variates + mu
end function rnorm_box_muller_vec
!
FUNCTION rnorm_ratio_uni() RESULT(fn_val) ! adapted from https://jblevins.org/mirror/amiller/ran_norm.f90
! Adapted from the following Fortran 77 code
! ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
! THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
! VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.
! The function random_normal() returns a normally distributed pseudo-random
! number with zero mean and unit variance.
! The algorithm uses the ratio of uniforms method of A.J. Kinderman
! and J.F. Monahan augmented with quadratic bounding curves.
real(kind=dp) :: fn_val
real(kind=dp), parameter :: s = 0.449871_dp, t = -0.386595_dp, a = 0.19600_dp, b = 0.25472_dp, &
r1 = 0.27597_dp, r2 = 0.27846_dp
! local variables
real(kind=dp) :: u, v, x, y, q
! Generate P = (u,v) uniform in rectangle enclosing acceptance region
DO
CALL RANDOM_NUMBER(u)
CALL RANDOM_NUMBER(v)
v = 1.7156_dp * (v - half)
! Evaluate the quadratic form
x = u - s
y = ABS(v) - t
q = x**2 + y*(a*y - b*x)
! Accept P if inside inner ellipse
IF (q < r1) EXIT
! Reject P if outside outer ellipse
IF (q > r2) CYCLE
! Reject P if outside acceptance region
IF (v**2 < -4.0_dp*LOG(u)*u**2) EXIT
END DO
! Return ratio of P's coordinates as the normal deviate
fn_val = v/u
END FUNCTION rnorm_ratio_uni
!
function rnorm_box_muller() result(variates) ! coded formulas from https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
! return two uncorrelated standard normal variates
real(kind=dp) :: variates(2)
real(kind=dp) :: u(2), factor, arg
do
call random_number(u)
if (u(1) > 0.0_dp) exit
end do
factor = sqrt(-2 * log(u(1)))
arg = two_pi*u(2)
variates = factor * [cos(arg),sin(arg)]
end function rnorm_box_muller
!
function rnorm_box_muller_single_variate() result(variate)
! return a standard normal variate
real(kind=dp) :: variate
real(kind=dp) :: u(2), factor, arg
call random_number(u)
factor = sqrt(-2 * log(u(1)))
arg = two_pi*u(2)
variate = factor * cos(arg)
end function rnorm_box_muller_single_variate
!
FUNCTION rnorm_polar() RESULT(fn_val) ! adapted from https://jblevins.org/mirror/amiller/rnorm.f90
! Generate a random normal deviate using the polar method.
! Reference: Marsaglia,G. & Bray,T.A. 'A convenient method for generating
! normal variables', Siam Rev., vol.6, 260-264, 1964.
IMPLICIT NONE
REAL(kind=dp) :: fn_val
! Local variables
REAL(kind=dp) :: u, sum
REAL(kind=dp), SAVE :: v, sln
LOGICAL, SAVE :: second = .FALSE.
REAL(kind=dp), PARAMETER :: one = 1.0_dp, vsmall = TINY( one )
IF (second) THEN
! If second, use the second random number generated on last call
second = .false.
fn_val = v*sln
ELSE
! First call; generate a pair of random normals
second = .true.
DO
CALL RANDOM_NUMBER( u )
CALL RANDOM_NUMBER( v )
u = SCALE( u, 1 ) - one
v = SCALE( v, 1 ) - one
sum = u*u + v*v + vsmall ! vsmall added to prevent LOG(zero) / zero
IF(sum < one) EXIT
END DO
sln = SQRT(-SCALE(LOG(sum),1)/sum)
fn_val = u*sln
END IF
END FUNCTION rnorm_polar
end module rnorm_mod
!
program main
use rnorm_mod, only: dp, rnorm_vec, rnorm_box_muller_vec
implicit none
integer , parameter :: n = 10**8, niter = 5, nmethods = 3
real(kind=dp) , parameter :: mu = 2.0_dp, sigma = 3.0_dp
character (len=*), parameter :: methods(nmethods) = ["polar ", &
"ratio_uni ","box_muller"]
integer :: ipow, iter, imethod
real(kind=dp) :: xmean, old_time, new_time, elapsed_time(nmethods)
real(kind=dp), allocatable :: x(:)
call random_seed()
elapsed_time = 0.0_dp
allocate (x(n))
print "(*(a10))", "n", "niter", "mu", "sigma"
print "(2i10,2f10.4)", n , niter, mu , sigma
do imethod=1,nmethods
print "(/,'method = ',a)",trim(methods(imethod))
print "('central moments',/,a10,*(i10))","mean",(ipow,ipow=1,4)
do iter=1,niter
call cpu_time(old_time)
x = rnorm_vec(n,methods(imethod),mu,sigma)
call cpu_time(new_time)
elapsed_time(imethod) = elapsed_time(imethod) + new_time - old_time
xmean = sum(x)/n
x = x - xmean
print "(*(f10.4))",xmean,(sum(x**ipow)/n,ipow=1,4)
end do
print*,"theoretical"
print "(*(f10.4))",mu,0.0_dp,sigma**2,0.0_dp,3*sigma**4
end do
print "(/,*(a15))", "method",(trim(methods(imethod)),imethod=1,nmethods)
print "(a15,*(f15.4))", "cpu_time",elapsed_time
print "(/,'check Box-Muller for n odd or even')"
print "(*(f10.4))",rnorm_box_muller_vec(3),rnorm_box_muller_vec(4)
end program main
on WSL2 with gfortran -O3 I get timings
method polar ratio_uni box_muller
cpu_time 13.0271 15.3004 12.3667
and with ifort -O3
method polar ratio_uni box_muller
cpu_time 16.8491 20.3027 8.6018
Full output
gfortran -O3 rnorm_dp.f90
(base) /mnt/c/fortran/test$ time ./a.out
n niter mu sigma
100000000 5 2.0000 3.0000
method = polar
central moments
mean 1 2 3 4
2.0002 0.0000 9.0006 -0.0003 243.0160
2.0001 -0.0000 8.9974 0.0124 242.8701
2.0000 -0.0000 9.0003 -0.0028 242.9531
2.0001 0.0000 8.9992 -0.0094 242.9653
1.9999 0.0000 9.0008 -0.0128 243.0604
theoretical
2.0000 0.0000 9.0000 0.0000 243.0000
method = ratio_uni
central moments
mean 1 2 3 4
1.9999 -0.0000 9.0011 -0.0086 243.0712
2.0003 0.0000 8.9992 -0.0098 242.9281
1.9999 -0.0000 9.0001 0.0071 243.0048
1.9998 -0.0000 8.9997 -0.0047 242.9843
2.0002 -0.0000 9.0011 -0.0098 243.0690
theoretical
2.0000 0.0000 9.0000 0.0000 243.0000
method = box_muller
central moments
mean 1 2 3 4
2.0000 -0.0000 9.0007 0.0096 243.1083
1.9992 -0.0000 9.0012 0.0070 243.0418
2.0000 0.0000 9.0015 0.0044 243.0983
2.0000 -0.0000 8.9992 0.0017 242.9084
2.0005 0.0000 9.0005 -0.0068 243.0650
theoretical
2.0000 0.0000 9.0000 0.0000 243.0000
method polar ratio_uni box_muller
cpu_time 13.0271 15.3004 12.3667
check Box-Muller for n odd or even
-0.4509 0.6263 0.1658 -0.5268 0.2805 1.0316 -0.5751
real 0m58.591s
user 0m58.355s
sys 0m0.230s
ulimit -s unlimited
(base) /mnt/c/fortran/test$ ifort -O3 rnorm_dp.f90
(base) /mnt/c/fortran/test$ time ./a.out
n niter mu sigma
100000000 5 2.0000 3.0000
method = polar
central moments
mean 1 2 3 4
2.0002 -0.0000 8.9993 0.0057 242.9807
2.0002 0.0000 8.9988 -0.0029 242.8929
2.0006 0.0000 8.9994 0.0061 243.0622
2.0007 -0.0000 9.0016 -0.0039 243.1514
2.0002 0.0000 8.9988 -0.0034 242.8767
theoretical
2.0000 0.0000 9.0000 0.0000 243.0000
method = ratio_uni
central moments
mean 1 2 3 4
2.0001 -0.0000 9.0002 0.0071 242.9891
2.0000 0.0000 9.0012 0.0030 243.0431
2.0003 -0.0000 8.9997 0.0047 242.9036
1.9995 0.0000 9.0012 0.0052 243.0759
1.9997 -0.0000 9.0008 0.0046 242.9873
theoretical
2.0000 0.0000 9.0000 0.0000 243.0000
method = box_muller
central moments
mean 1 2 3 4
2.0002 -0.0000 8.9980 -0.0033 242.9163
2.0003 0.0000 8.9983 -0.0006 242.8713
1.9994 0.0000 8.9998 -0.0070 243.0213
1.9998 0.0000 9.0010 0.0029 243.1056
2.0002 0.0000 8.9998 -0.0025 243.0385
theoretical
2.0000 0.0000 9.0000 0.0000 243.0000
method polar ratio_uni box_muller
cpu_time 16.8491 20.3027 8.6018
check Box-Muller for n odd or even
-0.4093 -0.5503 -0.5804 0.6342 0.7505 -1.8514 -0.9049
real 1m30.322s
user 1m30.162s
sys 0m0.141s
2 Likes
One way to get normal deviates is to use the inverse error function, see here: Inverse error function - MATLAB erfinv.
This method is also discussed in the following paper:
I didn’t know that Fortran has an erfinv function. Can you be more specific? I only know how to use the erfc function to get normcdf.
The suggested function is not a standard intrinsic function, but a user-provided function. One of the references provides approximations for erfinv() that can be coded fairly easily.
1 Like
The method from Giles I linked above has been designed to run fast on NVIDIA GPU’s. If I’m not mistaken, the author, Mike Giles, is somehow associated with NAG. The paper also states:
These approximations were originally developed as part of a commercial maths library providing Numerical Routines for GPUs.
There exist also older approximations, that perform well on CPU’s, such as
Blair, J. M., Edwards, C. A., & Johnson, J. H. (1976). Rational Chebyshev approximations for the inverse of the error function. Mathematics of computation, 30(136), 827-830.
A copy of that paper can be found here, but the tables of polynomial coefficients are missing. You can however orient yourself by the Julia implementantion: https://github.com/JuliaMath/SpecialFunctions.jl/blob/49be7e6add22057776b5ed8ae32e0dae84500203/src/erf.jl#L242
A draft Fortran attempt can be found here: List of special mathematical functions to include · Issue #179 · fortran-lang/stdlib · GitHub. Unfortunately, I haven’t done any proper testing of the accuracy or speed. The Julia special functions library provides both single and double precision versions.
I’ve located the routines from Giles in a third-party GitHub repository: lwafomps/erfinv.cpp at 1ff3fbbdb23d5772c1c45d201ba050e4819a23e8 · MersenneTwister-Lab/lwafomps · GitHub. A few (related) C codes and discussion can also be found in this thread: Inverse Error Function in C - Stack Overflow
I would also note v?ErfInv in the Intel oneAPI MKL Vector Mathematical Functions library.
Speaking of Intel MKL, the Statistical Functions include distribution generators including v?RngGaussian. Methods based on the Box-Muller transform or the inverse cumulative distribution function are both available.
If you don’t mind linking with the C++ standard library, there you have: normal_distribution - C++ Reference
Finally, if you’re using R, you could also use R’s internal pseudo random number generator as described in Writing R Extensions. The interfaces needed to call these from Fortran are very simple:
module r_random
implicit none
interface
subroutine r_getrngstate() bind(c,name="GetRNGstate")
end subroutine
subroutine r_putrngstate() bind(c,name="PutRNGstate")
end subroutine
real(c_double) function r_unif_rand() bind(c,name="unif_rand")
use, intrinsic :: iso_c_binding, only: c_double
end function
real(c_double) function r_norm_rand() bind(c,name="norm_rand")
use, intrinsic :: iso_c_binding, only: c_double
end function
real(c_double) function r_exp_rand() bind(c,name="exp_rand")
use, intrinsic :: iso_c_binding, only: c_double
end function
end interface
end module
You still have to figure out how to link your Fortran code with R’s standalone math library.
1 Like
I also use the Box-Muller method in my codes, but wondering if there are better (= more efficient) methods. In addition to the methods mentioned above, is the “ziggurat” method another option? (I’ve never used it myself, so not sure whether it is good or not…) GSL also seems to provide both routines:
https://www.gnu.org/software/gsl/doc/html/randist.html#the-gaussian-distribution
double gsl_ran_gaussian_ziggurat(const gsl_rng *r, double sigma)
This function computes a Gaussian random variate using
the alternative Marsaglia-Tsang ziggurat and Kinderman-Monahan-
Leva ratio methods. The Ziggurat algorithm is the fastest
available algorithm in most cases.My adaptation of Alan Miller’s ziggurat Fortran code is discussed at Pure random number generators and is at GitHub - Beliavsky/Ziggurat: generate uniform, normal, and exponential random deviates using algorithms from Marsaglia.
1 Like
I’ve found that both Box-Muller and Ziggurat methods are discussed in the book Random Number Generators – Principles and Practices (2018) by David Johnston (Principal Engineer, Intel Corporation).
The NAG Library routine g05skf (dist_normal) uses the inverse cumulative distribution function approach with an approximation by Wichura (1988). A PDF is available here (511 kb).
Searching through ACM publications shows lots of results:
- Algorithm 488: A Gaussian pseudo-random number generator (1974)
- Algorithm 712: A normal random number generator (1992)
- Fast pseudorandom generators for normal and exponential variates (1996)
- An Improved Ziggurat Method to Generate Normal Random Samples (2005)
- Gaussian random number generators (2007)
Some comments on Ziggurat can also be found in the stdlib thread: Proposal for statistical distributions · Issue #234 · fortran-lang/stdlib · GitHub. Another discussion is in the NumPy pull request: Patch with Ziggurat method for Normal distribution (Trac #1450) · Issue #2047 · numpy/numpy · GitHub
I’m surprised that we’ve all forgotten Fortran stdlib already has a function for normal variates: stats_distribution_normal – Fortran-lang/stdlib
program demo_normal_rvs
use stdlib_random, only: random_seed
use stdlib_stats_distribution_normal, only: norm => rvs_normal
implicit none
integer :: seed_put, seed_get
seed_put = 1234567
call random_seed(seed_put, seed_get)
print *, norm( ) ! single standard normal random variate
print *, norm(1.0, 2.0) ! normal random variate mu=1.0, sigma=2.0
print *, norm(0.0, 1.0, 10) ! an array of 10 standard norml random variates
end program
1 Like
Most algorithms for generating normal variates rely on uniform variates. It is possible to write the module with procedures for normal variates so that toggling a single line of code determines whether the intrinsic or user-defined random_number is used, as shown below. A benefit of the latter is that the normal variates will be the same across compilers.
module random_number_mod
implicit none
private
public :: dp, set_seed, random_number
integer, parameter :: dp = kind(1.0d0)
real(kind=dp), save :: ds = 123.0_dp
contains
!
subroutine set_seed(seed)
real(kind=dp), intent(in) :: seed
ds = seed
end subroutine set_seed
!
impure elemental subroutine random_number(harvest)
! adapted from function ggl in http://www.netlib.org/random/2drwtest.f
real(kind=dp), intent(out) :: harvest
real(kind=dp), parameter :: d2 = 2147483647.d0
ds = dmod(16807.d0*ds,d2)
! Generate U(0,1] distributed random numbers:
harvest = ds/d2
end subroutine random_number
end module random_number_mod
!
module random_normal_mod
! comment the line below to use the intrinsic random_number
use random_number_mod, only: random_number
implicit none
private
public :: dp, rnorm_polar, rnorm_ratio_uni, rnorm_vec, rnorm_polar_vec, rnorm_ratio_uni_vec, &
rnorm_box_muller, rnorm_box_muller_vec
integer, parameter :: dp = kind(1.0d0)
real(kind=dp), parameter :: half = 0.5_dp, two_pi = 6.28318530718_dp
contains
function rnorm_vec(n,method,mu,sigma) result(variates)
integer , intent(in) :: n ! # of normal variates
character (len=*), intent(in) :: method
real(kind=dp), intent(in), optional :: mu ! target mean
real(kind=dp), intent(in), optional :: sigma ! target standard deviation
real(kind=dp) :: variates(n) ! normal variates
select case (method)
case ("polar") ; variates = rnorm_polar_vec(n,mu,sigma)
case ("ratio_uni") ; variates = rnorm_ratio_uni_vec(n,mu,sigma)
case ("box_muller"); variates = rnorm_box_muller_vec(n,mu,sigma)
case default ; error stop trim(method) // " is invalid"
end select
end function rnorm_vec
!
function rnorm_polar_vec(n,mu,sigma) result(variates)
integer , intent(in) :: n ! # of normal variates
real(kind=dp), intent(in), optional :: mu ! target mean
real(kind=dp), intent(in), optional :: sigma ! target standard deviation
real(kind=dp) :: variates(n) ! normal variates
integer :: i
do i=1,n
variates(i) = rnorm_polar()
end do
if (present(sigma)) variates = sigma*variates
if (present(mu)) variates = variates + mu
end function rnorm_polar_vec
!
function rnorm_ratio_uni_vec(n,mu,sigma) result(variates)
integer , intent(in) :: n ! # of normal variates
real(kind=dp), intent(in), optional :: mu ! target mean
real(kind=dp), intent(in), optional :: sigma ! target standard deviation
real(kind=dp) :: variates(n) ! normal variates
integer :: i
do i=1,n
variates(i) = rnorm_ratio_uni()
end do
if (present(sigma)) variates = sigma*variates
if (present(mu)) variates = variates + mu
end function rnorm_ratio_uni_vec
!
function rnorm_box_muller_vec(n,mu,sigma) result(variates)
integer , intent(in) :: n ! # of normal variates
real(kind=dp), intent(in), optional :: mu ! target mean
real(kind=dp), intent(in), optional :: sigma ! target standard deviation
real(kind=dp) :: variates(n) ! normal variates
integer :: i,j
logical :: n_odd
n_odd = mod(n,2) /= 0
do i=1,n/2
j = 2*i - 1
variates(j:j+1) = rnorm_box_muller()
end do
if (n_odd) variates(n) = rnorm_box_muller_single_variate()
if (present(sigma)) variates = sigma*variates
if (present(mu)) variates = variates + mu
end function rnorm_box_muller_vec
!
FUNCTION rnorm_ratio_uni() RESULT(fn_val) ! adapted from https://jblevins.org/mirror/amiller/ran_norm.f90
! Adapted from the following Fortran 77 code
! ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
! THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
! VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.
! The function random_normal() returns a normally distributed pseudo-random
! number with zero mean and unit variance.
! The algorithm uses the ratio of uniforms method of A.J. Kinderman
! and J.F. Monahan augmented with quadratic bounding curves.
real(kind=dp) :: fn_val
real(kind=dp), parameter :: s = 0.449871_dp, t = -0.386595_dp, a = 0.19600_dp, b = 0.25472_dp, &
r1 = 0.27597_dp, r2 = 0.27846_dp
! local variables
real(kind=dp) :: u, v, x, y, q
! Generate P = (u,v) uniform in rectangle enclosing acceptance region
DO
CALL RANDOM_NUMBER(u)
CALL RANDOM_NUMBER(v)
v = 1.7156_dp * (v - half)
! Evaluate the quadratic form
x = u - s
y = ABS(v) - t
q = x**2 + y*(a*y - b*x)
! Accept P if inside inner ellipse
IF (q < r1) EXIT
! Reject P if outside outer ellipse
IF (q > r2) CYCLE
! Reject P if outside acceptance region
IF (v**2 < -4.0_dp*LOG(u)*u**2) EXIT
END DO
! Return ratio of P's coordinates as the normal deviate
fn_val = v/u
END FUNCTION rnorm_ratio_uni
!
function rnorm_box_muller() result(variates) ! coded formulas from https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
! return two uncorrelated standard normal variates
real(kind=dp) :: variates(2)
real(kind=dp) :: u(2), factor, arg
do
call random_number(u)
if (u(1) > 0.0_dp) exit
end do
factor = sqrt(-2 * log(u(1)))
arg = two_pi*u(2)
variates = factor * [cos(arg),sin(arg)]
end function rnorm_box_muller
!
function rnorm_box_muller_single_variate() result(variate)
! return a standard normal variate
real(kind=dp) :: variate
real(kind=dp) :: u(2), factor, arg
do
call random_number(u)
if (u(1) > 0.0_dp) exit
end do
factor = sqrt(-2 * log(u(1)))
arg = two_pi*u(2)
variate = factor * cos(arg)
end function rnorm_box_muller_single_variate
!
FUNCTION rnorm_polar() RESULT(fn_val) ! adapted from https://jblevins.org/mirror/amiller/rnorm.f90
! Generate a random normal deviate using the polar method.
! Reference: Marsaglia,G. & Bray,T.A. 'A convenient method for generating
! normal variables', Siam Rev., vol.6, 260-264, 1964.
IMPLICIT NONE
REAL(kind=dp) :: fn_val
! Local variables
REAL(kind=dp) :: u, sum
REAL(kind=dp), SAVE :: v, sln
LOGICAL, SAVE :: second = .FALSE.
REAL(kind=dp), PARAMETER :: one = 1.0_dp, vsmall = TINY( one )
IF (second) THEN
! If second, use the second random number generated on last call
second = .false.
fn_val = v*sln
ELSE
! First call; generate a pair of random normals
second = .true.
DO
CALL RANDOM_NUMBER( u )
CALL RANDOM_NUMBER( v )
u = SCALE( u, 1 ) - one
v = SCALE( v, 1 ) - one
sum = u*u + v*v + vsmall ! vsmall added to prevent LOG(zero) / zero
IF(sum < one) EXIT
END DO
sln = SQRT(-SCALE(LOG(sum),1)/sum)
fn_val = u*sln
END IF
END FUNCTION rnorm_polar
end module random_normal_mod
!
program main
use random_normal_mod, only: dp, rnorm_vec, rnorm_box_muller_vec
implicit none
integer , parameter :: n = 10**6, niter = 5, nmethods = 3
real(kind=dp) , parameter :: mu = 2.0_dp, sigma = 3.0_dp
character (len=*), parameter :: methods(nmethods) = ["polar ", &
"ratio_uni ","box_muller"]
integer :: ipow, iter, imethod
real(kind=dp) :: xmean, old_time, new_time, elapsed_time(nmethods)
real(kind=dp), allocatable :: x(:)
call random_seed()
elapsed_time = 0.0_dp
allocate (x(n))
print "(*(a10))", "n", "niter", "mu", "sigma"
print "(2i10,2f10.4)", n , niter, mu , sigma
do imethod=1,nmethods
print "(/,'method = ',a)",trim(methods(imethod))
print "('central moments',/,a10,*(i10))","mean",(ipow,ipow=1,4)
do iter=1,niter
call cpu_time(old_time)
x = rnorm_vec(n,methods(imethod),mu,sigma)
call cpu_time(new_time)
elapsed_time(imethod) = elapsed_time(imethod) + new_time - old_time
xmean = sum(x)/n
x = x - xmean
print "(*(f10.4))",xmean,(sum(x**ipow)/n,ipow=1,4)
end do
print*,"theoretical"
print "(*(f10.4))",mu,0.0_dp,sigma**2,0.0_dp,3*sigma**4
end do
print "(/,*(a15))", "method",(trim(methods(imethod)),imethod=1,nmethods)
print "(a15,*(f15.4))", "cpu_time",elapsed_time
print "(/,'check Box-Muller for n odd or even')"
print "(*(f10.4))",rnorm_box_muller_vec(3),rnorm_box_muller_vec(4)
end program main
subroutine rnorm(x, mu, sigma)
integer(8) :: n
real(8) :: x(:), mu, sigma
real(8), allocatable :: u1(:), u2(:)
n = size(x, kind = 8)
allocate(u1(n), u2(n))
call random_number(u1)
call random_number(u2)
u1 = 1d0 - u1
u2 = 1d0 - u2
x = mu + sigma*sqrt(-2d0*log(u1))*cos(2d0*pi*u2)
deallocate(u1, u2)
end subroutine
end
I followed your reference to write this subroutine, but it is slow. Any advice on how to make it faster?
My function rnorm_box_muller generates two normal variates rather than one without doing twice as many computations, by computing a sin as well as a cos.
Many math libraries have a subroutine that computes both sine and cosine values of the argument. In addition to the range reduction, some of the polynomial evaluation steps can be shared, making the total effort less than the sum of the effort for the two separate evaluations.
I wonder what is the typical value of n (the size of array x)? If it is very large (such that it requires 64-bit integer for the size), I guess allocating the two temporary arrays u1 and u2 may be very time-consuming…