//Template fuer numerische Loesung von "Maximum-Likelihood-Sch ̈tzung fur Normalverteilungen"

//Die folgende Konstante gibt die Groesse der Stichprobe an.
static const int N = 10;

//Die negative Log-Likelihood Funktion haengt ab von der Groesse der Stichprobe,
//sowie von der Standardabweichung und vom Mittelwert.
//Im 2D Fall werden Mittelwert und Standardabweichung als Parameter x uebergeben.
double NegativeLL2D(double* x, double par[N]){
  //------------------ Das folgende ist im Template modifiziert.
  float xx = x[0];
  float yy = x[1];
  double intermediate = 0.;
  for (int ii = 0; ii < N; ii++){
    intermediate += pow(par[ii] - xx, 2)/ (2. * pow(yy, 2));
  }
  return intermediate - (N/2.)*(TMath::Log(1/(2*TMath::Pi() * yy*yy)));
  //-------------------------------------------------------------
}

//1D Fall
double NegativeLLMean(double*x, double par[N+1]){
  float xx = x[0];
  double intermediate =0.;

  for (int ii = 0; ii < N; ii++){
    intermediate+= pow(par[ii] - xx, 2)/ (2. * pow(par[N], 2));
  }
  return intermediate - (N/2.)*(TMath::Log(1/(2*TMath::Pi() * par[N]*par[N])));
}


int ex13_lsg(float mean = 0., float sigma = 2){
  //Zuerst wird die Stichprobe generiert.
  float probe[N];

  //Zuerst wird nur eine Stichprobe fuer den graphischen Fall behandelt.
  for(int ii = 0; ii < N; ii++){
    probe[ii] = gRandom->Gaus(mean, sigma);
  //  cout << probe[ii] << endl;
  }

  TF2* myNegativeLL2D = new TF2("myNegativeLL2D", NegativeLL2D, mean - 5*sigma, mean  + 5*sigma, sigma / 4., sigma * 4., N);
  for(int ii =0; ii < N; ii++){
    myNegativeLL2D->SetParameter(ii, probe[ii]);
  }
  
  //Zur graphischen Analyse kann ein 2D Lego-Plot erstellt werden:
  myNegativeLL2D->SetNpx(100);
  myNegativeLL2D->SetNpy(100);
  myNegativeLL2D->Draw("lego");

  //Die Funktionalitaet, um das 2D Minimum zu bekommen ist
  //http://root.cern.ch/root/html/TF2.html#TF2:GetMinimumXY
  double minX = 100; // Das ist der Schaetzwert fuer den Mittelwert
  double minY = 100; // Das ist der Schaetzwert fuer die Standardabweichung
  myNegativeLL2D->GetMinimumXY(minX, minY);
  cout << minX << ", " << minY << " This is the Minimum " <<endl;

  //Nun die Methode mit parametrischem Bootstrapping.
  //Man auch mehr bins und eine groessere Zahl von Versuchen benutzen, aber diese Methode ist sehr rechenaufwaendig.
  TH1F* myhMean = new TH1F("myhMean", "Distribution for Estimation of Mean",25,  mean - 3.* sigma, mean + 3.*sigma);
  TH1F* myhSigma= new TH1F("myhSigma","Distribution for Estimation of Sigma", 25, sigma/4., 4.*sigma);
  double minX_tmp = 0.;
  double minY_tmp = 0.;
  for (int ii = 0; ii < 100; ii++){
    for (int jj = 0; jj < N; jj++){
      myNegativeLL2D->SetParameter(jj, gRandom->Gaus(minX, minY));
      myNegativeLL2D->GetMinimumXY(minX_tmp, minY_tmp);
      myhMean->Fill(minX_tmp);
      myhSigma->Fill(minY_tmp);
    } 
  }
  myhMean->SaveAs("Mean.root");
  myhSigma->SaveAs("Sigma.root");
  return 0;
}