###############################################################################
# CARGAMOS LAS LIBRERIAS DE FUNCIONES HECHAS QUE VAMOS A USAR
###############################################################################
import numpy as np #PARA HACER CUENTAS
import xlrd #PARA CARGAR ARCHIVOS DE EXCEL
import matplotlib.pyplot as plt #PARA GRAFICAR
import scipy as sc #PARA ANALISIS DE DATOS
plt.close('all') #CIERRA CUALQUIER GRAFICO QUE HAYA ABIERTO

###############################################################################
# DEFINIMOS NOSOTROS OTRAS DE LAS FUNCIONES QUE VAYAMOS A USAR
###############################################################################
def abrir_columna_excel(nombre_archivo,columna):
  wb = xlrd.open_workbook(nombre_archivo) 
  Hoja1 = wb.sheet_by_index(0) 
  N_datos=len(Hoja1.col_values(columna))
  datos=[]
  for i in range(1,N_datos):
    if type(Hoja1.cell_value(i, columna))==float:
      datos.append(Hoja1.cell_value(i, columna))
  return datos

def residuo(parametros,x,y_medida):
    y_teorica=gaussiana(x,parametros)
    diferencia=y_medida-y_teorica
    return diferencia
  
def gaussiana(x,parametros):
    mu,sigma,A=parametros
    y_teorica=A*np.exp(-(x-mu)**2/(2*sigma**2))
    return y_teorica

###############################################################################
# ESCRIBIMOS LA SECUENCIA DE PROCESOS DE NUESTRO PROGRAMA
###############################################################################
# CARGAMOS LA COLUMNA DE EXEL QUE QUERRAMOS ANALIZAR
nombre_archivo = ("C:\\Users\\Usuario\\Dropbox\\Clase Estadística\\mediciones.xlsx")
columna=2
datos=abrir_columna_excel(nombre_archivo,columna)
N_datos=len(datos)


#N_datos=100000
#datos=np.random.normal(10,1,N_datos)

# HAGO UN HISTOGRAMA CON BINEADO AUTOMATICO
#histo=plt.hist(datos,density=True)
histo=plt.hist(datos)

# DEFINO X e Y PARA AJUSTAR
x=[]
for i in range(len(histo[1])-1):
  x.append((histo[1][i+1]+histo[1][i])/2)
x=np.asarray(x)
y=histo[0]

# AJUSTO
parametros_iniciales=[10,1,1]  
parametros_ajuste=sc.optimize.leastsq(residuo, parametros_iniciales, args=(x,y), Dfun=None, full_output=0, col_deriv=0, ftol=1.49012e-08, xtol=1.49012e-08, gtol=0.0, maxfev=0, epsfcn=None, factor=100, diag=None)

# CALCULO X,Y TEORICOS CON LOS PARAMETROS DE LA GAUSSIANA
x_teorico=np.linspace(min(x),max(x),1000)
y_teorica=gaussiana(x_teorico,parametros_ajuste[0])

mu=parametros_ajuste[0][0]
sigma=parametros_ajuste[0][1]

# GRAFICO EL AJUSTE Y LOS INTERVALOS
plt.plot(x_teorico,y_teorica,label='ajuste')
plt.axvspan(mu-sigma,mu+sigma, alpha=0.25, color='orange',label='1 sigma (68.2 %)')
plt.axvspan(mu-2*sigma,mu+2*sigma, alpha=0.125, color='orange',label='2 sigma (95.4 %)')
plt.axvspan(mu-3*sigma,mu+3*sigma, alpha=0.0625, color='orange',label='3 sigma (99.6 %)')
plt.axvspan(mu-4*sigma,mu+4*sigma, alpha=0.03125, color='orange',label='4 sigma (99.8 %)')
plt.legend()