#!/usr/bin/env python # ---------------------------------------------------------------------------------- # # # Python/ROOT macro for illustrating the principle of maximum likelihood and a # likelihood fit. # # This is both an exercise, but also an attempt to illustrate four things: # - How to make a (binned and unbinned) Likelihood function/fit. # - The difference and a comparison between a Chi-square and a (binned) Likelihood. # - The difference and a comparison between a binned and unbinned Likelihood. # - What goes on behind the scenes in ROOT, when it is asked to fit something. # # In this respect, the exercise is more of an illustration rather than something to # be used directly, which is why it is followed by another exercise, where you can # test if you have understood the differences, and when to apply which fit method. # # The example uses 50 exponentially distributed random times, with the goal of # finding the best estimate of the lifetime (data is generated with lifetime, tau = 1). # Three estimates are considered: # - Chi-square fit (chi2) # - Binned Likelihood fit (bllh) # - Unbinned Likelihood fit (ullh) # # While the mean can be simply calculated, the three other methods are based on a scan # of values for tau in the range [0.5, 2.0]. For each value of tau, the chi2, bllh, and # llh are calculated (in the two likelihood cases, it is actually -2*log(likelihood) # which is calculated). # # Note that the unbinned likelihood is in principle the "optimal" fit, but also the # most difficult, since you are required to write the likelihood function yourself! # However, the two other constructions may be essentially the same, when there is # enough statistics (i.e. errors are Gaussian). # # Author: Troels C. Petersen (NBI) # Email: petersen@nbi.dk # Date: 27th of November 2016 # # ---------------------------------------------------------------------------------- # from ROOT import * # Import ROOT libraries (to use ROOT) from array import array # Import array for TGraphErrors... import math # Import MATH libraries (for extended MATH functions) gStyle.SetOptStat("emr") gStyle.SetOptFit(1111) # ---------------------------------------------------------------------------------- # def sqr(a) : return a*a # ---------------------------------------------------------------------------------- # # Parabola: # Note that the parabola is defined differently than normally. The Parameters are: # par[0]: Minimum value (i.e. constant) # par[1]: Minimum position (i.e. x of minimum) # par[2]: Quadratic term. # ---------------------------------------------------------------------------------- # def func_para(x, par) : return par[0] + par[2]*sqr(x[0]-par[1]) # ---------------------------------------------------------------------------------- # # Double parabola with different slopes on each side of the minimum: # Parameters are as follows: # par[0]: Minimum value (i.e. constant) # par[1]: Minimum position (i.e. x of minimum) # par[2]: Quadratic term on lower side # par[3]: Quadratic term on higher side # ---------------------------------------------------------------------------------- # def func_asympara(x, par) : if (x[0] < par[1]) : return par[0] + par[2]*sqr(x[0]-par[1]) else : return par[0] + par[3]*sqr(x[0]-par[1]) # ---------------------------------------------------------------------------------- # # Likelihood program: # ---------------------------------------------------------------------------------- # SavePlots = False # Determining if plots are saved or not verbose = True # Should the program print or not? veryverbose = False # Should the program print a lot or not? ScanChi2 = False # In addition to fit for minimum, do a scan... # Parameters of the problem: Ntimes = 50; # Number of time measurements. tau_truth = 1.0; # We choose (like Gods!) the lifetime. t = [] # List of times r = TRandom3() # Random numbers r.SetSeed(1) # We set the numbers to be random, but the same for all Nbins = 50 # Number of bins in histogram tmax = 10.0 # Maximum time in histogram dt = tmax / Nbins # Time step Hist_Exp = TH1F("Hist_Exp", "Hist_Exp;Decay time (s);Frequency", Nbins, 0.0, tmax) # ------------------------------------ # Generate data: # ------------------------------------ for i in range (Ntimes) : t.append(r.Exp(tau_truth)) # Exponential with lifetime tau. Hist_Exp.Fill(t[i]) if (veryverbose) : print " %2d: t = %5.2f"%(i,t[i]) # ------------------------------------ # Analyse data: # ------------------------------------ # Define the number of tau values and their range to test in Chi2 and LLH: # (As we know the "truth", namely tau = 1, the range [0.5, 1.5] seems fitting # for the mean. The number of bins can be increased at will...) Ntau_steps = 50 min_tau = 0.5 max_tau = 2.5 delta_tau = (max_tau-min_tau) / Ntau_steps # Loop over hypothesis for the value of tau and calculate Chi2 and LLH: # --------------------------------------------------------------------- chi2_minval = 999999.9 # Minimal Chi2 value found chi2_minpos = 0.0 # Position (i.e. time) of minimal Chi2 value bllh_minval = 999999.9 bllh_minpos = 0.0 ullh_minval = 999999.9 ullh_minpos = 0.0 tau = array( 'f', [0.0]*(Ntau_steps+1) ) chi2 = array( 'f', [0.0]*(Ntau_steps+1) ) bllh = array( 'f', [0.0]*(Ntau_steps+1) ) ullh = array( 'f', [0.0]*(Ntau_steps+1) ) # Now loop of POSSIBLE tau estimates: for itau in range (Ntau_steps+1) : tau_hypo = min_tau + itau*delta_tau # Scan in values of tau tau[itau] = tau_hypo # Calculate Chi2 and binned likelihood (from loop over bins in histogram): chi2[itau] = 0.0 bllh[itau] = 0.0 for ibin in range (1, Nbins) : # Note: The number of EXPECTED events is the intergral over the bin! xlow_bin = Hist_Exp.GetBinCenter(ibin) - dt/2.0 xhigh_bin = Hist_Exp.GetBinCenter(ibin) + dt/2.0 # Given the start and end of the bin, we calculate the INTEGRAL over the bin, # to get the expected number of events in that bin: nexp = Ntimes * (exp(-xlow_bin/tau_hypo) - exp(-xhigh_bin/tau_hypo)) # The observed number of events... that is just the data! nobs = Hist_Exp.GetBinContent(ibin) if (nobs > 0) : # For ChiSquare but not LLH, we need to require Nobs > 0 chi2[itau] += sqr(nobs-nexp) / nobs # Chi2 summation/function bllh[itau] += -2.0*log(TMath.Poisson(int(nobs), nexp)) # Binned LLH function if (veryverbose and itau == 0) : print " Nexp: %10.7f Nobs: %3.0f Chi2: %5.1f BLLH: %5.1f"%(nexp, nobs, chi2[itau], bllh[itau]) # Calculate Unbinned likelihood (from loop over events): ullh[itau] = 0.0 for iev in range ( Ntimes ) : ullh[itau] += -2.0*log(1.0/tau_hypo*exp(-t[iev]/tau_hypo)) # Unbinned LLH function if (verbose) : print " %3d: tau = %4.2f chi2 = %6.2f log(bllh) = %6.2f log(ullh) = %6.2f"%(itau, tau_hypo, chi2[itau], bllh[itau], ullh[itau]) # Search for minimum values of chi2, bllh, and ullh: if (chi2[itau] < chi2_minval) : chi2_minval = chi2[itau] chi2_minpos = tau_hypo if (bllh[itau] < bllh_minval) : bllh_minval = bllh[itau] bllh_minpos = tau_hypo if (ullh[itau] < ullh_minval) : ullh_minval = ullh[itau] ullh_minpos = tau_hypo print " Decay time of minimum found: chi2 = %7.4f s bllh = %7.4f s ullh = %7.4f s"%(chi2_minpos, bllh_minpos, ullh_minpos) # ------------------------------------ # Plot and fit results: # ------------------------------------ # Define range around minimum to be fitted: min_fit = 0.30 max_fit = 0.40 # ChiSquare: # ---------- c_chi2 = TCanvas("c_chi2", "ChiSquare", 80, 40, 600, 450) graph_chi2 = TGraphErrors(Ntau_steps+1, tau, chi2) fit_chi2 = TF1("fit_chi2", func_asympara, chi2_minpos-min_fit, chi2_minpos+max_fit, 4) fit_chi2.SetParameters(chi2_minval, chi2_minpos, 20.0, 20.0) fit_chi2.SetParNames("Minimum value", "Minimum position", "Quadratic Term (low)", "Quadratic Term (high)") fit_chi2.SetLineColor(kBlue) graph_chi2.SetTitle("ChiSquare") graph_chi2.GetXaxis().SetTitle("Value of #tau") graph_chi2.GetYaxis().SetTitle("Value of ChiSquare") graph_chi2.SetMaximum(25.0) graph_chi2.SetMarkerStyle(20) graph_chi2.SetMarkerSize(0.4) graph_chi2.Fit("fit_chi2","RQ") # Draw lines corresponding to minimum Chi2 and minimum Chi2+1, and draw errors: minval = fit_chi2.GetParameter(0) minpos = fit_chi2.GetParameter(1) err_lower = 1.0 / sqrt(fit_chi2.GetParameter(2)) err_upper = 1.0 / sqrt(fit_chi2.GetParameter(3)) # Define the range in the y-axis to draw in: graph_chi2.GetYaxis().SetRangeUser(minval-2.0, minval+15.0) graph_chi2.Draw("AP") line_min = TLine(minpos - err_lower, minval, minpos + err_upper, minval) line_min.Draw() line_min1 = TLine(minpos - err_lower, minval+1.0, minpos + err_upper, minval+1.0) line_min1.Draw() line_mid = TLine(minpos, minval+1.0, minpos, minval-2.0) line_mid.Draw() line_low = TLine(minpos-err_lower, minval+1.0, minpos-err_lower, minval-2.0) line_low.Draw() line_upp = TLine(minpos+err_upper, minval+1.0, minpos+err_upper, minval-2.0) line_upp.Draw() print " Chi2 fit gives: tau = %6.4f + %6.4f - %6.4f"%(minpos, err_lower, err_upper) c_chi2.Update() if (SavePlots) : c_chi2.SaveAs("FitMinimum_chi2.pdf") # Go through tau values to find minimum and +-1 sigma: # This assumes knowing the minimum value, and Chi2s above Chi2_min+1 if (ScanChi2) : if (((chi2[0] - chi2_minval) > 1.0) and ((chi2[Ntau_steps] - chi2_minval) > 1.0)) : found_lower = False found_upper = False for itau in range (Ntau_steps+1) : if ((not found_lower) and ((chi2[itau] - chi2_minval) < 1.0)) : tau_lower = tau[itau] found_lower = True if ((found_lower) and (not found_upper) and ((chi2[itau] - chi2_minval) > 1.0)) : tau_upper = tau[itau] found_upper = True print " Chi2 scan gives: tau = %6.4f + %6.4f - %6.4f"%(chi2_minpos, chi2_minpos-tau_lower, tau_upper-chi2_minpos) else : print "Error: Chi2 values do not fulfill requirements for finding minimum and errors!" # Binned Likelihood: # ------------------ c_bllh = TCanvas("c_bllh", "Binned Likelihood", 110, 60, 600, 450) graph_bllh = TGraphErrors(Ntau_steps+1, tau, bllh) fit_bllh = TF1("fit_bllh", func_asympara, bllh_minpos-min_fit, bllh_minpos+max_fit, 4) fit_bllh.SetParameters(bllh_minval, bllh_minpos, float(Ntimes), float(Ntimes)) fit_bllh.SetParNames("Minimum value", "Minimum position", "Quadratic Term (low)", "Quadratic Term (high)") fit_bllh.SetLineColor(kBlue) graph_bllh.SetTitle("Binned Likelihood") graph_bllh.GetXaxis().SetTitle("Value of #tau") graph_bllh.GetYaxis().SetTitle("Value of Binned Likelihood") graph_bllh.SetMarkerStyle(20) graph_bllh.SetMarkerSize(0.4) graph_bllh.Fit("fit_bllh","RQ") # Define the range in the y-axis to draw in: minval = fit_bllh.GetParameter(0) graph_bllh.GetYaxis().SetRangeUser(minval-2.0, minval+15.0) graph_bllh.Draw("AP") c_bllh.Update() if (SavePlots) : c_bllh.SaveAs("FitMinimum_bllh.pdf") # Unbinned Likelihood: # -------------------- c_ullh = TCanvas("c_ullh", "Unbinned Likelihood", 140, 80, 600, 450) graph_ullh = TGraphErrors(Ntau_steps+1, tau, ullh) fit_ullh = TF1("fit_ullh", func_asympara, ullh_minpos-min_fit, ullh_minpos+max_fit, 4) fit_ullh.SetParameters(ullh_minval, ullh_minpos, float(Ntimes), float(Ntimes)) fit_ullh.SetParNames("Minimum value", "Minimum position", "Quadratic Term (low)", "Quadratic Term (high)") fit_ullh.SetLineColor(kBlue) graph_ullh.SetTitle("Unbinned Likelihood") graph_ullh.GetXaxis().SetTitle("Value of #tau") graph_ullh.GetYaxis().SetTitle("Value of Unbinned Likelihood") graph_ullh.SetMarkerStyle(20) graph_ullh.SetMarkerSize(0.4) graph_ullh.Fit("fit_ullh","RQ") # Define the range in the y-axis to draw in: minval = fit_ullh.GetParameter(0) graph_ullh.GetYaxis().SetRangeUser(minval-2.0, minval+15.0) graph_ullh.Draw("AP") c_ullh.Update() if (SavePlots) : c_ullh.SaveAs("FitMinimum_ullh.pdf") # Fit the data using ROOT (both chi2 and binned likelihood fits), # based on the build in fits in ROOT (i.e. "Hist.Fit()): # --------------------------------------------------------------- c1 = TCanvas("c1", "Fit overlayed on data", 180, 110, 600, 450) Hist_Exp.SetXTitle("Lifetime [s]") Hist_Exp.SetYTitle("Frequency [ev/0.1s]") # Note this fitting function: # The constant for normalization is NOT left free, in order to have only ONE parameter! rootfit_chi2 = TF1("rootfit_chi2", "[1] * 0.20 / [0] * exp(-x/[0])", 0.0, 10.0) rootfit_chi2.SetParameters(1.0, float(Ntimes)) rootfit_chi2.FixParameter(1, float(Ntimes)) rootfit_chi2.SetLineColor(kBlue) Hist_Exp.Fit("rootfit_chi2","RE") # The constant for normalization is NOT left free, in order to have only ONE parameter! rootfit_bllh = TF1("rootfit_bllh", "[1] * 0.20 / [0] * exp(-x/[0])", 0.0, 10.0) rootfit_bllh.SetParameters(1.0, float(Ntimes)) rootfit_bllh.FixParameter(1, float(Ntimes)) rootfit_bllh.SetLineColor(kRed) Hist_Exp.Fit("rootfit_bllh","RLE+") Hist_Exp.SetLineWidth(2) Hist_Exp.Draw("e") # Legend: leg = TLegend( 0.45, 0.35, 0.89, 0.55 ) leg.SetFillColor(kWhite) leg.SetLineColor(kWhite) leg.AddEntry(rootfit_chi2, " ChiSquare fit", "L") leg.AddEntry(rootfit_bllh, " Binned Likelihood fit", "L") leg.Draw() c1.Update() if (SavePlots) : c1.SaveAs("ExponentialDist_Fitted.pdf") raw_input( ' ... Press enter to exit ... ' ) # ---------------------------------------------------------------------------------- # # # Make sure that you understand how the likelihood is different from the ChiSquare, # and how the binned likelihood is different from the unbinned. If you don't do it, # this exercise, and much of the course and statistics in general will be a bit lost # on you! :-) # # The binned likelihood resembels the ChiSquare a bit, only the evaluation in each bin # is different, especially if the number of events in the bin is low, as the PDF # considered (Poisson for the LLH, Gaussian for the ChiSquare) is then different. # # The unbinned likelihood uses each single event, and is thus different at its core. # This can make a different, if there are only few events and/or if each event has # several attributes, which can't be summarized in a simple histogram with bins. # # # Questions: # ---------- # 1) Consider the four plots produced by running the program, starting with the one # showing the data and two fits (Chi2 and BLLH) to it. Do the two methods give # similar results, or are they different? Consider both central value and # uncertainty. # Now consider the three curves showing the chi2, bllh, and ullh as a function # of lifetime. Do you see the relation between the curves and the fit result? # Again, consider both the central values and the uncertainty. # # 2) Try to decrease the number of exponential numbers you consider to say 10, # and see how things change. Does the difference between Chi2, BLLH, and ULLH # get bigger or not? # # 3) Try to increase the number of exponential numbers you consider to say 5000, # and see what happens to the difference between Chi2 and BLLH? # Also, does the errors become more symetric? Perhaps you will need to consider # a shorter range of the fit around the mimimal value, and have to also # increase the number of points you calculate the chi2/bllh/ullh (or decrease # the range you search!), and possibly change the ranges of your plotting. # # 4) Make (perhaps in a new program) a loop over the production of random data, # and try to see, if you can print (or plot) the Chi2 and BLLH results for each # turn. Can you spot any general trends? I.e. is the Chi2 uncertainty always # lower or higher than the (B/U)LLH? And are any of the estimators biased? # # Advanced Questions: # # 5) Make a copy of the program and put in a different PDF (i.e. not the exponential). # Run it, and see if the errors are still asymetric. For the function, try either # e.g. a Uniform or a Gaussian. # NOTE: If you want to try other more complex functions, then the method # "TH1::FillRandom" might come in handy! Look it up! We'll get to this. # # ---------------------------------------------------------------------------------- #