# -*- coding: utf-8 -*-
"""
Created on Wed Sep 15 11:32:26 2021

@author: Rene
"""

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

# implementing two-tailed t-test for the popultation mean

def one_sided_z_test(x, mu_0, the_sigma, alpha):
    n = len(x)
    x_bar = np.mean(x)
    z_quantile = norm.ppf(1-alpha,0,1)
    st_dev_x_bar = the_sigma/(n**0.5)
    z_statistic = (x_bar-mu_0)/st_dev_x_bar
    p_value = 1-norm.cdf(abs(z_statistic),0,1)
    results = np.zeros(3)
    results[0] = z_statistic
    results[1] = z_quantile
    results[2] = p_value
    return(results)

def TypeII_error_one_sided_z_test(x,alpha,mu_0,mu_1,the_sigma):
  n = len(x)
  st_dev_x_bar = the_sigma/(n**0.5)
  z_quantile = norm.ppf(1-alpha,mu_0,st_dev_x_bar)
  TypeII_error=norm.cdf(z_quantile,mu_1,st_dev_x_bar)
  return(TypeII_error)

# Evaluating x_values and y_values to fill for a tail

def evaluate_z_fill_coordinates(LB,UB,m,mu,stdev):
  x = np.linspace(LB,UB,m)
  x_trans = (x-mu)/stdev
  y = norm.pdf(x_trans,0,1)/stdev
  return(x,y)


data = np.array([123,122,122.5])
n = len(data)
x_bar = np.mean(data)

# H_0: mu = mu_0  H_1: mu not = mu_0

mu_0 = 120
alpha = 0.05
st_dev_known = 2
st_dev_x_bar = st_dev_known/(n**0.5)

# Evaluating the Type II Error for a specific-value of mu_1

mu_1 = 124
beta = TypeII_error_one_sided_z_test(data,alpha,mu_0,mu_1,st_dev_known)
power = 1-beta
The_File_Name = 'Dutch_Speed_example_Type_II_Error_5.png'

# Assuming H_0 is True and with estimate for sigma evaluating theoretical pdf

LB = 110
UB = 130
x_points = np.linspace(110,130,1000)
pdf_points_mu_0 = norm.pdf(x_points,mu_0,st_dev_x_bar)
pdf_points_mu_1 = norm.pdf(x_points,mu_1,st_dev_x_bar)

off_set_x = 0.5
off_set_y = 0.02
max_y=0.425

# Evaluating the pdf at the critical value if population mean equals mu_0
crit_value = norm.ppf(1-alpha,mu_0,st_dev_x_bar)
pdf_crit_value = norm.pdf(crit_value,mu_0,st_dev_x_bar)

# Evaluating the pdf at the critical value if population mean equals mu_1
pdf_mu_1_UB = norm.pdf(crit_value,mu_1,st_dev_x_bar)

# Setting the attributes for the Figure     
plt.rc('font', family='serif', size='10')
Figure = plt.figure()
Figure.set_figwidth(9)
Figure.set_figheight(5)
Figure.subplots_adjust(hspace=0.4, wspace=0.4)
x_lims = (LB-off_set_x,UB+off_set_x)
y_lims = (0-off_set_y,max_y+off_set_y)
plt.axvline(0, lw=2, ls = '-',color = 'lightgray',alpha=0.5)
plt.axhline(0, lw=2, ls = '-',color = 'lightgray',alpha=0.5)
plt.xlim(x_lims)
plt.ylim(y_lims)
plt.xlabel(r'$\bar{x}$')
plt.ylabel(r'pdf $\bar{X}$')

plt.axvline(mu_0, lw=1, ls = '--',color = 'blue',alpha=0.5)
plt.axvline(mu_1, lw=1, ls = '--',color = 'darkgreen',alpha=0.5)
plt.axvline(crit_value, lw=1, ls = '--',color = 'red',alpha=0.5)

plt.plot(x_points, pdf_points_mu_0, color = 'blue', linestyle="-",lw=1)
plt.plot(x_points, pdf_points_mu_1, color = 'darkgreen', linestyle="-",lw=1)

t_x_values = [crit_value, crit_value]
t_y_values = [0, pdf_crit_value]
plt.plot(t_x_values, t_y_values, color = 'red', linestyle="--",lw=1,alpha = 0.5)

t_x_values, t_y_values = evaluate_z_fill_coordinates(crit_value,UB,100,mu_0,st_dev_x_bar)
plt.fill_between(t_x_values, np.zeros(100), t_y_values, color = 'indianred',alpha=0.25)

t_x_values, t_y_values = evaluate_z_fill_coordinates(110,crit_value,100,mu_1,st_dev_x_bar)
plt.fill_between(t_x_values, np.zeros(100), t_y_values, color = 'darkgreen',alpha=0.25)

text_str = r'$\mu_0 = $' + f'{mu_0:2.1f}'
plt.text(mu_0-0.2,max_y-0.025,text_str,horizontalalignment='right',color = 'blue')
   
text_str = r'$c_{\alpha} = $' + f'{crit_value:2.1f}'
plt.text(crit_value+0.2,max_y-0.025,text_str,horizontalalignment='left',color = 'red')

      
text_str = r'$\mu_1 = $' + f'{mu_1:2.1f}'
plt.text(mu_1+0.2,max_y-0.05,text_str,horizontalalignment='left',color = 'green')

off_set_x = 0.15
text_str = r'$H_0: \mu = \mu_0$'
plt.text(LB-off_set_x+0.25,max_y-0.05,text_str,horizontalalignment='left',color = 'black')

text_str = r'$H_1: \mu > \mu_0$'
plt.text(LB-off_set_x+0.25,0.925*max_y-0.05,text_str,horizontalalignment='left',color = 'black')

text_str = r'$Z\,=\,(\overline{X}-\mu_0)/(\sigma/\sqrt{n})$'
plt.text(LB-off_set_x+0.25,0.85*max_y-0.05,text_str,horizontalalignment='left',color = 'red')

text_str = r'$Z\sim N(0,1)$'
plt.text(LB-off_set_x+0.25,0.775*max_y-0.05,text_str,horizontalalignment='left',color = 'red')

text_str = r'$\overline{X}(\mu_0) = (\sigma/\sqrt{n}) \times Z + \mu_0$'
plt.text(LB-off_set_x+0.25,0.7*max_y-0.05,text_str,horizontalalignment='left',color = 'blue')

text_str = r'$\overline{X}(\mu_1) = (\sigma/\sqrt{n}) \times Z + \mu_1$'
plt.text(LB-off_set_x+0.25,0.625*max_y-0.05,text_str,horizontalalignment='left',color = 'darkgreen')

text_str = r'$Type\,I\,Error = Pr (\overline{X}(\mu_0) > c_{\alpha}) $'
plt.text(LB-off_set_x+0.25,0.475*max_y-0.05,text_str,horizontalalignment='left',color = 'indianred')

text_str = r'$Type\,I\,Error = \alpha = $' + f'{100*alpha:2.0f}'+'%'
plt.text(LB-off_set_x+0.25,0.4*max_y-0.05,text_str,horizontalalignment='left',color = 'indianred')

#text_str = r'$Type\,II\,Error = \hat{\beta}(\mu_1) = $' + f'{100*beta:2.2f}'+'%'
#plt.text(LB-off_set_x+0.25,0.325*max_y-0.05,text_str,horizontalalignment='left',color = 'darkgreen')

text_str = r'$Type\,II\,Error = Pr (\overline{X}(\mu_1) < c_{\alpha})$'
plt.text(LB-off_set_x+0.25,0.325*max_y-0.05,text_str,horizontalalignment='left',color = 'darkgreen')

text_str = r'$Type\,II\,Error = \beta(\mu_1) = $' + f'{100*beta:2.2f}'+'%'
plt.text(LB-off_set_x+0.25,0.25*max_y-0.05,text_str,horizontalalignment='left',color = 'darkgreen')

text_str = r'Dutch Speed Example: $T\sim N(\mu_0, \sigma/\sqrt{n}):\mu_0=120,\,n=3,\,\alpha=$ 'f'{100*alpha:2.0f}'+'%'
Figure.suptitle(text_str,size='14')

plt.savefig(The_File_Name, dpi=300)