这段代码为什么给我不正确的答案？

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使用阵列和这种替代是否更高效和更快