#!/usr/bin/env python # ----------------------------------------------------------------------------------- # # # Python script for illustrating the two methods for generating random numbers according # to a known PDF, starting with random uniformly distributed numbers: # - Transformation method (if function can be integrated and then inverted). # - Hit & Miss (or Accept-Reject) method (by Ulam Stanislav and John Von Neumann). # # The function used for illustration is: f(x) = 2x, in the interval [0,1]. # # For more information see: # G. Cowan: Chapter 3 # P. R. Bevington: page 81-84 # # Author: Troels C. Petersen (NBI) # Email: petersen@nbi.dk # Date: 25th of November 2017 # # ----------------------------------------------------------------------------------- # from __future__ import division, print_function import numpy as np import matplotlib.pyplot as plt from iminuit import Minuit from probfit import Chi2Regression, pdf from scipy import stats 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] # ----------------------------------------------------------------------------------- # # # ----------------------------------------------------------------------------------- # # Random numbers r = np.random r.seed(22) SavePlots = True # Set parameters: Nnumbers = 10000 # Number of random numbers produced. # Histogram settings (note the choice of 120 bins): Nbins = 120 xmin = -0.1 xmax = 1.1 x_trans = np.zeros(Nnumbers) x_hitmiss = np.zeros(Nnumbers) #---------------------------------------------------------------------------------- # Loop over process (generate numbers according to PDF f(x) = 2x in [0,1]): #---------------------------------------------------------------------------------- for i in range(Nnumbers) : # Transformation method: # ---------------------- # Integration gives the function F(x) = x^2, which inverted gives F^-1(x) = sqrt(x): x_trans[i] = np.sqrt(r.uniform()) # ...so we let x_trans equal sqrt() of the uniform number! # Hit & Miss method: # ------------------ # Generate two random numbers uniformly distributed in [0,1]x[0,2], until they # fulfill the "Hit requirement": x_hitmiss[i] = 0.0 y = 1.0 while (2.0*x_hitmiss[i] < y) : # ...so keep making new numbers, until this is fulfilled! x_hitmiss[i] = 1.0 * r.uniform() y = 2.0 * r.uniform() Hist_Trans, Hist_Trans_edges = np.histogram(x_trans, bins=Nbins, range=(xmin, xmax)) Hist_Trans_centers = 0.5*(Hist_Trans_edges[1:] + Hist_Trans_edges[:-1]) Hist_Trans_error = np.sqrt(Hist_Trans) Hist_Trans_indexes = Hist_Trans > 0 # Produce a new histogram, using only bins with non-zero entries Hist_HitMiss, Hist_HitMiss_edges = np.histogram(x_hitmiss, bins=Nbins, range=(xmin, xmax)) Hist_HitMiss_centers = 0.5*(Hist_HitMiss_edges[1:] + Hist_HitMiss_edges[:-1]) Hist_HitMiss_error = np.sqrt(Hist_HitMiss) Hist_HitMiss_indexes = Hist_HitMiss > 0 # Produce a new histogram, using only bins with non-zero entries #---------------------------------------------------------------------------------- # Plot histograms on screen: #---------------------------------------------------------------------------------- fig, ax = plt.subplots(figsize=(14, 8)) ax.set_xlabel("Random number") ax.set_ylabel("Frequency") # Fitting histogram with Transformation numbers: chi2_object_trans = Chi2Regression(pdf.linear, Hist_Trans_centers[Hist_Trans_indexes], Hist_Trans[Hist_Trans_indexes], Hist_Trans_error[Hist_Trans_indexes]) #m_truth = (xmax-xmin)/Nbins*Nnumbers*2 minuit_trans = Minuit(chi2_object_trans, pedantic=False) # If you want to help the fit: limit_c=(-1,1), limit_m=(m_truth*0.8,m_truth*1.2) ) minuit_trans.migrad() chi2_trans = minuit_trans.fval ndof_trans = chi2_object_trans.ndof prob_trans = stats.chi2.sf(chi2_trans, ndof_trans) m, c = minuit_trans.args x_fit = np.linspace(xmin, xmax, 1000) pdf_linear_vec = np.vectorize(pdf.linear) y_fit_simple = pdf_linear_vec(x_fit,m,c) ax.plot(x_fit, y_fit_simple, 'b-') chi2_object_hitmiss = Chi2Regression(pdf.linear, Hist_HitMiss_centers[Hist_HitMiss_indexes], Hist_HitMiss[Hist_HitMiss_indexes], Hist_HitMiss_error[Hist_HitMiss_indexes]) minuit_hitmiss = Minuit(chi2_object_hitmiss, pedantic=False) minuit_hitmiss.migrad() chi2_hitmiss = minuit_hitmiss.fval ndof_hitmiss = chi2_object_hitmiss.ndof prob_hitmiss = stats.chi2.sf(chi2_hitmiss, ndof_hitmiss) m, c = minuit_hitmiss.args x_fit = np.linspace(xmin, xmax, 1000) y_fit_simple = pdf_linear_vec(x_fit,m,c) ax.plot(x_fit, y_fit_simple, 'r-') names = ['Transformation:','Entries','Mean','RMS','Chi2 / ndof','Prob','y-intercept','slope','','Hit & Miss:','Entries','Mean','RMS','Chi2 / ndof','Prob','y-intercept','slope'] values = ["", "{:d}".format(len(x_trans)), "{:.3f}".format(x_trans.mean()), "{:.3f}".format(x_trans.std(ddof=1)), "{:.3f} / {:d}".format(chi2_trans, ndof_trans), "{:.3f}".format(prob_trans), "{:.3f} +/- {:.3f}".format(minuit_trans.values['c'], minuit_trans.errors['c']), "{:.3f} +/- {:.3f}".format(minuit_trans.values['m'], minuit_trans.errors['m']), "","", "{:d}".format(len(x_hitmiss)), "{:.3f}".format(x_hitmiss.mean()), "{:.3f}".format(x_hitmiss.std(ddof=1)), "{:.3f} / {:d}".format(chi2_hitmiss, ndof_hitmiss), "{:.3f}".format(prob_hitmiss), "{:.3f} +/- {:.3f}".format(minuit_hitmiss.values['c'], minuit_hitmiss.errors['c']), "{:.3f} +/- {:.3f}".format(minuit_hitmiss.values['m'], minuit_hitmiss.errors['m']), ] ax.text(0.05, 0.95, nice_string_output(names, values), family='monospace', transform=ax.transAxes, fontsize=10, verticalalignment='top') # Drawing histograms with errors in the same plot: ax.errorbar(Hist_Trans_centers[Hist_Trans_indexes], Hist_Trans[Hist_Trans_indexes], Hist_Trans_error[Hist_Trans_indexes], fmt='b.', capsize=2, capthick=2, label="Transformation") ax.errorbar(Hist_HitMiss_centers[Hist_HitMiss_indexes], Hist_HitMiss[Hist_HitMiss_indexes], Hist_HitMiss_error[Hist_HitMiss_indexes], fmt='r.', capsize=2, capthick=2, label="Hit & Miss") ax.set_ylim(bottom=0) # Legend: ax.legend(loc='lower right') plt.tight_layout() plt.show(block=False) if (SavePlots) : fig.savefig("Hist_TransVsHitMiss.pdf", dpi=600) # Calculate the probability that these two distributions are the same: # Note: Makes sure that at least one of the two histograms is represented (i.e. non-zero) in each bin! # Finally, ensure that the program does not termine (and the plot disappears), before you press enter: try: __IPYTHON__ except: raw_input('Press Enter to exit') #---------------------------------------------------------------------------------- # # Questions: # ---------- # 1) Make sure that you understand how the two methods works! What is the # efficiency of each of the methods? # # 2) Compare the two histograms - are they from the same distribution? Answer by # first comparing the fit parameters, and see if they overlap. Afterwards, # calculate the Chi2 between the two histograms (i.e. the difference divided # by the uncertainty bin by bin), and compare this to the number of degrees # of freedom. # # 3) Try to repeat the exercise with the function f(x) = x^3 in the range [1,C], # (perhaps in a new file, so that you save this macro), where C is a number, # which ensures that the PDF f(x) is normalized. Think about what changes. # Again compare the two histograms. # #----------------------------------------------------------------------------------