Générer des nombres aléatoires depuis une loi normale en fortran 90

Loi normale centrée réduite

Générer des nombres aléatoires depuis une loi normale centrée réduite (ou loi normale standard) en fortran 90 (source):

program test_random_number

implicit none

integer :: i

real(kind=8), external :: rn_std_normal_dist

!----------------------------------------------------------------------------------------!
! test standard normal distribution

open(1,file='rn_standard_normal_distribution.txt')

do i = 1, 100000

    write(1,*)  rn_std_normal_dist()

end do

close(1)

end program test_random_number

!----------------------------------------------------------------------------------------!
! source: http://sepwww.stanford.edu/sep/prof/geelib/random.f90; Author: Alan Miller

real(kind=8) function rn_std_normal_dist()

real :: half = 0.5

real :: s = 0.449871, t = -0.386595, a = 0.19600, b = 0.25472, &
        r1 = 0.27597, r2 = 0.27846, u, v, x, y, q

do

  call random_number(u)
  call random_number(v)
  v = 1.7156 * (v - half)

  x = u - s
  y = ABS(v) - t
  q = x**2 + y*(a*y - b*x)

  if (q < r1) exit

  if (q > r2) cycle

  if (v**2 < -4.0*log(u)*u**2) exit

end do

rn_std_normal_dist = v/u

end function rn_std_normal_dist

!----------------------------------------------------------------------------------------!

Loi normale

Si on sait comment générer des nombres aléatoires depuis une loi normale centrée réduite, on peut alors facilement générer des nombres aléatoires depuis une loi normale quelconque avec la formule $$X = Z * \sigma + \mu$$
ou Z est un nombres aléatoire généré depuis une loi normale centrée réduite, $\sigma$ la standard deviation et $\mu$ la moyenne.

!----------------------------------------------------------------------------------------!
! test normal distribution (mean: 10; sigma: 2)

open(1,file='rn_normal_distribution.txt')

do i = 1, 100000

    write(1,*)  rn_std_normal_dist() * 2.0 + 10

end do

close(1)
!----------------------------------------------------------------------------------------!

Tracer les résultats avec python 2.7

import matplotlib.pyplot as plt
import numpy as np
import math

#----------------------------------------------------------------------------------------#
# plot normal_distribution

def normal_distribution(x,mu,sigma):
    num = math.exp( - ( x - mu )**2 / ( 2.0 * sigma **2 ) )
    den = math.sqrt( 2 * math.pi * sigma **2)
    return  num / den

x_array = np.linspace(-5.0, 15.0, 100)

#y_array = np.asarray( [normal_distribution(x,0,1) for x in x_array] )
#plt.plot(x_array,y_array)

y_array = np.asarray( [normal_distribution(x,10,2) for x in x_array] )
plt.plot(x_array,y_array)

#----------------------------------------------------------------------------------------#
# plot data

#data1 = np.loadtxt('rn_standard_normal_distribution.txt')
#plt.hist(data1, bins=50,normed=1,color='lightblue')

data2 = np.loadtxt('rn_normal_distribution.txt')
plt.hist(data2, bins=50,normed=1,color="lightcoral")

plt.grid()

#plt.savefig("rn_standard_normal_distribution.png", bbox_inches='tight')
#plt.title('Random numbers from a standard normal distribution')

plt.savefig("rn_normal_distribution.png", bbox_inches='tight')
plt.title('Random numbers from a normal distribution (10,2)')

plt.show()

Références