import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy import stats

%matplotlib notebook

df = pd.read_csv('http://ws-cheetah.stm.fce.vutbr.cz/vyuka/CD004/priklady/tlakova_zkouska/data.txt', 
                 sep=';', header=None)
df

f_cube_data = df.values.flatten()
f_cube_data

f_cube = []
for i in range(100):
    f = pd.read_csv('http://ws-cheetah.stm.fce.vutbr.cz/vyuka/CD004/priklady/tlakova_zkouska/ld_curves/krivka_{:02d}.txt'.format(i),
                    sep=';', header=None, usecols=(1,), skiprows=2)
    f_cube.append(f.values.max())
f_cube = np.array(f_cube)
f_cube

n = len(f_cube)
n

min_ = f_cube.min()
min_

max_ = f_cube.max()
max_

mu = f_cube.mean()
mu

std = f_cube.std()
std

std = f_cube.std(ddof=1)
std

skew = stats.skew(f_cube)
skew

skew = stats.skew(f_cube, bias=False)
skew

kurt = stats.kurtosis(f_cube)
kurt

CoV = std / mu
CoV

params = stats.norm.fit(f_cube)
rv_norm = stats.norm(*params)
print('params =', params)
print('mu, std =', rv_norm.mean(), rv_norm.std())

params = stats.lognorm.fit(f_cube, floc=0)
rv_lognorm = stats.lognorm(*params)
print('params =', params)
print('mu, std =', rv_lognorm.mean(), rv_lognorm.std())

params = stats.weibull_min.fit(f_cube, floc=0)
rv_weib = stats.weibull_min(*params)
print('params =', params)
print('mu, std =', rv_weib.mean(), rv_weib.std())

fig, ax = plt.subplots(figsize=(6,3), tight_layout=True)
ax.hist(f_cube);

fig, ax = plt.subplots(figsize=(6,3), tight_layout=True)
ax.hist(f_cube, bins='sqrt', edgecolor='k', normed=True, label='data', alpha=.5);

fig, ax = plt.subplots(figsize=(6,3), tight_layout=True)
ax.hist(f_cube, bins='sqrt', ec='k', normed=True, label='data', alpha=.5)

x = np.linspace(20, 45, 100)
ax.plot(x, rv_norm.pdf(x), label='norm')

ax.plot(x, rv_lognorm.pdf(x), label='lognorm')

ax.plot(x, rv_weib.pdf(x), label='weibull')

ax.legend();

fig, ax = plt.subplots(figsize=(6,3), tight_layout=True)
ax.hist(f_cube, bins='sqrt', ec='k', normed=True, label='data', 
        alpha=.5, cumulative=True)

x = np.linspace(20, 45, 100)
ax.plot(x, rv_norm.cdf(x), label='norm')

ax.plot(x, rv_lognorm.cdf(x), label='lognorm')

ax.plot(x, rv_weib.cdf(x), label='weibull')

ax.legend();

def print_quantils(rv, name):
    print(name)
    print('\tf_0.05  =', rv.ppf(0.05))
    print('\tf_0.95  =', rv.ppf(0.95))
    print('\tf_0.001 =', rv.ppf(0.001))

print_quantils(rv_norm, 'rv_norm')

print_quantils(rv_lognorm, 'rv_lognorm')

print_quantils(rv_weib, 'rv_weib')