Random numbers

Yutaka Masuda

February 2020

Back to index.html.

Random numbers

Random numbers in Fortran

Random numbers are useful in Monte Carlo studies. The numbers are generated from a definite formula, so the sequence of numbers is not really random but nearly random (i.e., pseudo-random numbers). The formula is called (pseudo) random number generator.

Fortran has an intrinsic subroutine (random_number) to generate random numbers. It generates a real number ranging from 0 to 1 (i.e., following a uniform distribution). The implementation of the generator depends on a compiler; the other compiler should generate a different sequence! The built-in generator is not useful if you need more complex random numbers (e.g., normal deviate) or a portable generator. Therefore, in such a case, you have to find external functions to do it, or you should program it by yourself. I do not recommend writing the random number generators because a random number generator is (very) complicated.

In this chapter, first, I show how the intrinsic generator works, and then, I introduce a few algorithms to generate random numbers following statistical distributions. The aim of my code shown here is to demonstrate simple random-number generators, and you should not use them in your serious study. Please use solid generators available online for the purpose.

Intrinsic random number generator

Basic usage

An intrinsic subroutine, random_number, generates a real random-number \(r\) uniformly distributed in \(0\leq r<1\). The subroutine returns the same type of real value as the argument.

! a variable to receive a random number
real :: r
double precision :: s
...
call random_number(r)  ! for single precision
call random_number(s)  ! for double precision

If you give an array as an argument, the subroutine fills all the elements with the random numbers.

real :: r(5),v(10,10)

call random_number(r)  ! returns 5 numbers at the same time
call random_number(v)  ! returns 100 numbers at the same time

Note that the generated random-numbers depend on the compiler (and often on the computer that the program is running). Do not expect that this subroutine gives you a unique sequence of random numbers across environments.

Seeding a random number

Given a definite environment, random_number generates a defined sequence of random numbers. The subroutine keeps an internal state, which is the source of random numbers. Once you call random_number, it changes the internal state and generates a new random number. It means, when you start from a particular state, you should get the same sequence of random numbers. This reproducibility is useful to reproduce the same results in simulation.

The seed is integer numbers that define the initial state. The seed size depends on random number generators, and each compiler implements a different generator in random_number. It means that some compilers define it as an integer vector with 2 elements, but some may use a vector with 4 elements. So, in Fortran, a subroutine (random_seed) tells you how many integers are needed in the seed vector. Also, this subroutine puts a new seed and returns the current seed.

Here is the following syntax to inquire the number of elements in a seed vector. You need a variable to receive the number.

integer :: n

call random_seed(size=n)

An idiom to prepare the seed is to allocate a seed vector and to put arbitrary values to the vector. The following code puts the seed.

integer :: n
integer,allocatable :: seed(:)

call random_seed(size=n)
allocate(seed(n))
seed = 123456789    ! putting arbitrary seed to all elements
call random_seed(put=seed)
deallocate(seed)

If you want to know the current seed, option get= is useful.

integer :: n
integer,allocatable :: seed(:)

call random_seed(size=n)
allocate(seed(n))
call random_seed(get=seed)
print *,seed
dealloocate(seed)

Exercises

  1. You should hit a runtime-error when computing \(\log(r)\) with a random number \(r\). Explain why it should occur, and suggest a solution if any.
  2. Confirm that you get the same sequence of random numbers if you put the same seed.
  3. Make a subroutine to put a scholar to the seed vector.

Random numbers from particular distributions

In this section, I introduce some algorithms to generate random numbers following the uniform, normal, and gamma distributions. There are many books and Internet articles about random numbers; the following is based on Statistical Computing by Kennedy and Gentle.

Note that the algorithms should not be efficient but be simple. If you need high-performance random number generators, please use optimized subroutines prepared by experts.

Standard uniform function

Note that some algorithms use \(\ln[U(0,1)]\), but as seen in the exercise above, it fails if it contains 0. Hereafter, in a mathematical notation, we assume that \(u\sim U(0,1)\) where \(0<u\leq 1\). It is available by transforming a random number \(0\leq r<1\) as \(u=1-r\). First, we develop a subroutine to return \(u\).

subroutine random_stduniform(u)
   implicit none
   real,intent(out) :: u
   real :: r
   call random_number(r)
   u = 1 - r
end subroutine random_stduniform

Uniform distribution

Consider the uniform distribution \(x \sim U(a,b)\) where \(a<x\leq b\). The random number is available by a linear transformation of \(u\sim U(0,1)\), which is the standard uniform distribution implemented in random_number.

\[ x = (b-a)u + a \]

We can implement it in a subroutine.

! assuming a<b
subroutine random_uniform(a,b,x)
   implicit none
   real,intent(in) :: a,b
   real,intent(out) :: x
   real :: u
   call random_stduniform(u)
   x = (b-a)*u + a
end subroutine random_uniform

Standard normal distribution

There ia a famous set of formula to generate the standard-normal deviates based on the standard uniform random number. It is called the Box-Muller transformation. The formula generates two normal-deviates (\(x_1\) and \(x_2\)) at the same time based on two uniform numbers \(u_1\) and \(u_2\).

\[ \begin{split} x_1 &= \sqrt{-2\ln(u_1)}\cos(2\pi u_2) \\ x_1 &= \sqrt{-2\ln(u_1)}\sin(2\pi u_2) \end{split} \]

We need only one of them, and we use the first formula.

subroutine random_stdnormal(x)
   implicit none
   real,intent(out) :: x
   real,parameter :: pi=3.14159265
   real :: u1,u2
   call random_stduniform(u1)
   call random_stduniform(u2)
   x = sqrt(-2*log(u1))*cos(2*pi*u2)
end subroutine random_stdnormal

Gamma distribution

The gamma distribution has alternative forms by different parameterization. Here, I take the following form \(\mathrm{Gamma}(\alpha,\beta)\) with shape parameter \(a>0\) and scale parameter \(\beta>0\). Note that an article in Wikipedia uses \(k=\alpha\) and \(\theta=\beta\).

\[ f(x)=\frac{1}{\Gamma(\alpha)\beta^\alpha}x^{\alpha-1}e^{-\frac{x}{\beta}} \]

If \(\alpha\) is integer, the Gamma deviation is easily available using the inverse CDF (cumulative distribution function) method. However, for the general case with any \(\alpha>0\), the derivation of algorithm needs much mathematical background. Here we first look at an algorithm for \(\alpha\ge 1\) (see Fishman, 1976 and Kennedy and Gentle, 1980).

  1. Generate \(y=-\ln(u_1)\) where \(u_1\sim U(0,1)\).
  2. Generate \(u_2\sim U(0,1)\). If \(u_2\le \left[y/\exp(y-1)\right]^{\alpha-1}\), return \(x=\alpha y\). Otherwise, go to 1.

The corresponding code is as follows.

subroutine random_stdgamma_alpha_ge_1(alpha,x)
   implicit none
   real,intent(in) :: alpha
   real,intent(out) :: x
   real :: y,t,u1,u2
   do
      call random_stduniform(u1)
      y = -log(u1)
      t = (y/exp(y-1))**(alpha-1)
      call random_stduniform(u2)
      if(u2 <= t) then
         x = alpha*y
         exit
      end if
   end do
end subroutine random_stdgamma_alpha_ge_1

If \(0<\alpha<1\), we can still use the above algorithm by shifting the parameter as \(\alpha+1\).

  1. Generate \(g\sim \mathrm{Gamma}(\alpha+1,1)\) using the above algorithm.
  2. Generate \(u\sim U(0,1)\), and return \(x=gu^{1/\alpha}\).

The complete subroutine for \(\alpha>0\) is available by combining the above algorithms.

subroutine random_stdgamma(alpha,x)
   implicit none
   real,intent(in) :: alpha
   real,intent(out) :: x
   real :: g,u
   if(alpha<=0) then
      stop "alpha<=0"
   else if(alpha<1) then
      call random_stdgamma_alpha_ge_1(alpha+1.0,g)
      call random_stduniform(u)
      x = g*u**(1.0/alpha)
   else
      call call random_stdgamma_alpha_ge_1(alpha,x)
   end if
end subroutine random_stdgamma

For the general gamma deviate with an arbitrary \(\beta>0\), you multiply \(\beta\) to the random number generated with the above subroutine.

subroutine random_gamma(alpha,beta,x)
   implicit none
   real,intent(in) :: alpha,beta
   real,intent(out) :: x
   call random_stdgamma(alpha,x)
   x = x*beta
end subroutine random_gamma

Back to index.html.