- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
Link Copied
- « Previous
-
- 1
- 2
- Next »
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
!-----------------------------------------------------------------------
!
! 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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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?
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
![](/skins/images/7B0AB6865064EAE32B30EC9A4E94B48A/responsive_peak/images/icon_anonymous_message.png)
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page
- « Previous
-
- 1
- 2
- Next »