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

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。Besieds ，令我感到困惑的是，高斯从 -5000 到 +50000 的积分只有 0.2。根据 mathematica，这个积分大约为 1。所以一定是有问题）

Answer 1

刚刚用你的例程在-2和2之间生成了1000个点，得到了高斯分布。

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

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

Creating a normal distribution with rejection method yields wrong prefactor

fortran

normal-distribution