#include <string>

const int N = 10;


void ml_poisson()
{
  double nu_true = 4.5;
  double n[N] = { 0 };
  char dummy;
  TF1* flh  = new TF1( "flh", poisson_likelihood, 0, 10, N ); 
  TF1* fllh = new TF1( "fllh", poisson_neg_log_likelihood, 0, 10, N ); 

  TLatex* txt = 0;
  for( int i = 0; i < N; ++i ) {

    if( txt != 0 ) delete txt;

    // draw Poisson random number
    n[i] = static_cast<double>( gRandom->Poisson( nu_true ) );
    cout << "Zufallszahl #" << i+1 << ":\t" << n[i] << endl;
    flh->SetParameters( n );
    //cout << "ML-Schaetzwert fuer nu: " << flh->GetMaximumX( 0, 10 ) << endl;

    fllh->SetParameters( n );
    double estimator = fllh->GetMinimumX( 0, 10 );
    cout << "ML-Schaetzwert fuer nu: " << estimator << endl;

    fllh->Draw();
    fllh->GetHistogram()->SetXTitle( "#nu" );
    fllh->GetHistogram()->SetYTitle( "-ln L(#nu)" );
    
    txt = new TLatex();
    txt->SetTextAlign( 22 );
    txt->DrawLatex( 5, 0.9*fllh->GetHistogram()->GetMaximum(), 
		    Form( "Zufallszahl: %2.0f -- Schaetzwert: %f", 
			  n[i], estimator ) );

    gPad->Modified();
    gPad->Update();

    // wait for click into canvas
    gPad->WaitPrimitive();
  }

}

double poisson_likelihood( double* x, double* par )
{
  double nu = x[0];
  double val = 1;
  for( int i = 0; i < N; ++i ) {
    if( par[i] < 1e-6 ) continue;
    val *= pow( nu, par[i] ) / TMath::Factorial( int(par[i]) ) * exp( -nu );
  }

  return val;
}

double poisson_neg_log_likelihood(double* x, double* par )
{
  double nu = x[0];
  if( nu < 0 ) return 0;

  double val = 0;
  for( int i = 0; i < N; ++i ) {
    if( par[i] < 1e-6 ) continue;
    val -= par[i] * log( nu );
    val += nu;
  }
  return val;
}