c f77 MonteCarlo.f -o MonteCarlo
c This program computes the uncertainty in the mean value of 15 numbers
c that do have individual uncertainties.  The uncertainty is computed
c via N Monte Carlo simulations.

      parameter (Npts=15,Ntrials=10)
      real A(Npts),A_trial(Npts),mean(Ntrials),eA(Npts)
      integer slot(Npts)
      data A/7.1,8.3,10.5,4.0,7.7,6.3,12.8,3.3,8.2,10.0,9.2,9.3,6.6,
     :       7.0,9.6/
      external gasdev
      idum=-479
      
      ! initialize how many times each slot is selected
      do j=1,Npts
       slot(j)=0
       eA(j)=0.25*A(j)
      enddo
      
      do n=1,Ntrials
       
       ! tweak all 15 data points
       do j=1,Npts
        A_trial(j)=A(j)+gasdev(idum)*eA(j)
       enddo
       
       call moment(A_trial,Npts,ave,adev,sdev,var,skew,curt)
       mean(n)=ave
             
      enddo
      
      print*
      print*
      print'(a16)','The data points:'
      print*,'7.1 8.3 10.5 4.0 7.7 6.3 12.8 3.3 8.2 10.0 9.2 9.3 6.6 ',
     :       '7.0 9.6'
      call moment(A,Npts,ave,adev,sdev,var,skew,curt)
      print 100,'have 25% uncertainty, a mean of ',
     : ave,'+/-',sdev/sqrt(real(Npts-1)),
     : ', and a standard deviation of ',sdev
 100  format(a32,f4.2,a3,f4.2,a30,f4.2)
      print*
      call moment(mean,Ntrials,ave,adev,sdev,var,skew,curt)
      print'(i7,1x,a24)',Ntrials,'Monte Carlo simulations:'
      print '(a16)','Average  Std_Dev'
      print '(1x,f4.2,5x,f4.2)',ave,sdev
      print*
            
      end ! of main program
      include '/usr/local/Numerical_Recipes/recipes_f/recipes/ran1.f'
      include '/usr/local/Numerical_Recipes/recipes_f/recipes/moment.f'
      include '/usr/local/Numerical_Recipes/recipes_f/recipes/gasdev.f'