program discrete implicit none integer, parameter :: NMAX = 1000 integer :: Nrandom integer :: i,id,j,ios,N,Nd,k,kmin,kmax double precision :: x(NMAX),p(NMAX),F(0:NMAX),prand(NMAX) double precision :: u,psum,pxsum,xsum ! Parameters for statistics of the means integer :: Nmean double precision, allocatable :: meanstat(:),meanstatcum(:) double precision :: dmean,mean,var,meanmin,meanmean,meanmax double precision, external :: grnd ! ! Program to analyze error of DSA-like distributions obtained with poor statistics. ! ! Assumes x points are evenly distributed, does not check for that ! ! Algorithm steps as in lecture notes 6.1 ! 01 a. Read in known shape of distribution open(10,file='dsa_exercise.dat') N=0 do N=N+1 if (N > NMAX) STOP 'Increase NMAX' read(10,fmt=*,iostat=ios) x(N),p(N) ! Assume end of file reached on io error if (ios < 0) exit enddo close(10) N=N-1 print *,'Read in ',N,' data points in distribution' ! 01 b. Calculate reference 'accurate' true mean psum=0.0; pxsum=0.0; do i=1,N psum=psum+p(i) pxsum=pxsum+x(i)*p(i) enddo print *,'True mean is',pxsum/psum; ! 02 a. Select number of counts to generate Nrandom=300; ! 02 b. Select number of distributions to generate Nd=300000; ! 02 c (in slightly modified form) ! allocate data arrays for the means Nmean=4000; dmean=(x(N)-x(1))/(Nmean-1); allocate(meanstat(0:Nmean)); allocate(meanstatcum(0:Nmean)); print *,'Using Nmean dmean',Nmean,dmean ! 03 Form cumulative distribution Fk psum=0.0; F(0)=0.0; do i=1,N psum=psum+p(i) F(i)=psum enddo ! Normalize F(i) by area = psum do i=1,N F(i)=F(i)/psum enddo ! 04 Set meanstat to 0 do i=0,Nmean meanstat(i)=0 meanstatcum(i)=0 enddo ! Initialize random number generator call sgrnd(888134) ! 05 loop over individual distributions do id=1,Nd ! ------------- Steps 06 - 14: as in exercise 2.3 --------------- ! Empty result array do i=1,N prand(i)=0 enddo ! Generate random numbers and do statistics do i=1,Nrandom ! This routine does not work right if u==0.0 ! Hence this silly little loop is needed. do u=grnd(); if (u > 0.0d0) exit enddo ! binary search to find correct index in F kmin=0; kmax=N; do k=(kmax+kmin)/2 if (F(k) < u) then kmin=k else kmax=k endif !print *,u,kmin,kmax,k,F(k) if (kmin >= kmax-1) exit enddo !print *,u,F(kmin),F(kmax) ! Now x(kmax) is the random number wanted ! but we do not actually need it in explicit form prand(kmax)=prand(kmax)+1.0 enddo ! ------------- End of steps 06 - 11: ---------------------------- ! 12 Calculate weighted mean and Gaussian error of this distribution psum=0.0; mean=0.0; do i=1,N psum=psum+prand(i); mean=mean+x(i)*prand(i); enddo mean=mean/psum var=0.0 do i=1,N var=var+prand(i)*(x(i)-mean)**2 enddo if (id < 10) then ! print this for first 10 cases print *,'mean and Gaussian error',mean,sqrt(var/Nrandom)/sqrt(1.0*Nrandom) endif ! 13 For position j sum up statistics j=nint((mean-x(1))/dmean); if (j<0 .or. j > Nmean) STOP 'Impossible j' meanstat(j)=meanstat(j)+1; !print *,'intermediate mean',mean,j,x(j) enddo ! 14 End loop over distributions ! Output the final generated distribution just as an example open(10,file='dsa_nonegative.rand') do i=1,N write(10,*) x(i),prand(i) enddo close(10) ! 15 Form the cumulative mean distribution function and normalize it xsum=0.0; meanstatcum(0)=0.0; do i=1,Nmean xsum=xsum+meanstat(i) meanstatcum(i)=xsum enddo do i=1,Nmean meanstatcum(i)=meanstatcum(i)/xsum enddo ! 16 find points where meanstatcum is 0.16 and 0.50 and and 0.84 do i=1,Nmean ! Use linear interpolation in determining exact value if (meanstatcum(i-1) < 0.16 .and. meanstatcum(i)>=0.16) then meanmin=x(1)+dmean*(i-1)+dmean*(0.16-meanstatcum(i-1))/(meanstatcum(i)-meanstatcum(i-1)) endif if (meanstatcum(i-1) < 0.5 .and. meanstatcum(i)>=0.5) then meanmean=x(1)+dmean*(i-1)+dmean*(0.50-meanstatcum(i-1))/(meanstatcum(i)-meanstatcum(i-1)); endif if (meanstatcum(i-1) < 0.84 .and. meanstatcum(i)>=0.84) then meanmax=x(1)+dmean*(i-1)+dmean*(0.84-meanstatcum(i-1))/(meanstatcum(i)-meanstatcum(i-1)); endif enddo print *,'Mean from distribution is',meanmean print *,'1 sigma limits are',meanmin,meanmax print *,'Uncertainty is + ',meanmax-meanmean,' - ',meanmean-meanmin ! Print out distributions open(11,file='meanstat.dat') open(12,file='meanstatcum.dat') do i=1,Nmean write (11,*) x(1)+dmean*i,meanstat(i) write (12,*) x(1)+dmean*i,meanstatcum(i) enddo close(11) close(12) end program discrete !************************************************************************ ! Mersenne twister random number generator !************************************************************************ subroutine sgrnd(seed) ! implicit integer(a-z) ! ! Period parameters parameter(N = 624) ! dimension mt(0:N-1) ! the array for the state particle common /block/mti,mt save /block/ ! ! setting initial seeds to mt[N] using ! the generator Line 25 of Table 1 in ! [KNUTH 1981, The Art of Computer Programming ! Vol. 2 (2nd Ed.), pp102] ! mt(0)= iand(seed,-1) do 1000 mti=1,N-1 mt(mti) = iand(69069 * mt(mti-1),-1) 1000 continue ! return end !************************************************************************ double precision function grnd() ! implicit integer(a-z) ! ! Period parameters parameter(N = 624) parameter(N1 = N+1) parameter(M = 397) parameter(MATA = -1727483681) ! constant vector a parameter(UMASK = -2147483648) ! most significant w-r bits parameter(LMASK = 2147483647) ! least significant r bits ! Tempering parameters parameter(TMASKB= -1658038656) parameter(TMASKC= -272236544) ! dimension mt(0:N-1) ! the array for the state vector common /block/mti,mt save /block/ data mti/N1/ ! mti==N+1 means mt[N] is not initialized ! dimension mag01(0:1) data mag01/0, MATA/ save mag01 ! mag01(x) = x * MATA for x=0,1 ! TSHFTU(y)=ishft(y,-11) TSHFTS(y)=ishft(y,7) TSHFTT(y)=ishft(y,15) TSHFTL(y)=ishft(y,-18) ! if(mti.ge.N) then ! generate N words at one time if(mti.eq.N+1) then ! if sgrnd() has not been called, call sgrnd(4357) ! a default initial seed is used endif ! do 1000 kk=0,N-M-1 y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK)) mt(kk)=ieor(ieor(mt(kk+M),ishft(y,-1)),mag01(iand(y,1))) 1000 continue do 1100 kk=N-M,N-2 y=ior(iand(mt(kk),UMASK),iand(mt(kk+1),LMASK)) mt(kk)=ieor(ieor(mt(kk+(M-N)),ishft(y,-1)),mag01(iand(y,1))) 1100 continue y=ior(iand(mt(N-1),UMASK),iand(mt(0),LMASK)) mt(N-1)=ieor(ieor(mt(M-1),ishft(y,-1)),mag01(iand(y,1))) mti = 0 endif ! y=mt(mti) mti=mti+1 y=ieor(y,TSHFTU(y)) y=ieor(y,iand(TSHFTS(y),TMASKB)) y=ieor(y,iand(TSHFTT(y),TMASKC)) y=ieor(y,TSHFTL(y)) ! if(y.lt.0) then grnd=(dble(y)+2.0d0**32)/(2.0d0**32-1.0d0) else grnd=dble(y)/(2.0d0**32-1.0d0) endif ! return end