# Random numbers

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

- 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.
- Confirm that you get the same sequence of random numbers if you put the same seed.
- 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).

- Generate \(y=-\ln(u_1)\) where \(u_1\sim U(0,1)\).
- 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\).

- Generate \(g\sim \mathrm{Gamma}(\alpha+1,1)\) using the above algorithm.
- 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.