// 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 <iostream>
#include <sstream>
#include <fstream>

#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; i<N; i++) {
    res[i]=1+(int)gRandom->Uniform(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;
}