#! /usr/bin/env python3
"""test_mlFit.py
   Maximum likelihood fit to unbinned data with package phyFit and iminiut

.. moduleauthor:: Guenter Quast <g.quast@kit.edu>

"""

import numpy as np, matplotlib.pyplot as plt
# from PhyPraKit.phyFit import mFit
from PhyPraKit import mFit

## some useful distributions
#Binomial distribution
def fBinomial(x,n,p):
  k=np.around(x)
  return sp.binom(n,k) * p**k * (1.-p)**(n-k)

# Poisson distribution
def fPoisson(x,mu):
  k=np.around(x)
  return (mu**k)/np.exp(mu)/sp.gamma(k+1.)

# Gauss distribution
def fGauss(x, mu=10., sigma=0.5):
    return (np.exp(-(x-mu)**2/2./sigma**2)/np.sqrt(2.*np.pi)/sigma)


if __name__ == "__main__": # --------------------------------------  
  #
  # **unbinned maximum-likelihodd fit** of pdf to unbinned data
  #
  
  #  real data from measurement with a Water Cherenkov detector ("Kamiokanne")
  #
  #    numbers represent time differences (in µs) between the passage of a muon
  #    and the registration of a second pulse, often caused by an electron from
  #    the decay of the stopped muon. As such events are rare, histogramming the
  #    data prior to fitting would introduce shifts and biases, and therefore the

  # set Model to use
  modelPDF=fGauss

  # genereta data (here: drwan from Gaussian distribution)
  mu0=2.
  sig0=0.5
  np.random.seed(314159)  # initialize random generator
  Data = mu0 + sig0 * np.random.randn(100)
  
  rdict = mFit( modelPDF,
          data = Data,   # data - if not None, normalised PDF is assumed as model  
          p0=None,                 # initial guess for parameter values 
        #  constraints=[['tau', 2.2, 0.01], # Gaussian parameter constraints
          neg2logL = True,         # use  -2 * ln(L)
          plot=True,               # plot data and model
          plot_band=True,          # plot model confidence-band
          plot_cor=True,          # plot profiles likelihood and contours
          showplots=False,         # show / don't show plots
          quiet=False,             # suppress informative printout if True
          axis_labels=['data x', 
                       'Probability Density  pdf(x; *p)'], 
          data_legend = 'data',    
          model_legend = 'Gauß' )
          
  plt.suptitle("Unbinned ML to Gaussian data",
                size='xx-large', color='darkblue')

  # Print results
#  pvals, perrs, cor, gof, pnams = rdict.values()
#  print('\n*==* unbinned ML Fit Result:')
#  print(" parameter names:       ", pnams)
#  print(" parameter values:      ", pvals)
#  print(" neg. parameter errors: ", perrs[:,0])
#  print(" pos. parameter errors: ", perrs[:,1])
#  print(" correlation matrix : \n", cor)  
## new version
  print('\n*==* unbinned ML Fit Result:')
  for key in rdict:
    print("{}\n".format(key), rdict[key])
  
  plt.show()