import numpy as np import matplotlib.pyplot as plt plt.rcParams.update({'font.size': 20}) from scipy.optimize import curve_fit

xvals = np.random.exponential(scale=0.5,size=1000)

# histogram the simulated data # documentation: https://numpy.org/doc/stable/reference/generated/numpy.histogram.html BINS = np.linspace(0.5,3,20) print ('bin edges are : ',BINS) binned_data, bin_edges = np.histogram(xvals,bins=BINS)

# plot the obtained distribution, assuming sqrt(N) for the obtained measurements fig = plt.figure(figsize=(6,6)) bin_centers = 0.5*(bin_edges[1:]+bin_edges[:-1]) binned_data = binned_data.astype(float) error = np.sqrt(binned_data) plt.errorbar(bin_centers,binned_data,yerr=error,fmt='o',color='b',lw=2,label='observation') plt.xlabel('x : observable',fontsize=20) plt.ylabel('event rate',fontsize=20) plt.legend(loc='best') plt.show()

def exponential(x,a,k,b=0): return a * np.exp( - x * k) + b

popt, pcov = curve_fit(exponential,bin_centers,binned_data,p0=(1,1e-6,1))

# print fit results and their uncertainty # stored in popt and diagonal entries of pcov (which gives the covariance matrix) in the order # in which they appear in the function exponential A_fit = popt[0] A_err = np.sqrt(np.diag(pcov)[0]) lambda_fit = popt[1] lambda_err = np.sqrt(np.diag(pcov)[1]) print ('best fit A : %.02f +/- %.02f'%(A_fit,A_err)) print ('best fit lambda : %.02f +/- %.02f'%(lambda_fit,lambda_err))

# plot the data and the fit result fig = plt.figure(figsize=(6,6)) bin_centers = 0.5*(bin_edges[1:]+bin_edges[:-1]) binned_data = binned_data.astype(float) error = np.sqrt(binned_data) # fine sampling of x-axis range xvalsfine = np.linspace(0.5,3,1000) plt.errorbar(bin_centers,binned_data,yerr=error,fmt='o',color='b',lw=2,label='observation') plt.plot(xvalsfine,exponential(xvalsfine,A_fit,lambda_fit),'r--',lw=2,label='fit result') plt.plot([],[],' ',label="decay constant ={}+-{}".format(lambda_fit,lambda_err)) plt.plot([],[],' ',label="number of events ={}+-{}".format(A_fit,A_err)) plt.xlabel('x : observable',fontsize=20) plt.ylabel('event rate',fontsize=20) plt.legend(loc=(1.04,0),frameon=True) plt.show()

# create 3 datasets following a normal distribution, with mean = 0, sigma = 1. x_vals = np.random.normal(0,1.0,1000) y_vals = np.random.normal(0,1.0,1000) z_vals = np.random.normal(0,1.0,1000) # add quantites in quadrature, where q_i = (x_i^2 + y_i^2 + z_i^2) q_vals = np.square(x_vals) + np.square(y_vals) + np.square(z_vals)

# histogram the simulated data BINS = np.linspace(-10,10,25) print ('bin edges are : ',BINS) #normalize binned_data, bin_edges = np.histogram(q_vals,bins=BINS) ntotal = np.sum(binned_data) bin_centers = 0.5*(bin_edges[1:]+bin_edges[:-1]) binned_data = binned_data.astype(float) error = np.sqrt(binned_data) binned_data /= ntotal error /= ntotal

import scipy.stats as stats # plot the obtained distribution fig = plt.figure(figsize=(6,6)) plt.errorbar(bin_centers,binned_data,yerr=error,fmt='o',color='b',lw=2,label='observation') plt.xlabel('x : observable',fontsize=20) plt.ylabel('event rate',fontsize=20) plt.legend(loc=(1.04,0),frameon=True) plt.show()

import scipy.stats as stats # plot the obtained distribution fig = plt.figure(figsize=(6,6)) # plot the chi-square distributions on top of the graph, with DOF = 2,3 and 4 plt.errorbar(bin_centers,binned_data,yerr=error,fmt='o',color='b',lw=2,label='observation') plt.xlabel('x : observable',fontsize=20) plt.ylabel('event rate',fontsize=20) plt.plot(bin_centers,stats.chi2.pdf(bin_centers,df=2),color='r',lw=2,label='2 DOF') plt.plot(bin_centers,stats.chi2.pdf(bin_centers,df=3),color='g',lw=2,label='3 DOF') plt.plot(bin_centers,stats.chi2.pdf(bin_centers,df=4),color='y',lw=2,label='4 DOF') plt.legend(loc=(1.04,0),frameon=True) plt.show()

chisquare(binned_data,ddof=[2,3,4])