//////////////////////////////////////////////////////////////////////
// an example code for esttimating the value of mean and sigma
// of a gaussian distribution from data, using MCMC for likelihood
// maximization  
// written for THEP SERC 2015
//////////////////////////////////////////////////////////////////////

void fit_gaus(int np=1000){
  const int N = 10;
  const int chain_size = np;

  double xg[N],mean=0.0; //array to store the data, N normal distributed random numbers

  double mu_true = 5.0, sigma_true = 1.0, mu_max = 0.0, sigma_max = 0.0;
  double ll_max=-10000000;

  //generate data. N gaussian distributed random numbers
  for (int ii = 0; ii < N; ii++){
    xg[ii]=gRandom->Gaus(mu_true, sigma_true);
    mean +=xg[ii];  

    cout<<ii<<"th data is "<<xg[ii]<<endl;
  }

  //cout<<"llmax in the beginning is "<<llmax<<endl;

  //cout<<mean/double(N)<<endl;

  double mu0, mu1=0., sigma0, sigma1=0.;
  //range of search: mu in (4.,6.), sigma in (0.5, 1.5)
  double mu_low = 4.0, mu_up = 6.0;
  double sigma_low = 0.5, sigma_up = 1.5;

  //widths of the components of the proposal function Q
  double del_mu = (mu_up - mu_low)/20.; //5% of the range of mu
  double del_sigma = (sigma_up - sigma_low)/20.; //5% of the range of mu

  // step 1: throw a random point in this range, this is the start point theta0 = (mu0,sigma0)
  mu0    = gRandom->Uniform(mu_low, mu_up);
  sigma0 = gRandom->Uniform(sigma_low, sigma_up);

  // step 2: calculate the likelihood function p(D|theta0) for the point theta0 = (mu0,sigma0)
  double ll0 = 0;
  for (int ii=0; ii<N; ii++)
    ll0 += -log(sigma0) - 0.5*log(2.*3.1415) + (-0.5*(xg[ii] - mu0)*(xg[ii] - mu0)/(sigma0*sigma0));

  //create vectors to store information of the chain
  int            npoint = 1; //number of successful moves
  vector<double> freq_theta;
  vector<double> mu;
  vector<double> sigma;
  vector<double> ll;

  freq_theta.push_back(1.);
  mu.push_back(mu0);
  sigma.push_back(sigma0);
  ll.push_back(ll0);

  for (int jj = 0; jj < chain_size; jj++){

    // step 3: generate a new random point using the proposal function
    mu1    = gRandom->Gaus(mu0,del_mu);
    //if the increment is outside the range of search in mu, generate again
    while ( mu1 > mu_up || mu1 < mu_low ) 
      mu1    = gRandom->Gaus(mu0,del_mu);

    //similarly take a step in sigma
    sigma1 = gRandom->Gaus(sigma0, del_sigma);
    while (sigma1 > sigma_up || sigma1 < sigma_low)
      sigma1 = gRandom->Gaus(sigma0, del_sigma);

    //now the new point is theta1 = (mu1,sigma1) = (mu0+del_mu,sigma0+del_sigma)

    
    // step 4: calculate p(D|theta1), actually the log of it 
    double ll1 = 0;
    for (int ii=0; ii<N; ii++)
      ll1 += -log(sigma1) - 0.5*log(2.*3.1415) + (-0.5*(xg[ii] - mu1)*(xg[ii] - mu1)/(sigma1*sigma1));
    
    // step 5: calculate the acceptance condition:
    //    if the likelihood of the new point is more than the old point, accept it 
    //    otherwise 
    //    accept it with a probability equal to the ratio of the two likelihoods
    //    to implement this algorithm,
    //    generate a uniform random no r in (0,1)  
    //    if log(r) < log((p(D|theta1)/p(D|theta0)) then,
    //      accept the move and theta1 becomes the new theta0
    //    else 
    //      increase the frequency of the old theta 
    if (log(gRandom->Uniform(0.,1.))<(ll1 - ll0)){

      //print the old point data:
      cout<<freq_theta[npoint-1]<<"  mu = "<<mu[npoint-1]<<" sigma = "<<sigma[npoint-1]<<" ll = "<<ll[npoint-1]<<endl;
      //increase the number of success
      npoint++; 

      //update mu0, sigma0 and ll0      
      mu0    = mu1;
      sigma0 = sigma1;
      ll0    = ll1;

      //create a new element in the vectors
      freq_theta.push_back(1.);
      mu.push_back(mu0);
      sigma.push_back(sigma0);
      ll.push_back(ll0);
      
    } else
      //increase the frequency of the old point
      freq_theta[npoint-1] +=1;
  }
  
  for (int ii =0; ii<ll.size(); ii++){
    if (ll[ii]>ll_max) {
      ll_max = ll[ii];
      mu_max = mu[ii];
      sigma_max = sigma[ii];
    }
  }
    cout<<"llmax = "<<ll_max<<" occurs at "<<mu_max<<","<<sigma_max<<endl;
// TH1F * hmu = new TH1F("hmu","Posterior for mu",30, 2.0,8.0);
// TH1F * hsi = new TH1F("hsigma","Posterior for sigma",30, 0.0,2.0);
// TH1F * hll = new TH1F("hl.","Posterior for mu",30, 0.0 ,20.0);
}