이 코드가 나에게 잘못된 대답을주는 이유는 무엇입니까?

Question

상자 뮬러 알고리즘을 사용하여 더 많은 난수를 계산하는 함수로 lcg를 사용하는 프로그램을 작성하려고합니다. lcg를 작동 시켰지만 상자 뮬러 알고리즘을 사용하는 함수가 잘못된 값을 제공합니다.

module rng 
    implicit none 

    integer, parameter :: dp = selected_real_kind(15,300) 
    real(kind=dp) :: A=100, B= 104001, M = 714025 

contains 

function lcg(seed) 

    integer, optional, intent(in) :: seed 
    real(kind=dp) :: x = 0, lcg 

    if(present(seed)) x = seed 
    x = mod(A * x + B, M) 
    lcg = x/714025 

end function 


function muller(seed) 
    integer, parameter :: dp = selected_real_kind(15,300) 
    integer, optional, intent(in) :: seed 
    real(kind = dp) :: y1, y2, mean = 0.49, sd = 0.5, muller1, muller2, 
muller, x1, x2, pi = 4.0*ATAN(1.0) 
    integer :: N = 0 

! I had to add the do while loop to ensure that this chunk of code would 
only execute once 

do while (N<1) 

    x1 = lcg() 
    x2 = lcg() 
    N = N + 1 
    y1 = sd * SQRT(-2.0*LOG(x1)) * COS(2*pi*(x2)) + mean 
    y2 = sd * SQRT(-2.0*LOG(x1)) * SIN(2*pi*(x2)) + mean 

    print *, y1, y2, x1, x2 ! Printing x1 and x2 to allow me to use a 
calculator to check program is working correctly 
end do 

end function 

end module 

program lcgtest 
    use rng 
    implicit none 
    integer :: N 

    real(kind=dp) :: lcgmean, ttl = 0, sumof, lcgvar, dev1, muller1, muller2, 
lcgerr, lcgdev 
    real, dimension(10000) :: array 


do N = 1, 10000 

    ttl = ttl + lcg() 
    dev1 = lcg() - lcgmean 
    sumof = sumof + dev1 

end do 
    muller1 = muller() 
    muller2 = muller() 
    lcgmean = ttl/10000 
    lcgvar = ((sumof)**2)/10000 
    lcgdev = SQRT((sumof)**2)/10000 
    lcgerr = lcgdev/100 
    print *, lcg(), "mean=", lcgmean, "variance=", lcgvar, lcgerr 


end program 

의 핵심 부분은 뮬러 기능 섹션 :

여기 내 코드입니다. 계산기로 얻은 값을 확인한 후에, y1과 y2에 대한 답이 다른 것을 볼 수 있습니다.

도움을 주시면 감사하겠습니다.

Answer 1

이 프로그램의 결과로 기대되는 것이 무엇인지 알 수 없습니다. 그러나 그것을 읽으면 쉽게 따라 할 수있는 논리를 얻을 수 있습니다. 나는 두 가지 사실적인 실수를 느낀다. 제 개선을 위해 아래 코드를 읽으십시오.

module rng 
    implicit none 
    integer, parameter :: dp = selected_real_kind(15,300) 
    real(kind=dp) :: A=100, B= 104001, M = 714025 

contains 

    function lcg(seed) 

     integer, optional, intent(in) :: seed 
     real(kind=dp) :: x = 0, lcg 

     if(present(seed)) x = seed 
     x = mod(A * x + B, M) 
     lcg = x/714025 

    end function lcg ! Note 'lcg' here @ the end 


    function muller(seed) 
     integer, parameter :: dp = selected_real_kind(15,300) 
     integer, optional, intent(in) :: seed 
     real(kind = dp) :: y1, y2, mean = 0.49, sd = 0.5, muller1, & 
          muller2, muller, x1, x2, pi = 4.0*ATAN(1.0) 
     integer :: N = 0 

     ! I had to add the do while loop to ensure that this chunk 
     ! of code would only execute once 

     do while (N<1) 
      x1 = lcg() 
      x2 = lcg() 
      N = N + 1 
      y1 = sd * SQRT(-2.0*LOG(x1)) * COS(2*pi*(x2)) + mean 
      y2 = sd * SQRT(-2.0*LOG(x1)) * SIN(2*pi*(x2)) + mean 

      ! Printing x1 and x2 to allow me to use a 
      ! calculator to check program is working correctly 
      print *, y1, y2, x1, x2 
     enddo 
    end function muller ! note the function name @ the end here 

end module rng ! Note 'rng' here added. 

program lcgtest 
    use rng 
    implicit none 
    integer :: N 

    real(kind=dp) :: lcgmean, ttl = 0, sumof, lcgvar, dev1, muller1, & 
        muller2, lcgerr, lcgdev 
    real, dimension(10000) :: array 


    ! In the original code the variables 'lcgmean' and 'dev1' were  
    ! *undefined* before they were used in the do-loop. This will cause the 
    ! compiler to set them some random garbage values, and it will 
    ! inevitably leads to unexpected result or error in most cases. 

    ! In, order to avoid this by setting them. 
    ! For example, lcgmean = 1.0 and dev1 = 0.1 
    ! We'll then have the following: 
    lcgmean = 1.0 
    dev1 = 0.1 
    do N = 1, 10000 
     ttl = ttl + lcg() 
     dev1 = lcg() - lcgmean 
     sumof = sumof + dev1 
    end do 

    muller1 = muller() 
    muller2 = muller() 
    lcgmean = ttl/10000 
    lcgvar = ((sumof)**2)/10000 
    lcgdev = SQRT((sumof)**2)/10000 
    lcgerr = lcgdev/100 
    print *, lcg(), "mean=", lcgmean, "variance=", lcgvar, lcgerr 

end program 

추가 제안

항상 블록을-의 코드에 가까운 의미 같이 (예를 들어) 나는 종종 매우 유용 찾을 :

real function func_name(arg1, arg2, ....) 
    implicit none 

    .... 
end function func_name 

subroutine sub_name(arg1, arg2, ...) 
    implicit none 

    ... 
end subroutine sub_name 

비 변수 seed은 함수 muller에서 사용되지 않습니다. 어쩌면 필요하지 않을 수도 있습니다. 나는 here을 읽은 것을 감안할 때

Answer 2

, 동시에 Y1 및 Y2을 계산할 수있는 서브 루틴 그래서 의해 기능 뮬러를 대체 할 더 나은 것 같다. 사실, "블록"*의 뮬러 "목적은 두 개의 이전에 생성 된 의사 난수 X1 및 X2하여 프로그램 구조에 따른.이어서

에에 기초하여 두 개의 다른 의사 난수를 생성하는 주요 프로그램은 대신 기능로 뮬러를 호출하면 해당 위치에 전화 뮬러를 작성하여 서브 루틴으로 호출해야합니다. 그러나, 기능 아니라 이상의 서브 루틴의 사용이 여전히 가능하지만 두 값 y1 및 y2을 반환하려면 v (1) = y1 인 벡터 v를 반환 할 수 있습니다. v (2) = y2.

원래의 프로그램은 다음과 같이 될 것이다 :

module rng 
    implicit none 
    integer, parameter :: dp = selected_real_kind(15,300) 
    real(kind=dp) :: A=100, B= 104001, M = 714025 

contains 

    function lcg(seed) 
     implicit none  
     integer, optional, intent(in) :: seed 
     real(kind=dp) :: x = 0, lcg 

     if(present(seed)) x = seed 
     x = mod(A * x + B, M) 
     lcg = x/714025 

    end function lcg ! Note 'lcg' here @ the end 

    !------------------------------------------------------------------- 
    ! Subroutine muller. 
    ! Here, we supply 4 arguments *y1, y2* and the optional 
    ! argaument *seed* which apparently is not required since it is not 
    ! used (but can be used in order to have a better pseudo number  
    ! generator. 
    !------------------------------------------------------------------- 
    subroutine muller(y1, y2, seed) 
     implicit none 
     real(kind=dp), intent(out)  :: y1, y2 
     integer, optional, intent(in) :: seed 

     ! Local variables 
     real(kind=dp)     :: x1, x2 
     real(kind=dp)     :: mean = 0.49, sd = 0.5 
     real(kind=dp)     :: pi = 4.0*ATAN(1.0) 
     integer      :: N = 0 

     ! The **do while** loop is not needed 
     ! do while (N<1) 
     x1 = lcg() 
     x2 = lcg() 
     N = N + 1 
     y1 = sd * SQRT(-2.0*LOG(x1)) * COS(2*pi*(x2)) + mean 
     y2 = sd * SQRT(-2.0*LOG(x1)) * SIN(2*pi*(x2)) + mean 

     ! display to the screen the values of x1, x2, y1, y2 
     print *, y1, y2, x1, x2 
     ! enddo 
    end subroutine muller 
end module rng 

program lcgtest 
    use rng 
    implicit none 
    integer :: N 
    real(kind=dp) :: lcgmean, ttl = 0, sumof, lcgvar, dev1 
    real(kind=dp) :: lcgerr, lcgdev 

    ! Because the variable **array** is not used, I will comment it out 
    !real, dimension(10000) :: array 
    real(kind=dp) :: out_lcg 

    ! In the original code the variables 'lcgmean' and 'dev1' were  
    ! *undefined* before they were used in the do-loop. This will cause the 
    ! compiler to set them some random garbage values, and it will 
    ! inevitably leads to unexpected result or error in most cases. 

    ! In, order to avoid this by setting them. 
    ! For example, lcgmean = 1.0 and dev1 = 0.1 
    ! We'll then have the following: 
    lcgmean = 1.0 
    dev1 = 0.1 
    ! The above is true for the variables **sumof** 
    sumof = 0.0 
    do N = 1, 10000 
     ttl = ttl + lcg() 
     dev1 = lcg() - lcgmean 
     sumof = sumof + dev1 
    enddo 


    call muller(y1, y2) 
    call muller(y1, y2) 
    lcgmean = ttl/10000 
    lcgvar = ((sumof)**2)/10000 
    lcgdev = SQRT((sumof)**2)/10000 
    lcgerr = lcgdev/100 
    out_lcg = lcg() 
    print *, out_lcg, "mean=", lcgmean, "variance=", lcgvar, lcgerr 

end program 

나는 프로그램이 이상 정확히 원하는 것을하고는 거리가 멀다 확신 해요. 그러나 그것은 당신이 마주 친 문제를 해결합니다.

참고 : 나는 당신이 아마 새로운 생성 된 의사 라돔 번호를 액세스하려는 생각 때문에
은 내가 서브 루틴 뮬러에 Y1 및 Y2을 공급했다. 또한 변수 배열을 사용할 수있는 많은 공간을 보았습니다. 마지막으로, 어쩌면 당신이 배열의 사용을 incorporrate 수있는 방법을 참조 알고리즘 및 lcgmean, lcgvar, lcgdev 및 lcgerr와의 마지막 단계에서 수행 계산을 확인하는 것이 좋습니다 그리고이 대안은 더 효율적인지 여부 및 더 빠름