#!/usr/bin/env python # ----------------------------------------------------------------------------------- # # Python/ROOT macro for illustration of the Central Limit Theorem. # # Random numbers from different distributions are added to a sum, and many such sums # are plottet. As it turns out (and as dictated by the CLT), their distribution turns # about to be Gaussian. # The example also illutrates how widths (and therefore uncertainties) are added in # quadrature, as one has to divide the sum by the square root of the number of random # numbers that went into the sum in order to get a Gaussian of unit width (when using # random numbers of unit width). # # For more information on the Central Limit Theorem, see: # R. Barlow: page 49 # G. Cowan: page 33 # http://en.wikipedia.org/wiki/Central_limit_theorem # # Run this macro by: # From prompt: > ./CentralLimit.py # # Result in file results.root, displayed by: # In prompt: > root -l results.root # Followed by (in ROOT): TBrowser a (use mouse from there) # # Author: Troels C. Petersen (NBI) # Email: petersen@nbi.dk # Date: 22nd of August 2014 # # ----------------------------------------------------------------------------------- // from ROOT import * import math gROOT.Reset() # Set the showing of statistics and fitting results (0 means off, 1111 means all on): gStyle.SetOptStat(1111) gStyle.SetOptFit(1111) r = TRandom3() # Random number generator Nexperiments = 1000 # Number of sums produced Nuniform = 10 # Number of uniform numbers used in sum Nexponential = 0 # Number of exponential numbers used in sum Ncauchy = 0 # Number of cauchy numbers used in sum pi = 3.14159265358979323846264338328 # Selfexplanatory!!! pi = TMath.Pi() # Another way of doing it - probably better! verbose = True # Print some numbers or not? Nverbose = 10 # If so, how many? SavePlots = False # Save the plots produced to file(s)? #---------------------------------------------------------------------------------- # Make histograms (name, title, number of bins, minimum, maximum): Hist_Sum = TH1F("Hist_Sum" , "Hist_Sum" , 240, -6.0, 6.0) Hist_Uniform = TH1F("Hist_Uniform" , "Hist_Uniform" , 240, -6.0, 6.0) Hist_Exponential = TH1F("Hist_Exponential", "Hist_Exponential", 240, -6.0, 6.0) Hist_Cauchy = TH1F("Hist_Cauchy" , "Hist_Cauchy" , 240, -6.0, 6.0) #---------------------------------------------------------------------------------- # Loop over process: #---------------------------------------------------------------------------------- N3sigma = 0 # Counter for the number of produced sums, that fall outside +-3 sigma for iexp in range( Nexperiments ) : if (iexp % 100 == 0) : print "At iexp : ", iexp # Show progress! sum = 0.0 # sum is the number we are going to add random numbers to! # According to the CLT, it should be Gaussianly distributed. # Generating uniform numbers (with mean 0, and RMS of 1): for i in range( Nuniform ) : x = math.sqrt(12.0) * (r.Uniform() - 0.5) # Uniform between +-sqrt(3) sum += x Hist_Uniform.Fill(x) if verbose and iexp == 0 and i < Nverbose : print " Uniform: %7.4f"%x # Generating exponential numbers (with mean 0, and RMS of 1): for i in range( Nexponential ) : x = -math.log( r.Uniform() ) - 1.0 # Exponential starting at -1 sum += x Hist_Exponential.Fill(x) if verbose and iexp == 0 and i < Nverbose : print " Exponential: %7.4f"%x # Generating numbers according to a Cauchy distribution (1 / (1 + x^2)): for i in range( Ncauchy ) : x = math.tan(pi * (r.Uniform() - 0.5)) # Cauchy with mean 0 sum += x Hist_Cauchy.Fill(x) if verbose and iexp == 0 and i < Nverbose : print " Cauchy : %7.4f"%x sum = sum / math.sqrt( Nuniform + Nexponential + Ncauchy ) Hist_Sum.Fill(sum) if not (-3.0 < sum < 3.0) : N3sigma += 1 print N3sigma #---------------------------------------------------------------------------------- # Draw output plots with corresponding fits to the screen: #---------------------------------------------------------------------------------- # Define a "canvas", i.e. a new window on the screen. # This is again a ROOT thing (hence the "T" in TCanvas), # and you give it name, title, x-position, y-position, x-size and y-size canvas = TCanvas( "canvas", "Distribution of sums", 50, 50, 1200, 600) Hist_Sum.GetXaxis().SetRangeUser(-6.1, 6.1) # From the histogram, we get the x-axis, and set the range. # Note that I plot beyond the histogram range to include under- and overflow bins! Hist_Sum.SetMinimum(0) # We also set the minimum. Hist_Sum.GetXaxis().SetTitle("Sum") # We give the x-axis a label. Hist_Sum.GetYaxis().SetTitle("Frequency") # We give the y-axis a label. Hist_Sum.SetMarkerColor(kBlue) # Set the histogram marker color to blue. Hist_Sum.SetMarkerStyle(20) # Set the histogram marker style to filled circles. # How could you know that 20 is filled circles? Well, ROOT is well documented and quite used, # so if you write "setmarkerstyle" in Google, the first hit is: http://root.cern.ch/root/html/TAttMarker.html # Here, you can find all the options, details, etc. # We define a function in ROOT as shown below. "T" from ROOT and "F1" is a 1D function. # Again it has a name (chosen to be the same as the function name), and then one has to # specify the function expression, here a predefined ROOT function "gaus" (guess what that is!). # Finally, the function range is set: fit_Sum = TF1("fit_Sum", "gaus", -5.0, 5.0) fit_Sum.SetParameters(0.0, 1.0) # Set the starting values of [0] and [1] to 0 and 1. fit_Sum.SetLineColor(kRed) # Set the line color to red. fit_Sum.SetLineWidth(3) # Set tne line width to 4. # So now we ask that the histogram is fitted by the fitting function. # The options "L" is for Likelihood fit (ChiSquare is defalt) Hist_Sum.Fit("fit_Sum", "L") # Draw the histogram with the option "e" for showing the errors in the plot. Hist_Sum.Draw("e") if (SavePlots): canvas.SaveAs("Plot_CentralLimit.pdf") # Save plot (format follow extension name) #---------------------------------------------------------------------------------- # Draw Cauchy number distribution, if used: #---------------------------------------------------------------------------------- canvas2 = TCanvas( "canvas2", "Distribution of Cauchy numbers", 20, 80, 1200, 500) canvas2.Divide(3,1) canvas2.cd(1); Hist_Uniform.Draw() canvas2.cd(2); Hist_Exponential.Draw() canvas2.cd(3); Hist_Cauchy.Draw() #---------------------------------------------------------------------------------- # Write histograms to file: #---------------------------------------------------------------------------------- outfile = TFile("results.root", "RECREATE", "Histograms") Hist_Sum.Write() Hist_Uniform.Write() Hist_Exponential.Write() Hist_Cauchy.Write() outfile.Write() outfile.Close() # Finally, ensure that the program does not termine (and the plot disappears), before you press enter: raw_input('Press Enter to exit') # ---------------------------------------------------------------------------------- # # First make sure that you understand what the Central Limit Theorem (CLT) states! # Then, acquaint yourself with the program. Do you understand why the uniform distribution # needs to go from +-sqrt(3) in order to give a distribution with a width of one (i.e. unit) # and why you subtract one from the exponential distribution (and how this works at all)? # # Then try to see what the result of adding 10 uniform random numbers is? Which # distribution does it give (surprice!!!), and how well does it resemble it? # # # Questions: # ---------- # 0) Why is there a "/sqrt(N)" in line 93 (when summing up the various contributions to sum)? # # 1) Using the sum of 10 uniform random numbers with mean 0 and width 1, # what is the probability of obtaining a number beyond 3 sigma? # What would you expect from a true Gaussian distribution? # And what about the same question for 3.5 sigma? # Perhaps increase the number of experiments run... # # 2) Now try to add 10 uniform and 10 cauchy numbers. Does that give something Gaussian? # If not Gaussian, why do the Cauchy distribution "ruin" the Gaussian distribution? # And is this in conflict with the Central Limit Theorem? # # 3) If the original distributions which are added into a sum did not have a mean # of zero and unit width, but instead each had a different mean and width, what # would be the result be? Would the Central Limit Theorem still apply? And what # would the width of the resulting distribution be? # # 4) Try to combine 10 exponential numbers into sums 1000 times, and do a ChiSquare fit # with a Gaussian to the resulting distribution. What is the Chi2 and Ndof of this? # Is that a good fit? # # 5) How few/many uniform random numbers needs to be added, before the probability # for the sum to follow a Gaussian distribution is greater than 1% (on average) # when using 1000 sums (i.e. Nexperiments = 1000)? # # # Advanced questions: # ------------------- # 1) If one used a trunkated mean of 10 Cauchy numbers (throwing away the top and bottom e.g. 10%), # will the truncated mean of 1000 Cauchy numbers then converge to a Gaussian? # # 2) Try to add a triangular PDF, which again has mean 0 and width 1. # By which method(s) is this possible? (we'll get to that in Week2) # # # ---------------------------------------------------------------------------------- # # Key questions: # -------------- # https://docs.google.com/forms/d/1uStLkhCSoeccE53ByF0KP6EfSmT4_C0DwFX3pbbKQ8E/viewform # * Using 10 uniform numbers in the sum, what fraction of sums were beyond 3.0 and 3.5 sigma? # And what fraction should have been, if the numbers were truly Gaussian distributed? # * What is the Chi2 and Ndof of the Gaussian Chi2 fit to 10 exponentials 1000 times? # * How many uniform numbers were required in problem 5? # # ----------------------------------------------------------------------------------