fortran - 使用拒绝方法创建正态分布会产生错误的前因子

Question

我正在尝试使用拒绝方法从 Fortran 中的均匀分布值创建正态分布。它实际上或多或少工作得很好，但我没有得到我想要的结果。

我用这段代码生成正态分布

    function generator result(c)
            implicit none
            integer, dimension(2) :: clock
            double precision :: c,d
            call System_clock(count=clock(1))
            call random_seed(put=clock)
            !initialize matrix with random values
            call random_number(c)
    end function

    subroutine Rejection(aa,bb,NumOfPoints)
            implicit none
            double precision :: xx, yy, cc
            integer :: ii, jj, kk
            integer, intent(in) :: NumOfPoints
            double precision, intent(in) :: aa, bb
            cc=1
            xx=generator()
            allocate(rejectionArray(NumOfPoints))
            do ii=1, NumOfPoints

            call random_number(xx)
            xx=aa+(bb-aa)*xx
            call random_number(yy)
                    do while(cc*yy>1/sqrt(pi)*exp(-xx**2))
                            call random_number(xx)
                            xx=aa+xx*(bb-aa)
                            call random_number(yy)
                    end do
                    rejectionArray(ii)=xx
            end do
    end subroutine

由于我使用函数 1/pi *exp(-x^2)，所以我认为我获得的正态分布也应该给出前置因子 1/pi^(1/2) 的分布，但它才不是。如果我创建一个直方图并用正态分布拟合该直方图，我会得到大约 0.11 作为前置因子。

这怎么可能？我究竟做错了什么？

编辑:这就是我创建直方图的方式

    implicit none
    double precision ::  aa, bb
    integer :: NumOfPoints, ii, kk, NumOfBoxes, counter, CounterTotal,counterTotal2
    logical :: exists
    character(len=15) :: frmat
    double precision :: Intermediate
    %read NumOfPoints (Total amount of random numbers), NumOfBoxes
    %(TotalAmountofBins)
    open(unit=39, action='read', status='old', name='samples.txt')
            read(39,*) NumOfPoints, aa, bb, NumOfBoxes
    close(39)
    % number of Counts will be stored temporarily in 'counter'
    counter=0
    open(unit=39, action='write', status='replace', name='distRejection.txt')
    Call Rejection(aa,bb,NumOfPoints)

    do ii=1, NumOfBoxes
            counter=0
            %calculate the middle of the bin
            Intermediate=aa+(2*ii-1)*((bb-aa)/NumOfBoxes)/2
             %go through all the random numbers and check if they are within
             % one of the bins. If they are in one bin -->increase Counter
             % by one
            do kk=1, size(rejectionArray,1)
                    if(abs(RejectionArray(kk)-intermediate).le.((bb-aa)/NumOfBoxes/2)) then
                            counter=counter+1
                   end if
            end do
            %save Points + relative number of Counts in file
            write(39,100)intermediate,dble( counter)/dble(NumOfPoints)
            100 format (f10.3,T20,f10.3,/)

    end do
    close(39)

这是我获得的直方图:

现在的前置因子是 0.056，即 1/sqrt(pi)*1/10。这是我期望的前置因子的 1/10 倍。问题是，如果我扩大集成函数的区域，这个前置因素不会变得更好。这意味着，如果使用此代码创建从 -5000 到 + 5000 的分布，那么我仍然获得相同的前置因子，即使该函数从 -5000 到 5000 的积分导致我使用的分布为 0.2。 (我把随机分布的值放入matlab中，用这些值计算从-5000到5000的数值积分，得到0.2。这意味着这里积分的前置因子应该是1/pi*1/5。除此之外，我对高斯从 -5000 到 +50000 的积分只有 0.2 感到困惑。根据 mathematica，这个积分大约为 1。所以一定有问题)

Answer 1

我刚刚使用您的例程生成了 -2 和 2 之间的 1000 个点，并获得了高斯分布。

如何生成直方图？可以使用函数 N exp(-x**2)/sqrt(pi) * dx 绘制非标准化直方图其中 N 是点数，dx 是分箱间隔。