Box Muller Fortran

[JohnNichols]

The Statistical Techniques Research Group was the home for Muller of Box-Muller fame.  Box was the director. 

Muller published a research paper in 1958 for the Army (attached) about the method.   I have been using the MKL version of the Box-Muller method for some Monte Carlo Analysis.  

If you read the Muller paper and it is quite interesting, he talks about the issues with the method, in some detail.  

I found a simple implementation of the method from Japan.   Muller notes the real problems using this method for the extreme values, but for Monte Carlo the extremes are important,  we are looking at engineering failures so an error of 4 instead of 1 in a million is of interest in the work.  

Muller makes use of polynomials to approximate X = X(U). I have never used the Chebyshev polynomials, has anyone used these polynomials?  

I am not going to stop using the MKL procedures, this is just plain interest in understanding a dense paper. 

!****************************************************
!  正規乱数生成（by ボックス＝ミューラー法）
!  Normal random number generation (by box = Mueller method)
!  0 〜 20 の正規乱数をボックス＝ミューラー法により
!: Normal random numbers from 0 to 20 by the Box-Muller method
!  計算し、ヒストグラムを出力する。
!  Compute and output a histogram.
!
!  date          name            version
!  2018.12.03    mk-mode.com     1.00 新規作成
!
! Copyright(C) 2018 mk-mode.com All Rights Reserved.
!****************************************************
!
module const
  implicit none
  ! SP: 単精度(4), DP: 倍精度(8)
  ! SP: single precision (4), DP: double precision (8)
  integer,     parameter :: SP = kind(1.0)
  integer(SP), parameter :: DP = selected_real_kind(2 * precision(1.0_SP))
  integer(SP), parameter :: M  = 10                     ! 平均 ! average
  real(DP),    parameter :: S  = 2.5_DP                 ! 標準偏差 ! standard deviation
  integer(SP), parameter :: N  = 10000                  ! 発生させる乱数の個数 ! Number of random numbers to generate
  real(DP),    parameter :: PI = 4.0_DP * atan(1.0_DP)  ! 円周率 ! Pi 
  real(DP),    parameter :: SC = N / 100.0_DP           ! ヒストグラム用スケール ! Histogram Scale
end module const

module box_muller
  use const, only : SP, DP, S, M, PI
  implicit none
  private
  public :: rnd_seed, rnd

    contains
  ! 乱数の種の設定
    ! Set random number seed
  subroutine rnd_seed
    implicit none
    integer(SP) :: seed_size, clock
    integer(SP), allocatable :: seed(:)

    call system_clock(clock)
    call random_seed(size=seed_size)
    allocate(seed(seed_size))
    seed = clock
    call random_seed(put=seed)
    deallocate(seed)
  end subroutine rnd_seed

  ! 正規乱数生成
  ! normal random number generator
  ! :param(out) integer(4) r(2)
  subroutine rnd(r)
    integer(SP), intent(out) :: r(2)
    real(DP) :: r_u(2)  ! [0, 1] の一様乱数 2 個 ! 2 uniform random numbers in [0, 1]

    call random_number(r_u(1))
    call random_number(r_u(2))
    r(1) = int(S * sqrt(-2 * log(r_u(1))) * cos(2 * PI * r_u(2)) + M)
    r(2) = int(S * sqrt(-2 * log(r_u(1))) * sin(2 * PI * r_u(2)) + M)
  end subroutine rnd
end module box_muller

program rndnum_bm
  use const, only : SP, N, M, SC
  use box_muller
  implicit none
  integer(SP) :: i, j, r(2), hist(0:M * 2)

  ! 乱数の種の設定
  ! (正規乱数生成に使用する一様乱数の種)
  ! Set random number seed
  ! (uniform random number seed for normal random number generation)
  call rnd_seed

  ! 正規乱数の生成
  ! generate normal random numbers
  hist(:) = 0
  do i = 1 , N
    call rnd(r)
    hist(r) = hist(r) + 1
  end do

  ! 結果出力
  ! result output
  do i = 0, M * 2
    write (*, '(I3, ":", I4, " | ", A)') &
      & i, hist(i), repeat("*", int(hist(i) / SC))
  end do
end program rndnum_bm

[JohnNichols]

    !-----------------------------------------------------------------------
    !
    !   Returns the probability from [0 to 1] of Gaussian from beta value.
    !
    !-----------------------------------------------------------------------
    Real Function ProbGauss(prob)
    implicit none
    real, INTENT(IN) :: prob
    
    probGauss = 10.0**(-(0.3602 + 0.2052*prob + 0.2069*prob*prob))
    end Function ProbGauss

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

@Arjen_Markus  beta to probability value.  

[JohnNichols]

Intel  Corporation (2020). Intel® Fortran Compiler 19.0 User and Reference Guide Santa Clara, California, Intel  Corporation.

if I am referencing oneapi, what number do I give the manuals for a paper reference?

[Steve_Lionel]

See the title of this page.

Intel® Fortran Compiler Classic 2021.1 and Intel® Fortran Compiler...

[JohnNichols]

Thank you.  

I solved one of the standard steel problems from Gorenc and Tinyou  with Monte Carlo today. The probability of failure was 1e-17.  

About as long as humans have been on the planet for an event say one a day or a minute or thereabouts.  

Same chances I will meet Julia Roberts or Steve in person.