# coding: utf-8

# In[1]:

get_ipython().magic('matplotlib inline')
import numpy as np
import matplotlib.pyplot as plt


# In[2]:

hd = np.loadtxt("hd188041.dat", usecols=(7,), unpack=True)
hr = np.loadtxt("hr7575.dat", usecols=(7,), unpack=True)


# In[3]:

def histogram(data, bins=10):
    m = np.min(data)
    M = np.max(data)
    w = (M-m)/bins
    x = np.linspace(m, M, bins+1)
    N = np.zeros(bins)
    for d in data:
        mn = d <= x
        n = np.argmax(mn)
        if n > 0: n -= 1 
        N[n] += 1
    plt.bar(x[:-1], N, width=w)


# In[4]:

#vykresleni klasickeho histogramu pro srovnani
plt.figure(figsize=(16,4))
plt.subplot(121), plt.title("hd188041.dat"), histogram(hd)
plt.subplot(122), plt.title("hr7575.dat"), histogram(hr);


# In[5]:

plt.figure(figsize=(16,4))
plt.subplot(121), plt.title("hd188041.dat"), plt.hist(hd);
plt.subplot(122), plt.title("hr7575.dat"), plt.hist(hr);


# In[6]:

def mad(x):
    """Median absolut deviation"""
    return np.median(np.abs(x-np.median(x)))


# In[7]:

def Lkurt(x):
    n = len(x)
    return 0.22 * np.sqrt(n) * ((n+1) * np.mean((x-np.mean(x))**4)
                             /  (n-1) / np.std(x)**4 - 3)


# In[8]:

def Lmad(x):
    n = len(x)
    return np.sqrt(n) * (1.084 - (1.608 * mad(x)) / np.std(x))


# In[9]:

Lkurt(hd), Lmad(hd)


# In[10]:

Lkurt(hr), Lmad(hr)


# In[11]:

Norm = np.random.normal(scale=0.01, size=1000)


# In[12]:

histogram(Norm)


# In[13]:

Lkurt(Norm), Lmad(Norm)


# In[ ]:

#promenna hvezda je hd188041