Creating a normal distribution with rejection method yields wrong prefactor

Question

I'm trying to create normal distribution out of evenly distributed values in Fortran with the rejection method. It actually works more or less well, but I don't get exactly the result that I want.

I generate the normal distribution with this segment of code

    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

Since I am using as function 1/pi *exp(-x^2), I thought that the normal distribution that I obtain should also give a distribution with the prefactor 1/pi^(1/2), but it does not. If I create a histogram and fit this histogram with a normal distribution, I obtain as prefactor approximately 0.11.

How is this possible? What am I doing wrong?

EDIT: This is how I create the histogram

    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)

This is what I obtain as histogram:

The prefactor now is 0.056 what is 1/sqrt(pi)*1/10. This is 1/10 times my desired prefactor. The Problem is that this prefactor does not get any better if I enlarge the Region over which I integrate the function. This means that if create with this Code a Distribution from -5000 to + 5000 then I still obtain the same prefactor even though the integral from -5000 to 5000 of this function leads to 0.2 with the Distribution that I used. (I took the randomly distributed values and put them into matlab and calculated the numerical integral from -5000 to 5000 with those values and obtained 0.2. This means that here the prefactor of the integral should be 1/pi*1/5. Besieds that, I am puzzled by the fact that the Integral from -5000 to +50000 of a gaussian is only 0.2. According to mathematica this integral is approximately 1. So something has got to be wrong)

Answer 1

I just used your routine to generate 1000 points between -2 and 2 and obtained a gaussian distribution.

How do you generate your histogram? The un-normalized histogram can be plotted with the function N exp(-x**2)/sqrt(pi) * dx where N is the number of points and dx is the binning interval.