#Necessary imports

# lib for numeric calculations
import numpy as np

# standard lib for python plotting
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

# seaborn lib for more option in DS
import seaborn as sns

# so to obtain pseudo-random numbers
import random 
                                     
# fits a curve to data
from scipy.optimize import curve_fit                                     

# Lets say that a physics process gives us a gaussian distribution with a
# mean value somewhere between 3 nad 5
# and a sigma with a value around 1
# We do not know exactly their precise number and this is what we want to figure out with the fit.
mu =  random.randrange(3, 5, 1)
sigma = 0.1 * random.randrange(8, 12, 1)

# Create the data
data = np.random.normal(mu, sigma, size=10000)

bin_values, bin_edges = np.histogram(data, density=False, bins=100)
bin_centres = (bin_edges[:-1] + bin_edges[1:])/2

# Define model function to be used to fit to the data above:
def gauss(x, *p):
    A, mu, sigma = p
    return A*np.exp(-(x-mu)**2/(2.*sigma**2))

# p0 is the initial guess for the fitting coefficients (A, mu and sigma above)
p0 = [1., 0., 1.]

coeff, var_matrix = curve_fit(gauss, bin_centres, bin_values, p0=p0)

# Get the fitted curve
hist_fit = gauss(bin_centres, *coeff)

#plt.plot(bin_centres, hist, label='Test data')
plt.plot(bin_centres, hist_fit, label='Fitted data', color='red', linewidth=2)
_ = plt.hist(data, bins=100, color='blue', alpha=.3)

# Finally, lets get the fitting parameters, i.e. the mean and standard deviation:
print 'Fitted mean = ', coeff[1]
print 'Fitted standard deviation = ', coeff[2]

plt.show()

"""
Explanation: The typical problem we have to solve:
Lets say we want to study the properties of a distribution e.g. physics process.
Typically we see the histogram and we want to see to which class of distribution belongs to.
we can obtain this by fitting our data and checking the residuals (goodness of fit)
End of explanation
"""


# Lets say that a physics process gives us a gaussian distribution with a
# mean value somewhere between 3 nad 5
# while it has a broad sigma with a value around 1
mu =  random.randrange(3, 5, 1)
sigma = 0.1 * random.randrange(8, 12, 1)

# Lets create now the distribution
data = np.random.normal(mu, sigma, 1000)
 
# Lets visualize it
_ = plt.hist(data, bins=100, alpha=.5)

"""
Explanation: Create the distribution and visualize it
End of explanation
"""


# get the binned data
# density = True => histo integral = 1, normalized
hist_data = np.histogram(data, density=True, bins=100) 

# Define model function to be used to fit to the data above, we assume here that we know that is a gaussian function:
def gauss(x, *p):
    A, mu, sigma = p
    return A * np.exp( - (x - mu) ** 2 / (2. * sigma ** 2))


# get distribution information
bin_edges = hist_data[1]
bin_centres = (bin_edges[:-1] + bin_edges[1:]) / 2
bin_values = hist_data[0]

# p0 is the initial guess for the fitting coefficients (A, mu and sigma above)
p0 = [1., 4., 1.]

# Fit the data
coeff, var_matrix = curve_fit(gauss, bin_centres, bin_values, p0=p0)
print "The fit coeff. are:", coeff

# Obtain the final fit function
hist_fit = gauss(bin_centres, *coeff)

hist_data = plt.hist(data, bins=100, alpha = 0.5)
plt.plot(bin_centres, hist_fit, label='Fitted data')



#hist_data = plt.hist(s, bins=100)

"""
Explanation: Fit a function to the distribution and obtain its properties
End of explanation
"""





# Get the fitted curve
hist_fit = gauss(bin_centres, *coeff)



hist_data = plt.hist(s, bins=100, alpha = 0.5)
plt.plot(bin_centres, hist_fit, label='Fitted data')

coeff

"""
Explanation: Now do the fit
End of explanation
"""


