#include #include #include "TFile.h" #include "TH1D.h" #include "TProfile.h" #include "TCanvas.h" #include "TString.h" #include "TStyle.h" #include "TChain.h" #include "TH2.h" #include "TH1.h" #include "TF1.h" #include "TTree.h" #include "TKey.h" #include "Riostream.h" #include "TMath.h" #include "TRandom.h" #include "TRandom3.h" #include "TVirtualFitter.h" #include #include // This macro computes the variance of the distribution of Bootstrapped data // in the bins of a histogram of event counts drawn from a Poisson distribution. // // The Bootstrap drawing of events makes the bin contents behave as a compound Poisson // distribution, which has variance lambda*mu*(1+mu), where mu is the average number // of copies of each event drawn in the Bootstrap replica, and lambda is the expected number of // different events appearing in each bin. // void Bootstrap_variance (double Ndata=10000, int Nrep=100, double fracBoot=1.0) { // Ndata = Expectation value of number of events in original histogram // Nrep = Number of Bootstrap replicas drawn // fracBoot = fraction of Ndata drawn in Bootstrapped sets // This is the number of events in the Bootstrap replica samples double NdataB=Ndata*fracBoot; const int Nbins = 100; // We fix the number of bins to 100, but this is inessential if (Ndata>100000) { cout << "Too much data per sample, reduce to <100000. Exiting..." << endl; return; } // Repeat many times to get average of variance over replicas double sumvar =0; double Average_content=0; for (int i=0; iPoisson(Ndata); for (int j=0; jUniform(0.,(double)Nbins); data[j]= x; } // Create Bootstrap sample double Bdata[100000]; thisdata = gRandom->Poisson(NdataB); Average_content+=thisdata; for (int j=0; jUniform(0.,Ndata); if (index==Ndata) index=Ndata-1; // need this as Uniform can return max value Bdata[j]=data[index]; } // Study statistical properties of Bdata in each bin by computing the bin-by-bin variances int Contents[Nbins]; double sum=0; double sum2=0; for (int k=0; k=(double)k && Bdata[j]<(double)k+1.) { Contents[k]++; } } sum+= Contents[k]; sum2+= Contents[k]*Contents[k]; } double var = sum2/Nbins-pow(sum/Nbins,2); sumvar +=var; } Average_content = Average_content/Nrep; double Average_variance = sumvar/Nrep; cout << endl; cout << " Average variance in bootstrapped sets is " << Average_variance << endl; cout << " Expectation for compound Poisson is " << NdataB/Nbins*(1+Average_content/Ndata) << endl; cout << " Variance for a Poisson distribution is " << NdataB/Nbins << endl; cout << endl; }