Sobol sequence by 바죠

Sobol sequence
https://en.wikipedia.org/wiki/Sobol_sequence
https://en.wikipedia.org/wiki/Low-discrepancy_sequence
https://en.wikipedia.org/wiki/Ilya_M._Sobol


랜덤 넘버, 막수는 컴퓨터 프로그램에서 매우 중요한 것이다. 특히, 유전 알고리듬, 몬테칼로 (Monte Carlo) 알고리듬에서와 마찬가지로 매우 중요한 의미를 가지고 있다.
따라서, 많은 경우, 랜덤 넘버 생성기는 컴퓨터 언어자체에서 내장함수 형식으로 제공된다.  
랜덤 넘버는 통상 함수로 만들어져 있으며, 초기화가 가능하다. 함수를 부를 때마다 새로운 값을 얻을 수 있다.
0.에서 1.0 사이의 값이 생성되게 된다. 몬테칼로 알고리듬에서 필요로하는 확률이라고 볼 수 있다.
특정한 행위가 받아들여질 확률이 바로 랜덤 넘버와 비교되는 것이다.

   x  ---------> x',
V(x) ----> V(x')

===>

exp[-{V(x')-V(x)}/(kBT)]  > tmp


포트란 90 언어의 경우,
call random_seed()라고 실행하면 랜덤 넘버 생성기를 초기화할 수 있다.
call random_number(tmp)처럼 부함수를 부르면 부를 때마다 새로운 랜덤 넘버를 얻어낼 수 있다.

일차원 또는 이차원 문제의 경우 매우 많은 시도로 부터 균일한 분포를 얻을 수 있다.
다차원으로 가게 되면 우리가 생각하고 있는 균일한 그러한 분포를 가질 수 있을까?  훨씬 더 많은 시도가 있어야하는 것은 자명하다. 유한한 갯수의 시도로 균일한 분포를 이야기하는 것 자체가 매우 까다롭고 신중하게 취급해야만 하는 것이다. 특히, 다차원의 문제풀이에서는 더욱더 신중해야만 할 것이다.  
수렴정도를 체크하기가 쉽지 않기 때문이다. 아래의 그림 두 개를 보자. 뚜렷한 차이를 볼 수 있다. 이차원 보기임에도 불구하고 확연히 다른 랜덤 넘버들의 분포를 볼 수 있다.  







http://people.sc.fsu.edu/~jburkardt/f_src/sobol_dataset/sobol_dataset.html
http://people.sc.fsu.edu/~jburkardt/f_src/sobol/sobol.f90

다운로드
받은 루틴은 체크를 하고 사용해야 한다. 한 가지 문제점을 발견했다.
integer type 정의가 잘못된 곳이 있었다.
아래의
컴파일러 옵션으로 쉽게 찾을 수 있었다.

ifort  -CB -check all  -warn interface -assume  realloc_lhs

아래에 표시된 적분을 한 번 실행해
보자. 정답이 알려진 2차원 정적분이다.
컴퓨터 알고리듬을 사용하여 실행하는 것을 목표로 한다.
두 가지 방식의 랜덤 넘버들을
이용하고 결과를 비교한다.
차이점을 확인할 수 있다.


많은 응용을 기대할 수 있다.


http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.473.490&rep=rep1&type=pdf


!234567890
       program  main
       implicit none
       integer mdim,n,iskip
       integer j,i,k,is,kase,ksgn,k1,k2
       real*8 tmp,vec(3)
       real ( kind = 8 ), allocatable, dimension ( :, : ) :: r
       call random_seed()
       call random_number(tmp)
       mdim=3
       mdim=2
       n=120000
       iskip=n*2**4+tmp*10000
       allocate(r(1:mdim,1:n))
       call i8_sobol_generate(mdim,n,iskip,r)
       if(mdim == 2)then
       r=r-0.5d0
                    endif
       if(mdim == 3)then
       r=r-0.5d0
                    endif
       tmp=0.d0
       do j=1,n
       tmp=tmp+4.d0*r(1,j)**2+3.d0*r(2,j)
       enddo
       write(6,*) tmp/dble(n)
       do j=1,n
       call random_number(tmp)
       r(1,j)=tmp-0.5d0
       call random_number(tmp)
       r(2,j)=tmp-0.5d0
       enddo
       tmp=0.d0
       do j=1,n
       tmp=tmp+4.d0*r(1,j)**2+3.d0*r(2,j)
       enddo
       write(6,*) tmp/dble(n)
       if(allocated(r)) deallocate(r)
       stop
       end
  0.333351784496133    
  0.338188940039489
    



http://incredible.egloos.com/4832215




핑백

덧글

댓글 입력 영역

최근 포토로그



MathJax