random number 활용하기 by 바죠

컴퓨터 시뮬레이션 분야에서는 랜덤 넘버(난수, 막수) 생성기를 활용한 여러 가지 작업들을 수행한다.  랜덤 넘버들을 어떻게 만들어서 사용하는지는 꽤 까다로운 과학적 주제이다. 많은 논의가 있을 수 있다.  원론적으로 그렇다는 것이다. 병렬 계산에서는 또 다른 이쓔를 만들어 낼 수 있다. 하지만, 대체로 유사 랜덤넘버를 사용한다. 많은 경우 각종 컴퓨터 언어가 지원하는 함수를 사용하면 별 문제가 없다. 대부분의 경우가 그렇다는 뜻이다.  

cf.
http://incredible.egloos.com/4832215
http://incredible.egloos.com/4837144

http://incredible.egloos.com/7378157



대부분의 경우, 다른 사람들이 만들어둔 랜덤 넘버들을 테스트하고 잘 사용할 수 있는 수준에 이르면 된다.  사용하기전에 몇 가지 테스트를 해 볼 수 있을 것이다. 통계물리학 분야에서는 소위 주사위만 잘 던져도 물리적인 해석이 가능하다고 말할 수 있다.


통계역학에서처럼 통계적인 확률에 의해서 측정 가능한 물리량이 결정된다. 따라서, 컴퓨터 시늉에서도 동일한 방법으로 이러한 물리량 계산이 가능하다. 가장 잘 알려진 것이 Monte Carlo 시늉이다. 예를 들어, x-->x' 으로 변환되는 것을 가정할 경우, 이러한 변환에 대한 확률이 알려진다면 P(x'), P(x)가 각각 알려진다면, 우리는 랜덤 함수를 이용하여, x-->x' 변환을 매우 많이 수행할 수 있다. 변환이 받아들여지고 안 받아들여지고를 확률적으로 판단한다는 뜻이다. 이러한 시늉이 굉장히 많은 수에 대해서 적용된다면 충분히 수렴된 물리량도 계산을 할 수 있다.

그림 출처: google.co.kr 이미지



진짜 랜덤 넘버는 꿰 까다로운 조건들을 만족해야 한다.  실제로 사용되는 랜덤 넘버 생성기를 통상 슈도 랜덤넘버라고 한다. 즉, 상당한 수준으로 랜덤 넘버에 가깝다는 뜻이다.  진짜는 아니지만, 상당히 높은 수준이라는 뜻이다.

항상 다시 만들어 낼 수 있는 랜덤 넘버들을 만들기 위해서 씨드를 사용한다. 이 씨드를 알고 있으면 항상 동일한 랜덤 넘버들을 만들어 낼 수 있다.

!234567890
       implicit none
       integer :: values(1:8), k
       integer, dimension(:), allocatable :: iseed
       real(8) :: r

       call date_and_time(values=values)

       call random_seed(size=k)
       allocate(iseed(1:k))
       iseed(:) = values(8)
       iseed(:)=0
       iseed(1)=124
       call random_seed(put=iseed)
       call random_number(r)
       print *,r

       stop
       end


random integer generation
random sequence generation
Gaussian generation
http://en.wikipedia.org/wiki/Random_number_generator
http://www.random.org/

http://jblevins.org/mirror/amiller/random.f90
http://www.fortranlib.com/freesoft.htm

통계적으로 일을 처리할 때, 난수를 발생시켜서 사용하는 경우가 있다. 예를 들어, 주사위를 던지면서 나오는 숫자를 활용하는 시뮬레이션을 생각할 수 있다.

1,2,3,4,5,6 주사위의 경우는 아래와 같다.

i=(6.d0)*ranmar()+1.d0

이 때, 즐겨 사용하는 함수가 있다.  이 함수를 만드는 방법보다 사용하는 것이 더 중요하다. 만드는 방법은 굉장히 많이 있다. 만들어진 함수는  균일한 분포로 난수 숫자를 만들어낼 수 있어야 한다. 통상 0.0에서 1.0 사이의 숫자를 균등하게 불규칙하게 만들어 낼 수 있으면 된다. ranmar()가 0에서 1.0이하의 임의의 숫자라고 가정하자. [0,1) 이렇게 표시한다. 1.0은 생성하지 않는다는 뜻이다.

얼마나 좋은 랜덤 넘버 생성기인지는 아래와 같은 수식으로 테스트해 볼 수 있다.



아래와 같이 실행하면 

       i=dble(n)*ranmar()

i가 가질 수 있는 숫자는 아래와 같다.
       0,1,2,...n-1

다시 말해서, 아래와 같이 수행하면
        i=dble(n)*ranmar()+1
i가 가질 수 있는 숫자는 아래와 같다.
        1,2,3,....,n



만약, k0=10 이라고 하면 평균이 10이 면서, 표준 편차가 10이 되는 exponential 분포를 가지는 값 (kk) 하나를 만들 수도 있다.
        kk=-float(k0)*log(ranmar())


a exp(-a x)
http://mathworld.wolfram.com/ExponentialDistribution.html



g=sqrt(-2*log(ranmar()))*cos(2*pi*ranmar())




       implicit none
       integer nseed
       real*8 tmp,tmq,aa
       real*8 histo(0:100000000)
       real*8 avg,sig,tmr
       integer i

       nseed=222
       call random_seed(nseed)
       call random_number(tmp)

       histo=0.d0

       avg=0.d0
       sig=0.d0
       tmr=0.d0

       aa=1.d0
       aa=100.d0
       do i=1,1000000
       call random_number(tmp)
       tmq=-aa*log(tmp)
       avg=avg+tmq
       sig=sig+tmq**2
       tmr=tmr+1.d0
       histo(int(tmq*10))=histo(int(tmq*10))+1.d0
       enddo

       avg=avg/tmr
       sig=sig/tmr
       sig=sqrt(sig-avg**2)
       print*, avg,sig

       open(9,file='fort.9',form='formatted')
       do i=0,10000
       write(9,*) i, histo(i)
       enddo
       close(9)

       stop
       end


       implicit none
       integer nn
       real*8 pi
       real*8 x(1000000),xl,xl0,stddev
       real*8 avg,sig,tmp,tmq
       integer i,j
       integer iseed1,iseed2
       real ranmar

       iseed1=1101
       iseed2=1121
       call rmarin(iseed1,iseed2)

       pi=4.d0*atan(1.d0)
       xl0=0.d0
       stddev=1.d0


       nn=100000
       avg=0.d0
       do i=1,nn
       call gauss(stddev,xl0,xl)
       x(i)= xl
       avg=avg+x(i)
       enddo
       avg=avg/float(nn)
       print*, avg
       sig=0.d0
       do i=1,nn
       sig=sig+(x(i)-avg)**2
       enddo
       sig=sqrt(sig/float(nn))
       print*, sig

 

       stop
       end
!
!234567890
       subroutine gauss(sigma,xl0,xl)
!      Written by In-Ho Lee, KRISS, September 18 (2007).
       IMPLICIT NONE
       real*8 sigma,xl0,xl
       real*8 tmp,r,v1,v2
       real ranmar
       r=2.0d0
       do while (r >= 1.d0)
       v1=2.0*ranmar()-1.0
       v2=2.0*ranmar()-1.0
       r=v1**2+v2**2
       enddo
       xl=v1*sqrt(-2.d0*log(r)/r)
       xl=xl0+sigma*xl
       return
       end








--------------------------------------------------------------------------------------------------------------------

call random_seed()

초기화 하고

그 다음 부터는 아래와 같은 방식으로 부함수를 불러서 사용할 수 있다.

call random_number(tmp)

--------------------------------------------------------------------------------------------------------------------




아래에 있는 랜덤 넘버 생성기는 약간 독특하다.
우선 초기화를 위해서 아래의 서브루틴을 부른다. 이 때, 씨드 정수로 두 개가 필요하다.
초기화 하고 나면, real*8이 아닌 real형으로 랜덤 넘버를 생성시켜준다.

!
!234567890
      subroutine rmarin(ij,kl)
!  This subroutine and the next function generate random numbers. See
!  the comments for SA for more information. The only changes from the
!  orginal code is that (1) the test to make sure that RMARIN runs first
!  was taken out since SA assures that this is done (this test didn't
!  compile under IBM's VS Fortran) and (2) typing ivec as integer was
!  taken out since ivec isn't used. With these exceptions, all following
!  lines are original.

! This is the initialization routine for the random number generator
!     RANMAR()
! NOTE: The seed variables can have values between:    0 <= IJ <= 31328
!                                                      0 <= KL <= 30081
      real u(97), c, cd, cm
      integer i97, j97
      common /raset1/ u, c, cd, cm, i97, j97
      if( ij .lt. 0  .or.  ij .gt. 31328  .or.  &
          kl .lt. 0  .or.  kl .gt. 30081 ) then
          print '(a)', ' The first random number seed must have a value  &
    & between 0 and 31328'
          print '(a)',' The second seed must have a value between 0 and  &
    & 30081'
            stop
      endif
      i = mod(ij/177, 177) + 2
      j = mod(ij    , 177) + 2
      k = mod(kl/169, 178) + 1
      l = mod(kl,     169)
      do 2 ii = 1, 97
         s = 0.0
         t = 0.5
         do 3 jj = 1, 24
            m = mod(mod(i*j, 179)*k, 179)
            i = j
            j = k
            k = m
            l = mod(53*l+1, 169)
            if (mod(l*m, 64) .ge. 32) then
               s = s + t
            endif
            t = 0.5 * t
  3      continue
         u(ii) = s
  2   continue
      c = 362436.0 / 16777216.0
      cd = 7654321.0 / 16777216.0
      cm = 16777213.0 /16777216.0
      i97 = 97
      j97 = 33
      return
      end
!
!234567890
      function ranmar()
      real u(97), c, cd, cm
      integer i97, j97
      common /raset1/ u, c, cd, cm, i97, j97
         uni = u(i97) - u(j97)
         if( uni .lt. 0.0 ) uni = uni + 1.0
         u(i97) = uni
         i97 = i97 - 1
         if(i97 .eq. 0) i97 = 97
         j97 = j97 - 1
         if(j97 .eq. 0) j97 = 97
         c = c - cd
         if( c .lt. 0.0 ) c = c + cm
         uni = uni - c
         if( uni .lt. 0.0 ) uni = uni + 1.0
         ranmar = uni
      return
      end


http://incredible.egloos.com/4779646
http://incredible.egloos.com/4768344


       implicit none
       integer nseed
       real*8 avg,sig
       real*8 tmp,tmq
       real*8 histo(100)
       integer i,j,k
       real*8 xmu(3)
       integer nsize
       integer, allocatable :: iseed(:)

       nseed=111
       call random_seed(size=nsize)
       if(.not. allocated(iseed)) allocate(iseed(nsize))
       do i=1,nsize
       iseed(i)=i+nseed
       enddo
       call random_seed(put=iseed)
       if(.not. allocated(iseed)) deallocate(iseed)


       do i=1,10
       call random_number(tmp)
       print*, tmp
       enddo



포트란 90에서는 랜덤 넘버 생성기를 지원한다. 아래와 같이 사용하면 된다. random_seed, random_number라는 두 루틴만 사용할 수 있으면 된다.  이를 위해서, 정수, 정수 배열이 필요하다.   씨드를 공급하고, 다시 사용하게될  씨드를 확인할 수 있으면 된다.
random_seed()에서는
size, put, get이 각각 가능하다.
size를 알아내고, seed를 선택해서 공급하고, 다시 사용할 seed를 알아 낼 수 있다.

포트란 90에서 지원하는 random 넘버를 사용하는 방법을 알아 둘 필요가 있다.


implicit none

integer :: i,size
integer, allocatable :: seed(:)
real :: r

call random_seed(size)
allocate (seed(size))
write(*,*)'give ',size,' random seeds '
read(*,*)seed
write(*,*)

call random_seed(put=seed)
do i=1,10
call random_number(r)
write(*,*)r
end do

write(*,*)
call random_seed(get=seed)
write(*,*)seed

end



!--------------------------------------------------------!
!Random number generator based on algorithm 'mzran' by !
!Marsiglia and Zaman, Computers in Physics, 8, 107 (1994)!
!--------------------------------------------------------!
!rand() generates a random number in the range [0,1) !
!initrand(w) reads 4 integer seeds from 'seed.in', !
!writes new random seeds to the same file if w/=0. !
!--------------------------------------------------------!

!-----------------------!
real(8) function rand()
!-----------------------!
implicit none

integer :: iir,jjr,kkr,nnr,mzran
common/brand/iir,jjr,kkr,nnr,mzran

mzran=iir-kkr
if (mzran < 0) mzran=mzran+2147483579
iir=jjr; jjr=kkr; kkr=mzran
nnr=69069*nnr+1013904243
mzran=mzran+nnr
rand=0.5d0+mzran*0.23283064d-9

end function rand
!-----------------!

!----------------------!
subroutine initrand(w)
!----------------------!
implicit none

integer :: iir,jjr,kkr,nnr,mzran
common/brand/iir,jjr,kkr,nnr,mzran

integer :: w
real(8) :: rand

open(1,file='seed.in',status='old')
read(1,*)iir
read(1,*)jjr
read(1,*)kkr
read(1,*)nnr
close(1)
iir=abs(iir)+1; jjr=abs(jjr)+1; kkr=abs(kkr)+1
if (w /= 0) then
open(1,file='seed.in',status='replace')
write(1,*)abs(int((rand()-0.5d0)/0.23283064d-9))
write(1,*)abs(int((rand()-0.5d0)/0.23283064d-9))
write(1,*)abs(int((rand()-0.5d0)/0.23283064d-9))
write(1,*)abs(int((rand()-0.5d0)/0.23283064d-9))
close(1)
end if

end subroutine initrand
!-----------------------!


http://arxiv.org/pdf/1005.4117v1.pdf

E(N,n) = 1/N \sum_{i=1}^{N} x_{i}x_{i+n} -E^2
E=1/N \sum_{i=1}{N} x_i

M(N,k)=abs(1/N\sum x_i^k-1/(k+1))


!234567890
module ran_mod
! ran1 returns a uniform random number between 0-1
! the seed is optional and used to reset the generator
contains
function ran1(my_seed)
implicit none
real(8) ran1,r
integer, optional ,intent(in) :: my_seed ! optional argument not changed in the routine
integer,allocatable :: seed(:)
integer the_size,j
if(present(my_seed))then ! use the seed if present
call random_seed(size=the_size) ! how big is the intrisic seed?
allocate(seed(the_size)) ! allocate space for seed
do j=1,the_size ! create the seed
seed(j)=abs(my_seed)+(j-1) ! abs is generic
enddo
call random_seed(put=seed) ! assign the seed
deallocate(seed) ! deallocate space
endif
call random_number(r)
ran1=r
end function ran1
end module

랜덤 넘버 씨드를 넣는 것도 귀찮은 경우:
그 때 그 때 달라지는 씨드의 구현은 아래와 같이 하면 된다.  게임과 같은 소프트웨어에서 이러한 작업이 필요할 것이다.

      integer Count(1)   !Current value of system clock

!     initialize random number generator
      call System_Clock(Count(1))
      call Random_Seed(Put=Count)

      call Random_Number(x)  !   returns random x >= 0 and <1  


       character*8 date
       character*10 time
       character*5 zone
       integer ivalues(8)
       integer i,j
       integer isize,kount(1)
       integer, allocatable :: iseed(:)

       call date_and_time(date,time,zone,ivalues)
!      write(6,*) date,time,zone
       write(6,*) ivalues
       call random_seed(isize)
       allocate(iseed(isize))
       iseed=0
       call system_clock(kount(1),i,j)
       kount(1)=ivalues(8)
!      write(6,*) kount(1)

       i=min(isize,8)
       iseed(1:i)=iseed(1:i)+ivalues(1:i)
       do i=1,isize
       iseed(i)=iseed(i)+i
       enddo
       call random_seed(put=iseed)
       deallocate(iseed)







http://lib.stat.cmu.edu/apstat/183
         real function random()
c
c     Algorithm AS 183 Appl. Statist. (1982) vol.31, no.2
c
c     Returns a pseudo-random number rectangularly distributed
c     between 0 and 1.   The cycle length is 6.95E+12 (See page 123
c     of Applied Statistics (1984) vol.33), not as claimed in the
c     original article.
c
c     IX, IY and IZ should be set to integer values between 1 and
c     30000 before the first entry.
c
c     Integer arithmetic up to 30323 is required.
c
      integer ix, iy, iz
      common /randc/ ix, iy, iz
c
      ix = 171 * mod(ix, 177) - 2 * (ix / 177)
      iy = 172 * mod(iy, 176) - 35 * (iy / 176)
      iz = 170 * mod(iz, 178) - 63 * (iz / 178)
c
      if (ix .lt. 0) ix = ix + 30269
      if (iy .lt. 0) iy = iy + 30307
      if (iz .lt. 0) iz = iz + 30323
c
c     If integer arithmetic up to 5212632 is available, the preceding
c     6 statements may be replaced by:
c
c     ix = mod(171 * ix, 30269)
c     iy = mod(172 * iy, 30307)
c     iz = mod(170 * iz, 30323)
c
      random = mod(float(ix) / 30269. + float(iy) / 30307. +
     +                        float(iz) / 30323., 1.0)
      return
      end
c
c
c
c
 real function uniform()
c
c Generate uniformly distributed random numbers using the 32-bit
c generator from figure 3 of:
c L'Ecuyer, P. Efficient and portable combined random number
c generators, C.A.C.M., vol. 31, 742-749 & 774-?, June 1988.
c
c The cycle length is claimed to be 2.30584E+18
c
c Seeds can be set by calling the routine set_uniform
c
c It is assumed that the Fortran compiler supports long variable
c names, and integer*4.
c
 integer*4 z, k, s1, s2
 common /unif_seeds/ s1, s2
 save /unif_seeds/
c
 k = s1 / 53668
 s1 = 40014 * (s1 - k * 53668) - k * 12211
 if (s1 .lt. 0) s1 = s1 + 2147483563
c
 k = s2 / 52774
 s2 = 40692 * (s2 - k * 52774) - k * 3791
 if (s2 .lt. 0) s2 = s2 + 2147483399
c
 z = s1 - s2
 if (z .lt. 1) z = z + 2147483562
c
 uniform = z / 2147483563.
 return
 end


 subroutine set_uniform(seed1, seed2)
c
c Set seeds for the uniform random number generator.
c
 integer*4 s1, s2, seed1, seed2
 common /unif_seeds/ s1, s2
 save /unif_seeds/

 s1 = seed1
 s2 = seed2
 return
 end


       character*8 date
       character*10 time
       character*5 zone
       integer ivalues(8)
       integer i,j
       integer isize,kount(1)
       integer, allocatable :: iseed(:)

       call date_and_time(date,time,zone,ivalues)
!      write(6,*) date,time,zone
       write(6,*) ivalues
       call random_seed(isize)
       allocate(iseed(isize))
       iseed=0
       call system_clock(kount(1),i,j)
       kount(1)=ivalues(8)
!      write(6,*) kount(1)

       i=min(isize,8)
       iseed(1:i)=iseed(1:i)+ivalues(1:i)
       do i=1,isize
       iseed(i)=iseed(i)+i
       enddo
       call random_seed(put=iseed)
       deallocate(iseed)


--------------------------------------------------------------------------------------------------------------------
http://incredible.egloos.com/4832215


실행 시간 기준으로 랜덤넘버 자동 씨드 만들기, 포트란 90

!234567890
      program test
      implicit none
      real*8 tmp
      integer i

      call init_seed()

      do i=1,10
      call random_number(tmp)
      write(6,*) tmp
      enddo

      end program test

      subroutine init_seed()
      implicit none
      integer n, ival(8), v(3), i
      integer, allocatable :: seed(:)
      call date_and_time(values=ival)
      v(1) = ival(8) + 2048*ival(7)
      v(2) = ival(6) + 64*ival(5) ! value(4) isn't really 'random'
      v(3) = ival(3) + 32*ival(2) + 32*8*ival(1)
      call random_seed(size=n)
      allocate(seed(n))
      call random_seed()          ! Give the seed an implementation-dependent kick
      call random_seed(get=seed)
      do i=1, n
      seed(i) = seed(i) + v(mod(i-1, 3) + 1)
      enddo
      call random_seed(put=seed)
      deallocate(seed)
      end subroutine
                       
--------------------------------------------------------------------------------------------------------------------
program test_random_number  
REAL(8) :: r  
CALL random_seed()  
CALL RANDOM_NUMBER(r)  
print *, r
end program

Reference https://www.physicsforums.com/threads/fortran-90-how-do-i-use-random_number-and-random_seed.551106/

--------------------------------------------------------------------------------------------------------------------
subroutine init_pseudorandom

    integer, parameter                 :: maxr = HUGE(0) - 1
    integer(kind=8)                    :: count,countrate,countmax,pid
    integer                            :: i,j,iters,size
    integer, allocatable, dimension(:) :: seed
    real(kind=dp)                      :: rr

    pid=getpid()

    ! ** Seed the random number generator
   
    call random_seed(size=size)
   
    allocate(seed(size))
   
    call random_seed(get=seed)

    do j=1,size
       call system_clock(count,countrate,countmax)
       iters = int(mod(count,pid),4)
       do i=1,iters+1
          call random_number(rr)
       end do      
       seed(j) = int(2*(rr-0.5_dp)*maxr)
    end do

    call random_seed(put=seed)

    deallocate(seed)
   
  end subroutine init_pseudorandom




핑백

덧글

  • 바죠 2019/10/09 12:01 # 답글

    !234567890
    implicit none
    real*8 tmp
    integer i,j,k

    call random_seed()
    do i=1,10
    call random_number(tmp)
    j=tmp*10+1
    write(6,*) j
    enddo
    stop
    end
  • 바죠 2019/10/09 12:17 # 답글

    !234567890
    subroutine getij(ii,jj,npop,isel)
    implicit none
    integer ii,jj,npop,isel
    real*8 tmp
    integer i,j,k,m

    if(isel == 1)then
    call random_number(tmp)
    ii=tmp*npop+1
    endif
    if(isel == 2)then
    call random_number(tmp)
    k=tmp*npop+1
    call random_number(tmp)
    m=tmp*npop+1
    i=min(k,m)
    22 continue
    call random_number(tmp)
    k=tmp*npop+1
    call random_number(tmp)
    m=tmp*npop+1
    j=min(k,m)
    if(i == j) goto 22
    ii=i ; jj=j
    endif
    if(isel == 3)then
    call random_number(tmp)
    i=npop/2 ; i=-i*log(tmp+1.d-33)
    call random_number(tmp)
    if(i <= 0 .or. i > npop) i=npop*tmp+1
    33 continue
    call random_number(tmp)
    j=npop/2 ; j=-j*log(tmp+1.d-33)
    call random_number(tmp)
    if(j <= 0 .or. j > npop) j=npop*tmp+1
    if(i == j) goto 33
    ii=i ; jj=j
    endif
    end


    implicit none
    integer ii,jj,npop,isel
    npop=50

    call random_seed()

    isel=1
    call getij(ii,jj,npop,isel)
    write(6,*) ii

    isel=2
    call getij(ii,jj,npop,isel)
    write(6,*) ii,jj

    isel=3
    call getij(ii,jj,npop,isel)
    write(6,*) ii,jj
    end
  • 바죠 2020/03/17 12:37 # 답글

    import random
    import time
    random.seed(time.time())
댓글 입력 영역

최근 포토로그



MathJax