import numpy as np

#generate random examples for computing the confusion matrix
true_mean=0.5
true_scale=0.5
y_true = np.random.normal(loc=true_mean, scale=true_scale, size=100)
print("non-bootstrapped point estimate mean", np.mean(y_true), "standard error", np.std(y_true)/np.sqrt(len(y_true)))

def bootstrap(data, statistic, boostrap_size, num_boostrap_samples=1000):
    indices=range(len(data))
    statistic_values=[]
    for i in range(num_boostrap_samples):
        choice= np.random.choice(indices, size = boostrap_size, replace = True)
        y_true = data[choice]+np.random.normal(loc=0, scale = 0.3, size=boostrap_size)
        statistic_value = statistic(y_true)
        statistic_values.append(statistic_value)
    return np.average(statistic_values, axis = 0), np.std(statistic_values, axis = 0), statistic_values
mean_statistic, std_statistic, statistic_values= bootstrap(y_true, np.mean, boostrap_size= int(len(y_true)))

print("bootstrapped estimate for the mean", mean_statistic)
print("bootsrapped estimate for the standard error of mean", std_statistic)

a=true_mean
b=true_scale/np.sqrt(len(y_true)) #np.std(y_true, ddof=1)/np.sqrt(len(y_true))

import matplotlib.pyplot as plt
plt.hist(statistic_values, density=True, bins=50, label = "(Smooth) bootstrap density of mean values")
xs=np.linspace(0.2,0.8, num = 100)
#plt.plot(xs, 1/(np.sqrt(2*np.pi)*std_statistic)*np.exp(-(xs-mean_statistic)**2/(2*std_statistic**2)), label = "(Smooth) bootstrap distribution of mean values")
plt.vlines(np.mean(y_true), 0,9, color = "red", label = "Estimate mean value")
plt.plot(xs, 1/(np.sqrt(2*np.pi)*b)*np.exp(-(xs-a)**2/(2*b**2)), label = "Ground truth normal distribution of sample mean values")
plt.xlabel("Mean value")
plt.ylabel("Proability density")
plt.legend()
plt.show()