# -*- coding: utf-8 -*-
#      diese Zeile legt die Codierung von Umlauten fest
##################################################################

# script negLOgLplot.py
''' 
   example of a negative Log Likelihood and its derivatives

.. author:: Günter Quast <g.quast@kit.edu> 
     für den Kurs Computergestützte Datenauswertung (CgDA)
'''
#--------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt


# Beispiele für Funktionsdefinitionen
#   Funktionen aus numpy verwenden, da sie auf array angewendet werden können

# define pdfs
def fgauss(x, mu=0., sigma=1.):
# Gauss distribution 
    return np.exp(-(x-mu)**2/2./sigma**2)/np.sqrt(2.*np.pi)/sigma 

def fexp(x, lam=1.):
# exponential
    return np.exp(-x*lam)*lam

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

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


##### ---- main Program starts here -----

# create figure object
fig = plt.figure('negLogLplot', figsize=(7.5, 10.) )
## numpy-array(s) fuer Stuetzstellen und axis-Objekt(e) erzeugen
ax0 = fig.add_subplot(3,1,1)
ax1 = fig.add_subplot(3,1,2, sharex=ax0)
ax2 = fig.add_subplot(3,1,3, sharex=ax0)

fig.subplots_adjust(hspace=0.15)
plt.setp([a.get_xticklabels() for a in fig.axes[:-1]], visible=False)


# draw random sample from exponential
ns = 300
xm = -np.log(np.random.rand(ns))

def negLogL(tau, xvals):
    nLL=0.
    for x in xvals:
       nLL+= -np.log(fexp(x, tau))
    return nLL

pmin, pmax=0.88, 1.1
npoints=256
par=np.linspace( pmin, pmax, npoints, endpoint=True)

# calculate negLogL and derivatives
nL=negLogL(par, xm)
dnL=np.diff(nL)*(npoints-1.)/(pmin-pmax)
d2nL=np.diff(dnL)*(npoints-2.)/(pmin-pmax)

# get minium and value of negLogL and 2nd derivative
Lmin=min(nL)
idx=nL.tolist().index(Lmin) # minimum of negLogL
phat=par[idx]
curv=min(d2nL)

## plots erzeugen
ax0.plot(par, nL, '-', color='steelblue', linewidth=3)
ax0.plot(par, Lmin+curv/2.*(par-phat)**2, '--', color='darkgreen', linewidth=2,
         alpha=0.5)
ax1.plot((par[:-1]+par[1:])/2., dnL, '-', color='steelblue', linewidth=3)
ax2.plot(par[1:-1], d2nL, '-', color='steelblue', linewidth=3)
ax0.set_xlim(pmin, pmax)
ax1.set_xlim(pmin, pmax)
ax2.set_xlim(pmin, pmax)
ax2.set_ylim(0., 1.5*max(d2nL))



# und Grafiken noch beschriften
ax0.set_title(' ')
ax2.set_xlabel(r'$a$', size='x-large')
ax0.set_ylabel(r'$F(a)=-{\rm ln}({\cal{L}}(a)$)', size='x-large')
ax1.set_ylabel(r'$\partial F(a)/{\partial a}$', size='x-large')
ax2.set_ylabel(r'$\partial^2 F(a)/{\partial a}^2$', size='x-large')

# add some grapical elements

plt.show() # show on screen
#fig.savefig('negLogLplot.pdf')