#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 estimates the load on a die from the result of repeated throws // The load is computed from the likelihood of the die throw results // // T.Dorigo, 2016 // ------------------------------------------------------------------------------------------------ void Die2(int N=100, int Nrep=100) { // N = number of die throws // Nrep = number of repetitions of the experiment if (N>100000) return; // do not allow more than 100000 throws per exp double bignumber=100000.; // number used to avoid division by zero in calculation of error below double ok=0; // this stores the number of repetitions when the load on the die is compatible with the true value // Start of cycle on repetitions of the experiment for (int rep=0; repUniform(0.,6.); if (res[i]==1) { n1++; } else if (res[i]<6) { n2++; } else { n3++; } } /* cout << "Die throwing" << endl; */ /* cout << "Results:" << endl; */ /* cout << "n1 = " << n1 << endl; */ /* cout << "n2 = " << n2 << endl; */ /* cout << "n3 = " << n3 << endl; */ /* cout << endl; */ // Hypothesize we have a bias t in the probability of the outcome of // the die throw: // p(1) = 1/6-t/2; // p(2)... p(5) = 1/6-t/8; // p(6) = 1/6+t; // let us fit for t with a explicit calculation of the max likelihood // Quadratic equation for max of L: coefficients double a = 18*(n1+n2+n3); double b = -3*(7*n1+n2+10*n3); double c = -(4*n1+n2-8*n3); double t1,t2; double discr=b*b-4*a*c; // Lower and upper bounds on t are determined by the fact that 0<=p(i)<=1 (see above) double tmin=-1./6.; double tmax=1./3.; // Tstar is estimate of bias from this experiment double tstar=0; double rms=0; if (discr<0) { cout << " No solution for max likelihood" << endl; // This never happens in the considered case } else { t1 = (-b-sqrt(discr))/(2*a); t2 = (-b+sqrt(discr))/(2*a); if (t1>=tmin && t1<=tmax) { if (t2>=tmin && t2<=tmax) { // NNBB when n1=0 there is always one solution at t=0.33333 //cout << "Two solutions ?? Picking one away from boundary" << endl; if (t1-tmin>tmax-t2) { cout << " Bias is estimated to be t = " << t1; tstar=t1; } else { cout << " Bias is estimated to be t = " << t2; tstar=t2; } } else { cout << " Bias is estimated to be t = " << t1; tstar=t1; } } else if (t2>=tmin && t2<=tmax) { cout << " Bias is estimated to be t = " << t2; tstar=t2; } // Determine error from inverse of second derivative of logL (Rao-Cramer-Frechet bound) double d2logl; if (tstar>-1/6. && tstar<1/3.) { d2logl=9*n1/pow(1-3*tstar,2)+ 9*n2/pow(4-3*tstar,2)+ 36*n3/pow(1+6*tstar,2); rms = sqrt(1/d2logl); } cout << " +- " << rms << endl; } // turn into p-values /* cout << endl; */ /* cout << "This means the following probabilities:" << endl; */ /* cout << "p(1) = " << 1./6.-tstar/2. << " +- " << rms/2. << endl; */ /* cout << "p(x) = " << 1./6.+tstar/8. << " +- " << rms/8. << endl; */ /* cout << "p(6) = " << 1./6.+tstar << " +- " << rms << endl; */ if (fabs(tstar)