用拒絕方法創建正態分佈產生錯誤的前因子

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），I認爲正態分佈，我獲得也應該給予前因子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日。這是我想要的prefactor的1/10。問題是如果我放大我整合函數的區域，這個前因子不會變得更好。這意味着如果使用此代碼創建-5000到+ 5000的分佈，那麼即使該函數的-5000到5000的積分導致我使用的分佈爲0.2，我仍然可以獲得相同的前因子。 （我把隨機分佈的值放到matlab中，並用這些值計算-5000到5000的數值積分，得到0.2，這意味着這裏積分的前因子應該是1/pi * 1/5。 ，我很困惑的事實是積分從-5000到+50000的高斯僅爲0.2，根據數學這個積分約爲1.所以有些事情是錯的）

Answer 1

我剛剛使用你的例程來在-2和2之間生成1000個點並獲得高斯分佈。

如何生成直方圖？未標準化的直方圖可以用函數N exp(-x**2)/sqrt(pi) * dx繪製，其中N是點數，dx是合併間隔。