import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import math



def normal_gaussian(x, mu, h):
    prob = 0
    xmu = (x - mu)
    prob = np.exp(-0.5 * np.power(xmu,2) / (h * h)) / (np.sqrt(2 * math.pi) * h)
    #prob=((xmu<=h/2)*(xmu>=-1*h/2))
    return prob


def calculate_h_optimal(data):
    std = np.sqrt(np.var(data))
    m = np.shape(data)[0]
    n = 1
    dem = (2 * n + 1) * m

    h = std * np.power(4 / dem,1 / (n + 4))
    return h


data = pd.read_csv("temp.txt", sep=',',header=None)
data = data.values


fig = plt.figure()
ax = fig.add_subplot('111')

h_opt = calculate_h_optimal(data)
min_x = np.min(data)
max_x = np.max(data)
number_of_axis_points = np.shape(data)[0]

x = np.linspace(min_x, max_x, number_of_axis_points)

y = np.zeros(number_of_axis_points)

for x_value in data:
    y = y + normal_gaussian(x,x_value,h_opt)

y = y / number_of_axis_points


plt.fill_between(x, y, color='b', alpha=0.2)
plt.plot(x, y, color='b', alpha=0.6)

plt.show()