Создание нормального распределения с методом отклонения дает неправильный префактор

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/пи * ехр (-х^2), я думал, что нормальное распределение, что я должен также получить дают распределение с префактором 1/pi^(1/2), но это не так. Если я создаю гистограмму и подгоню эту гистограмму с нормальным распределением, я получу в качестве префактора примерно 0,11.

Как это возможно? Что я делаю не так?

EDIT: Это, как я создаю гистограмму

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 (пи) * 1/10 , Это 1/10 раз мой желаемый префактор. Проблема в том, что этот префактор не улучшается, если я увеличиваю область, в которой я интегрирую функцию. Это означает, что если создать с этим кодом распространение с -5000 до + 5000, то я все равно получаю тот же префактор, даже если интеграл от -5000 до 5000 этой функции приводит к 0,2 с использованием распределения, которое я использовал. (Я взял случайные распределенные значения и поместил их в матлаб и вычислил численный интеграл от -5000 до 5000 с этими значениями и получил 0.2. Это означает, что здесь префактор интеграла должен быть 1/pi * 1/5. , Меня озадачивает тот факт, что интеграл от -5000 до +50000 гауссова составляет всего 0,2. Согласно математике этот интеграл равен примерно 1. Так что что-то должно быть не так)

Answer 1

Я просто использовал вашу рутину для генерировать 1000 точек между -2 и 2 и получать гауссовское распределение.

Как вы создаете свою гистограмму? Ненормированную гистограмму можно построить с помощью функции N exp(-x**2)/sqrt(pi) * dx, где N - количество точек, а dx - интервал биннинга.