#!/usr/bin/env python # ---------------------------------------------------------------------------------- # # # Python script for illustrating how to generate random numbers following a # specific PDF using uniformly distributed random numbers. Both the Accept-Reject # (Von Neumann) and transformation methods are used. # # Author: Troels C. Petersen (NBI) # Email: petersen@nbi.dk # Date: 27th of November 2017 # # ---------------------------------------------------------------------------------- # from __future__ import division, print_function import numpy as np import matplotlib.pyplot as plt from iminuit import Minuit # you will need these, from probfit import Chi2Regression, BinnedLH # if you want to fit plt.close('all') # function to create a nice string output def nice_string_output(names, values, extra_spacing = 2): max_values = len(max(values, key=len)) max_names = len(max(names, key=len)) string = "" for name, value in zip(names, values): string += "{0:s} {1:>{spacing}} \n".format(name, value, spacing = extra_spacing + max_values + max_names - len(name)) return string[:-2] plt.close("all") r = np.random r.seed(1) SavePlots = False # Determining if plots are saved or not Npoints = 10000 print(" Yes, you have to do some calculations and coding, before the correct results appear!") #----------------------------------------------------------------------------------- # # Problem 1: Produce random points following PDF f(x) ~ exp(-x/3), x in [0,inf] #----------------------------------------------------------------------------------- # # FIRST: # Which method should you use to produce random numbers according to this distribution? # Think about this, argue for one method and possibly discuss it with a fellow student. # NEXT: # Do so, and see that it fits, by drawing (not fitting) the function on top! Hist_exp = np.zeros(Npoints) for i in range (Npoints) : x = r.uniform() # This is of course where you should change things!!! Hist_exp[i] = x fig, ax = plt.subplots(figsize=(10, 6)) ax.set_xlabel("x (exponentially distributed)") ax.set_ylabel("Frequency") ax.set_xlim( 0.0, 20.0 ) ax.hist( Hist_exp, bins=100, range=(0.0,20.0), histtype='step', label='histogram' ) def exp_func(x) : # Ask yourself, if you can figure out the normalisation constant, knowing the number # of events produced and taking the bin width into account. normalization = 1 return normalization * np.exp(-x/3.0) func_x = np.linspace(0.0, 20.0, 1000) func_y = exp_func(func_x) ax.plot(func_x,func_y,'r-', label='function (not fitted)') names = ['Entries', 'Mean', 'STD Dev.'] values = ["{:d}".format(len(Hist_exp)), "{:.3f}".format(Hist_exp.mean()), "{:.3f}".format(Hist_exp.std(ddof=1)) ] ax.text(0.05, 0.95, nice_string_output(names, values), family='monospace', transform=ax.transAxes, fontsize=12, verticalalignment='top') ax.legend(loc='best') plt.tight_layout() plt.show(block=False) if (SavePlots) : fig.savefig("Hist_exp.pdf", dpi=600) #----------------------------------------------------------------------------------- # # Problem 2: Produce random points following f(x) ~ x*cos(x), x in [0,pi/2] #----------------------------------------------------------------------------------- # # FIRST: # Which method should you use to produce random numbers according to this distribution? # NEXT: # Do so, and see that it fits, by drawing (not fitting) the function on top! List_xcos = [] Nhit = 0 for i in range (Npoints) : x = r.uniform() # Use these to make a "box" around the range wanted. y = r.uniform() # Use these to make a "box" around the range wanted. if (y < x) : # Test if x-value is accepted or not (change yourself!) Nhit += 1 List_xcos.append(x) List_xcos = np.array(List_xcos) fig2, ax2 = plt.subplots(figsize=(10, 6)) ax2.set_xlabel("x") ax2.set_ylabel("Frequency") ax2.set_xlim( 0.0, np.pi/2.0 ) ax2.hist( List_xcos, bins=100, range=(0.0,np.pi/2.0), histtype='step', label='histogram' ) def xcos_func(x) : # Again, think about the normalisation: N = np.pi return N * x * np.cos(x) func2_x = np.linspace(0.0, np.pi/2.0, 1000) func2_y = xcos_func(func2_x) ax2.plot(func2_x, func2_y, 'r-', label='function (not fitted)') names = ['Entries', 'Mean', 'STD Dev.'] values = ["{:d}".format(len(List_xcos)), "{:.3f}".format(List_xcos.mean()), "{:.3f}".format(List_xcos.std(ddof=1)) ] ax2.text(0.03, 0.95, nice_string_output(names, values), family='monospace', transform=ax2.transAxes, fontsize=12, verticalalignment='top') ax2.legend(loc='best') plt.tight_layout() plt.show(block=False) if (SavePlots) : fig2.savefig("Hist_xcos.pdf", dpi=600) # ALSO: # Use the random numbers to calculate an integral and its uncertainty. # Integral of f(x) in the range [0,pi/2]: integral = Nhit / Npoints * np.pi/2.0 * 1.0 eintegral = -1.0 # Calculate this yourself! print(" Integral of f(x) = x*cos(x), x in [0,pi/2] is: {:7.4f} +- {:6.4f}".format(integral, eintegral)) #----------------------------------------------------------------------------------- # # Problem 3 for you: Try to combine the two methods yourself: #----------------------------------------------------------------------------------- # # # Estimate the integral of exp(-x/3.0)*cos(x)^2 in the interval [0,inf] using a # combination of Transformation and Hit-and-Miss method. Nhit = 0 for i in range ( Npoints ) : # Well, the inside of the loop is up to you! Nhit += 1 integral = -1.0 # Well, put your value here eintegral = -1.0 # ... print(" Integral of f(x) = exp(-x/3)*cos(x)^2, x in [0,inf] is: {:7.4f} +- {:6.4f}".format(integral,eintegral)) try: __IPYTHON__ except: raw_input('Press Enter to exit') # ----------------------------------------------------------------------------------- # # # The purpose of this exercise is understand how to produce random numbers both using # the Hit-and-Miss (Von Neumann) and the Transformation method, possibly to evaluate # an integral (here only in 1D, to make it easy/illustrative). # # Questions: # ---------- # 1) Find the integral of exp(-x/3.0)*cos(x)^2 in the interval [0,inf] using a # combination of the Transformation and Hit-and-Miss methods. # # ----------------------------------------------------------------------------------- #