# -*- coding: utf-8 -*-
"""
Created on Thu Sep  9 12:02:47 2021

@author: Rene
"""

import numpy as np
import Dist_Library as dl
import matplotlib.pyplot as plt
from scipy.stats import t
from scipy.stats import norm


plt.rc('font', family='serif', size='10')
T_Est_Figure = plt.figure()
T_Est_Figure.set_figwidth(10)
T_Est_Figure.set_figheight(6)
T_Est_Figure.subplots_adjust(hspace=0.4, wspace=0.4)

#********************************************************

# n is the sample size of the dataset, mu and sigma are arbitrary chosen values
n = 2
the_mu = 5
the_sigma = 2

# m is the the number of datasets randomly generated.
m = 100000
t_stat = np.zeros(m)
i_range = np.arange(0,m)
for i in i_range:
    sample = np.random.normal(the_mu,the_sigma,n)
    x_bar = np.mean(sample)
    st_dev_x_bar = np.std(sample,ddof=1)/(n**0.5)
    t_stat[i] = (x_bar-the_mu)/st_dev_x_bar


LB_x = -4
UB_x = 4
off_set_x = 0.5

LB_y = 0
UB_y = 0.45
off_set_y = 0.025

x_values = np.linspace(LB_x, UB_x, 1000) 
y_t_values = t.pdf(x_values,df = n-1)
y_n_values = norm.pdf(x_values,0,1)

x_lims = (LB_x-off_set_x,UB_x+off_set_x)
y_lims = (LB_y-off_set_y,UB_y+off_set_y)

N_bins = 25
bins = np.linspace(LB_x, UB_x, N_bins+1)
T_estimator_hist_data = dl.Estimate_empirical_histogram_table(bins,t_stat)
the_bounds = np.append(T_estimator_hist_data["LB"][0],T_estimator_hist_data["UB"])
the_density = np.append(0,T_estimator_hist_data["Bin PDF"])

Panel = T_Est_Figure.add_subplot(2,2,1)

plt.xlim(x_lims)
plt.ylim(y_lims)
plt.xlabel('Residuals')
plt.ylabel('Density')
plt.axvline(0, lw=2, ls = '-',color = 'lightgray',alpha=0.5)
plt.axhline(0, lw=2, ls = '-',color = 'lightgray',alpha=0.5)
plt.grid(b=None, which='major', axis='both',ls='--',lw=0.5)
dl.add_pmf_density_hist_to_figure_panel(Panel,the_bounds,the_density,'indianred')
plt.plot(x_values,y_n_values,color='green',lw=1,ls='-',label = 'N(0,1) pdf')
plt.plot(x_values,y_t_values,color='blue',lw=1,ls='-',label = 't(n-1) pdf')
plt.legend(loc='best')

text_str = r'$X_i \sim N(\mu,\sigma)$'
plt.text(LB_x,0.4,text_str,color='red', size =9)

text_str = r'$n = $'+str(n)
plt.text(LB_x,0.35,text_str,color='red', size = 9)

text_str = r'$T=\frac{\overline{X}_n-\mu}{s/\sqrt{n}}$'
plt.text(LB_x,0.3,text_str,color='red', size = 9)

# n is the sample size of the dataset, mu and sigma are arbitrary chosen values
n = 4
the_mu = 5
the_sigma = 2

# m is the the number of datasets randomly generated.
m = 100000
t_stat = np.zeros(m)
i_range = np.arange(0,m)
for i in i_range:
    sample = np.random.normal(the_mu,the_sigma,n)
    x_bar = np.mean(sample)
    st_dev_x_bar = np.std(sample,ddof=1)/(n**0.5)
    t_stat[i] = (x_bar-the_mu)/st_dev_x_bar


LB_x = -4
UB_x = 4
off_set_x = 0.5

LB_y = 0
UB_y = 0.45
off_set_y = 0.025

x_values = np.linspace(LB_x, UB_x, 1000) 
y_t_values = t.pdf(x_values,df = n-1)
y_n_values = norm.pdf(x_values,0,1)

x_lims = (LB_x-off_set_x,UB_x+off_set_x)
y_lims = (LB_y-off_set_y,UB_y+off_set_y)

N_bins = 25
bins = np.linspace(LB_x, UB_x, N_bins+1)
T_estimator_hist_data = dl.Estimate_empirical_histogram_table(bins,t_stat)
the_bounds = np.append(T_estimator_hist_data["LB"][0],T_estimator_hist_data["UB"])
the_density = np.append(0,T_estimator_hist_data["Bin PDF"])

Panel = T_Est_Figure.add_subplot(2,2,2)

plt.xlim(x_lims)
plt.ylim(y_lims)
plt.xlabel('Residuals')
plt.ylabel('Density')
plt.axvline(0, lw=2, ls = '-',color = 'lightgray',alpha=0.5)
plt.axhline(0, lw=2, ls = '-',color = 'lightgray',alpha=0.5)
plt.grid(b=None, which='major', axis='both',ls='--',lw=0.5)
dl.add_pmf_density_hist_to_figure_panel(Panel,the_bounds,the_density,'indianred')
plt.plot(x_values,y_n_values,color='green',lw=1,ls='-',label = 'N(0,1) pdf')
plt.plot(x_values,y_t_values,color='blue',lw=1,ls='-',label = 't(n-1) pdf')
plt.legend(loc='best')

text_str = r'$X_i \sim N(\mu,\sigma)$'
plt.text(LB_x,0.4,text_str,color='red', size =9)

text_str = r'$n = $'+str(n)
plt.text(LB_x,0.35,text_str,color='red', size = 9)

text_str = r'$T=\frac{\overline{X}_n-\mu}{s/\sqrt{n}}$'
plt.text(LB_x,0.3,text_str,color='red', size = 9)

# n is the sample size of the dataset, mu and sigma are arbitrary chosen values
n = 16
the_mu = 5
the_sigma = 2

# m is the the number of datasets randomly generated.
m = 100000
t_stat = np.zeros(m)
i_range = np.arange(0,m)
for i in i_range:
    sample = np.random.normal(the_mu,the_sigma,n)
    x_bar = np.mean(sample)
    st_dev_x_bar = np.std(sample,ddof=1)/(n**0.5)
    t_stat[i] = (x_bar-the_mu)/st_dev_x_bar


LB_x = -4
UB_x = 4
off_set_x = 0.5

LB_y = 0
UB_y = 0.45
off_set_y = 0.025

x_values = np.linspace(LB_x, UB_x, 1000) 
y_t_values = t.pdf(x_values,df = n-1)
y_n_values = norm.pdf(x_values,0,1)

x_lims = (LB_x-off_set_x,UB_x+off_set_x)
y_lims = (LB_y-off_set_y,UB_y+off_set_y)

N_bins = 25
bins = np.linspace(LB_x, UB_x, N_bins+1)
T_estimator_hist_data = dl.Estimate_empirical_histogram_table(bins,t_stat)
the_bounds = np.append(T_estimator_hist_data["LB"][0],T_estimator_hist_data["UB"])
the_density = np.append(0,T_estimator_hist_data["Bin PDF"])

Panel = T_Est_Figure.add_subplot(2,2,3)

plt.xlim(x_lims)
plt.ylim(y_lims)
plt.xlabel('Residuals')
plt.ylabel('Density')
plt.axvline(0, lw=2, ls = '-',color = 'lightgray',alpha=0.5)
plt.axhline(0, lw=2, ls = '-',color = 'lightgray',alpha=0.5)
plt.grid(b=None, which='major', axis='both',ls='--',lw=0.5)
dl.add_pmf_density_hist_to_figure_panel(Panel,the_bounds,the_density,'indianred')
plt.plot(x_values,y_n_values,color='green',lw=1,ls='-',label = 'N(0,1) pdf')
plt.plot(x_values,y_t_values,color='blue',lw=1,ls='-',label = 't(n-1) pdf')
plt.legend(loc='best')

text_str = r'$X_i \sim N(\mu,\sigma)$'
plt.text(LB_x,0.4,text_str,color='red', size =9)

text_str = r'$n = $'+str(n)
plt.text(LB_x,0.35,text_str,color='red', size = 9)

text_str = r'$T=\frac{\overline{X}_n-\mu}{s/\sqrt{n}}$'
plt.text(LB_x,0.3,text_str,color='red', size = 9)

# n is the sample size of the dataset, mu and sigma are arbitrary chosen values
n = 32
the_mu = 5
the_sigma = 2

# m is the the number of datasets randomly generated.
m = 100000
t_stat = np.zeros(m)
i_range = np.arange(0,m)
for i in i_range:
    sample = np.random.normal(the_mu,the_sigma,n)
    x_bar = np.mean(sample)
    st_dev_x_bar = np.std(sample,ddof=1)/(n**0.5)
    t_stat[i] = (x_bar-the_mu)/st_dev_x_bar


LB_x = -4
UB_x = 4
off_set_x = 0.5

LB_y = 0
UB_y = 0.45
off_set_y = 0.025

x_values = np.linspace(LB_x, UB_x, 1000) 
y_t_values = t.pdf(x_values,df = n-1)
y_n_values = norm.pdf(x_values,0,1)

x_lims = (LB_x-off_set_x,UB_x+off_set_x)
y_lims = (LB_y-off_set_y,UB_y+off_set_y)

N_bins = 25
bins = np.linspace(LB_x, UB_x, N_bins+1)
T_estimator_hist_data = dl.Estimate_empirical_histogram_table(bins,t_stat)
the_bounds = np.append(T_estimator_hist_data["LB"][0],T_estimator_hist_data["UB"])
the_density = np.append(0,T_estimator_hist_data["Bin PDF"])

Panel = T_Est_Figure.add_subplot(2,2,4)

plt.xlim(x_lims)
plt.ylim(y_lims)
plt.xlabel('Residuals')
plt.ylabel('Density')
plt.axvline(0, lw=2, ls = '-',color = 'lightgray',alpha=0.5)
plt.axhline(0, lw=2, ls = '-',color = 'lightgray',alpha=0.5)
plt.grid(b=None, which='major', axis='both',ls='--',lw=0.5)
dl.add_pmf_density_hist_to_figure_panel(Panel,the_bounds,the_density,'indianred')
plt.plot(x_values,y_n_values,color='green',lw=1,ls='-',label = 'N(0,1) pdf')
plt.plot(x_values,y_t_values,color='blue',lw=1,ls='-',label = 't(n-1) pdf')
plt.legend(loc='best')

text_str = r'$X_i \sim N(\mu,\sigma)$'
plt.text(LB_x,0.4,text_str,color='red', size =9)

text_str = r'$n = $'+str(n)
plt.text(LB_x,0.35,text_str,color='red', size = 9)

text_str = r'$T=\frac{\overline{X}_n-\mu}{s/\sqrt{n}}$'
plt.text(LB_x,0.3,text_str,color='red', size = 9)

# Adding the main title to the plot
Title_str = "T Estimator Histogram and PDF: Sample Size = " + f'{m:,}'
T_Est_Figure.suptitle(Title_str,size='14')

plt.savefig('Student_T_Normal_Sample.png', dpi=300)