// Assume 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 // // T.Dorigo, 2016 // --------------------------------------------------------------------------------- #include #include #include #include "TProfile.h" #include "TH1.h" #include "TF1.h" #include "TH1D.h" #include "TH2.h" #include "TStyle.h" #include "TCanvas.h" #include "TMath.h" #include "TROOT.h" #include "Riostream.h" #include "TMath.h" #include "TRandom.h" #include "TRandom3.h" using namespace std; void Die(int N=100) { int res[100000]; int n1=0, n2=0, n3=0; for (int i=0; iUniform(0.,6.); if (res[i]==1) { n1++; } else if (res[i]<6) { n2++; } else { n3++; } } cout << endl << " Die throwing results:" << endl; cout << " n1 = " << n1; cout << " n2 = " << n2; cout << " n3 = " << n3 << endl << endl; // 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 rms, t1, t2; double discr=b*b-4*a*c; double tmin=-1./6., tmax=1./3., tstar=0; if (discr<0) { cout << " No solution for max likelihood" << endl; } 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 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 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; } // Compute corresponding p-values cout << endl; cout << " This corresponds to 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 << endl; }