這段代碼爲什麼給我不正確的答案？

Question

我正在嘗試編寫一個程序，該程序使用lcg作爲一個函數，以使用方塊計算器計算更多的隨機數。我已經獲得了lcg的工作，但使用box muller算法的函數正在給出錯誤的值。

這裏是我的代碼：

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中使用。也許這不是必需的。

Answer 2

鑑於我讀here，似乎最好是由子程序所以更換功能穆勒，它可以同時計算Y1和Y2。事實上，「塊」 *研磨機」的宗旨是基於兩個先前生成的僞隨機數X1和X2根據您程序結構。

然後，在產生其它兩個僞隨機數主程序，而不是調用穆勒爲功能，你應該在相應的位置寫呼叫穆勒稱其爲子程序。然而，使用功能而非子程序仍可能的，但要返回兩個值y1和y2，您可以返回矢量v，使v（1）= y1; 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到子程序穆勒，因爲我想，你可能要訪問新生成的僞拉多姆數。另外，我看到很多可以使用變量數組的空間。最後，也許是明智的檢查你的算法，並在最後階段進行lcgmean，lcgvar，lcgdev和lcgerr和計算看看如何incorporrate使用陣列和這種替代是否更高效和更快