plot_material.py

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

###############################################################################################################    
###############################################################################################################
def reg_corr_plot():    
    class LinearRegression: 
        def __init__(self, beta1, beta2, error_scale, data_size):
            self.beta1 = beta1
            self.beta2 = beta2
            self.error_scale = error_scale
            self.x = np.random.randint(1, data_size, data_size)
            self.y = self.beta1 + self.beta2*self.x + self.error_scale*np.random.randn(data_size)

        def x_y_cor(self):
            return np.corrcoef(self.x, self.y)[0, 1]

    fig, ax = plt.subplots(nrows = 2, ncols = 4,figsize=(24, 12))

    beta1, beta2, error_scale,data_size =  2, .05, 1, 100
    lrg1 = LinearRegression(beta1, beta2, error_scale, data_size)
    ax[0, 0].scatter(lrg1.x, lrg1.y)
    ax[0, 0].plot(lrg1.x, beta1 + beta2*lrg1.x, color = '#FA954D', alpha = .7)
    ax[0, 0].set_title(r'$Y={}+{}X+{}u$'.format(beta1, beta2, error_scale))
    ax[0, 0].annotate(r'$\rho={:.4}$'.format(lrg1.x_y_cor()), xy=(0.1, 0.9), xycoords='axes fraction')

    beta1, beta2, error_scale,data_size =  2, -.6, 1, 100
    lrg2 = LinearRegression(beta1, beta2, error_scale, data_size)
    ax[0, 1].scatter(lrg2.x, lrg2.y)
    ax[0, 1].plot(lrg2.x, 2 - .6*lrg2.x, color = '#FA954D', alpha = .7)
    ax[0, 1].set_title(r'$Y={}+{}X+{}u$'.format(beta1, beta2, error_scale))
    ax[0, 1].annotate(r'$\rho={:.4}$'.format(lrg2.x_y_cor()), xy=(0.1, 0.9), xycoords='axes fraction')

    beta1, beta2, error_scale,data_size =  2, 1, 1, 100
    lrg3 = LinearRegression(beta1, beta2, error_scale, data_size)
    ax[0, 2].scatter(lrg3.x, lrg3.y)
    ax[0, 2].plot(lrg3.x, beta1 + beta2 * lrg3.x, color = '#FA954D', alpha = .7)
    ax[0, 2].set_title(r'$Y={}+{}X+{}u$'.format(beta1, beta2, error_scale))
    ax[0, 2].annotate(r'$\rho={:.4}$'.format(lrg3.x_y_cor()), xy=(0.1, 0.9), xycoords='axes fraction')

    beta1, beta2, error_scale,data_size =  2, 3, 1, 100
    lrg4 = LinearRegression(beta1, beta2, error_scale, data_size)
    ax[0, 3].scatter(lrg4.x, lrg4.y)
    ax[0, 3].plot(lrg4.x, beta1 + beta2 * lrg4.x, color = '#FA954D', alpha = .7)
    ax[0, 3].set_title(r'$Y={}+{}X+{}u$'.format(beta1, beta2, error_scale))
    ax[0, 3].annotate(r'$\rho={:.4}$'.format(lrg4.x_y_cor()), xy=(0.1, 0.9), xycoords='axes fraction')

    beta1, beta2, error_scale,data_size =  2, 3, 3, 100
    lrg5 = LinearRegression(beta1, beta2, error_scale, data_size)
    ax[1, 0].scatter(lrg5.x, lrg5.y)
    ax[1, 0].plot(lrg5.x, beta1 + beta2 * lrg5.x, color = '#FA954D', alpha = .7)
    ax[1, 0].set_title(r'$Y={}+{}X+{}u$'.format(beta1, beta2, error_scale))
    ax[1, 0].annotate(r'$\rho={:.4}$'.format(lrg5.x_y_cor()), xy=(0.1, 0.9), xycoords='axes fraction')

    beta1, beta2, error_scale,data_size =  2, 3, 10, 100
    lrg6 = LinearRegression(beta1, beta2, error_scale, data_size)
    ax[1, 1].scatter(lrg6.x, lrg6.y)
    ax[1, 1].plot(lrg6.x, beta1 + beta2 * lrg6.x, color = '#FA954D', alpha = .7)
    ax[1, 1].set_title(r'$Y={}+{}X+{}u$'.format(beta1, beta2, error_scale))
    ax[1, 1].annotate(r'$\rho={:.4}$'.format(lrg6.x_y_cor()), xy=(0.1, 0.9), xycoords='axes fraction')

    beta1, beta2, error_scale,data_size =  2, 3, 20, 100
    lrg7 = LinearRegression(beta1, beta2, error_scale, data_size)
    ax[1, 2].scatter(lrg7.x, lrg7.y)
    ax[1, 2].plot(lrg7.x, beta1 + beta2 * lrg7.x, color = '#FA954D', alpha = .7)
    ax[1, 2].set_title(r'$Y={}+{}X+{}u$'.format(beta1, beta2, error_scale))
    ax[1, 2].annotate(r'$\rho={:.4}$'.format(lrg7.x_y_cor()), xy=(0.1, 0.9), xycoords='axes fraction')


    beta1, beta2, error_scale,data_size =  2, 3, 50, 100
    lrg8 = LinearRegression(beta1, beta2, error_scale, data_size)
    ax[1, 3].scatter(lrg8.x, lrg8.y)
    ax[1, 3].plot(lrg3.x, beta1 + beta2 * lrg3.x, color = '#FA954D', alpha = .7)
    ax[1, 3].set_title(r'$Y={}+{}X+{}u$'.format(beta1, beta2, error_scale))
    ax[1, 3].annotate(r'$\rho={:.4}$'.format(lrg8.x_y_cor()), xy=(0.1, 0.9), xycoords='axes fraction')
    
###############################################################################################################    
###############################################################################################################

def central_limit_theorem_plot():    
    fig, ax = plt.subplots(4, 3, figsize = (20, 20))

    ########################################################################################
    x = np.linspace(2, 8, 100)
    a = 2 # range of uniform distribution
    b = 8
    unif_pdf = np.ones(len(x)) * 1/(b-a)

    ax[0, 0].plot(x, unif_pdf, lw = 3, color = 'r')
    ax[0, 0].plot([x[0],x[0]],[0, 1/(b-a)], lw = 3, color = 'r', alpha = .9) # vertical line
    ax[0, 0].plot([x[-1],x[-1]],[0, 1/(b-a)], lw = 3, color = 'r', alpha = .9)
    ax[0, 0].fill_between(x, 1/(b-a), 0, alpha = .5, color = 'r')

    ax[0, 0].set_xlim([1, 9])
    ax[0, 0].set_ylim([0, .4])
    ax[0, 0].set_title('Uniform Distribution', size = 18)
    ax[0, 0].set_ylabel('Population Distribution', size = 12)

    ########################################################################################
    ss = 2 #sample size
    unif_sample_mean = np.zeros(1000)
    for i in range(1000):
        unif_sample = np.random.rand(ss)
        unif_sample_mean[i] = np.mean(unif_sample)
    ax[1, 0].hist(unif_sample_mean, bins = 20, color = 'r', alpha = .5)
    ax[1, 0].set_ylabel('Sample Distribution， $n = 2$', size = 12)
    ########################################################################################
    ss = 10 #sample size
    unif_sample_mean = np.zeros(1000)
    for i in range(1000):
        unif_sample = np.random.rand(ss)
        unif_sample_mean[i] = np.mean(unif_sample)
    ax[2, 0].hist(unif_sample_mean, bins = 30, color = 'r', alpha = .5)
    ax[2, 0].set_ylabel('Sample Distribution, $n = 10$', size = 12)

    ########################################################################################
    ss = 1000 #sample size
    unif_sample_mean = np.zeros(1000)
    for i in range(1000):
        unif_sample = np.random.rand(ss)
        unif_sample_mean[i] = np.mean(unif_sample)
    ax[3, 0].hist(unif_sample_mean, bins = 40, color = 'r', alpha = .5)
    ax[3, 0].set_ylabel('Sample Distribution, $n = 1000$', size = 12)

    ########################################################################################
    a = 6
    b = 2
    x = np.linspace(0, 1, 100)
    beta_pdf = sp.stats.beta.pdf(x, a, b)
    ax[0, 1].plot(x, beta_pdf, lw = 3, color = 'g')
    ax[0, 1].set_ylim([0, 6])
    ax[0, 1].fill_between(x, beta_pdf, 0, alpha = .5, color = 'g')
    ax[0, 1].set_title('Beta Distribution', size = 18)
    ########################################################################################

    ss = 2 #sample size
    beta_sample_mean = np.zeros(1000)
    for i in range(1000):
        beta_sample = sp.stats.beta.rvs(a, b, size = ss)
        beta_sample_mean[i] = np.mean(beta_sample)
    ax[1, 1].hist(beta_sample_mean, color = 'g', alpha = .5)

    ########################################################################################

    ss = 10 #sample size
    beta_sample_mean = np.zeros(1000)
    for i in range(1000):
        beta_sample = sp.stats.beta.rvs(a, b, size = ss)
        beta_sample_mean[i] = np.mean(beta_sample)
    ax[2, 1].hist(beta_sample_mean, color = 'g', bins = 20, alpha = .5)


    ########################################################################################

    ss = 100000 #sample size
    beta_sample_mean = np.zeros(1000)
    for i in range(1000):
        beta_sample = sp.stats.beta.rvs(a, b, size = ss)
        beta_sample_mean[i] = np.mean(beta_sample)
    ax[3, 1].hist(beta_sample_mean, color = 'g', bins = 30, alpha = .5)

    ########################################################################################
    a = 6
    x = np.linspace(0, 25, 100)

    gamma_pdf = sp.stats.gamma.pdf(x, a)
    ax[0, 2].plot(x, gamma_pdf, lw = 3, color = 'b')
    ax[0, 2].set_ylim([0, 0.34])
    ax[0, 2].fill_between(x, gamma_pdf, 0, alpha = .5, color = 'b')
    ax[0, 2].set_title('Gamma Distribution', size = 18)

    ########################################################################################
    ss = 2 #sample size
    gamma_sample_mean = np.zeros(1000)
    for i in range(1000):
        gamma_sample = sp.stats.gamma.rvs(a, size = ss)
        gamma_sample_mean[i] = np.mean(gamma_sample)
    ax[1, 2].hist(gamma_sample_mean, color = 'b', alpha = .5)

    ########################################################################################
    ss = 10 #sample size
    gamma_sample_mean = np.zeros(1000)
    for i in range(1000):
        gamma_sample = sp.stats.gamma.rvs(a, size = ss)
        gamma_sample_mean[i] = np.mean(gamma_sample)
    ax[2, 2].hist(gamma_sample_mean, bins = 20, color = 'b', alpha = .5)
    ########################################################################################
    ss = 1000 #sample size
    gamma_sample_mean = np.zeros(1000)
    for i in range(1000):
        gamma_sample = sp.stats.gamma.rvs(a, size = ss)
        gamma_sample_mean[i] = np.mean(gamma_sample)
    ax[3, 2].hist(gamma_sample_mean, bins = 30, color = 'b', alpha = .5)
    ########################################################################################
    plt.show()
    
##########################################################################################################
##########################################################################################################
def type12_error():
    x = np.linspace(-6, 9, 200)
    null_loc, alter_loc = 0, 3
    y_null = sp.stats.norm.pdf(x, loc = null_loc)
    y_alter =  sp.stats.norm.pdf(x, loc = alter_loc)
    fig, ax = plt.subplots(figsize = (18, 6))
    ax.plot(x, y_null, x, y_alter)
    ax.annotate('Null', (null_loc-.2, max(y_null)/2), size = 15)
    ax.annotate('Alternative', (alter_loc-.6, max(y_alter)/2), size = 15)
    ax.annotate('Type I Error', (2, max(y_alter)/30), size = 15)
    ax.annotate('Type II Error', (0, max(y_alter)/30), size = 15)
    ax.fill_between(x[-98:], y_null[-98:])
    ax.fill_between(x[:103], y_alter[:103])
    ax.set_ylim([0, .5])
    plt.show()
    
##########################################################################################################
##########################################################################################################   
# def reject_region(): 
#     data = pd.read_csv('500_Person_Gender_Height_Weight_Index.csv')
    
#     male_mean = data[data['Gender']=='Male']['Height'].mean()
#     male_std = data[data['Gender']=='Male']['Height'].std(ddof=1)
#     male_std_error = male_std/np.sqrt(len(data[data['Gender']=='Male']))
#     male_null = 172
    
#     df = len(data[data['Gender']=='Male'])-1
#     t_975 = sp.stats.t.ppf(.975, df=df)
#     t_025 = sp.stats.t.ppf(.025, df=df)

#     x = np.linspace(male_null-5, male_null+5, 200)
#     df = len(data[data['Gender']=='Male'])-1
#     y_t = sp.stats.t.pdf(x, df = df, loc = male_null)

#     fig, ax = plt.subplots(2, 1, figsize = (18,8))

#     ax[0].plot(x, y_t, color = 'tomato', lw = 3)

#     rejection_lower = male_null - t_975*male_std_error
#     x_rej_lower = np.linspace(rejection_lower-3, rejection_lower, 30)
#     y_rej_lower = sp.stats.t.pdf(x_rej_lower, df = df, loc = male_null)
#     ax[0].fill_between(x_rej_lower, y_rej_lower, color = 'tomato', alpha = .7)

#     rejection_upper = male_null + t_975*male_std_error
#     x_rej_upper = np.linspace(rejection_upper, rejection_upper+3, 30)
#     y_rej_upper = sp.stats.t.pdf(x_rej_upper, df = df, loc = male_null)
#     ax[0].fill_between(x_rej_upper, y_rej_upper, color = 'tomato', alpha = .7)

#     ax[0].set_ylim([0, .45])

#     x = np.linspace(-5, 5, 200)
#     y_t = sp.stats.t.pdf(x, df = df, loc = 0)

#     ax[1].plot(x, y_t, color = 'tomato', lw = 3)

#     x_rej_lower = np.linspace(t_025-3, t_025, 30)
#     y_rej_lower = sp.stats.t.pdf(x_rej_lower, df = df)
#     ax[1].fill_between(x_rej_lower, y_rej_lower, color = 'tomato', alpha = .7)

#     x_rej_lower = np.linspace(t_975+3, t_975, 30)
#     y_rej_lower = sp.stats.t.pdf(x_rej_lower, df = df)
#     ax[1].fill_between(x_rej_lower, y_rej_lower, color = 'tomato', alpha = .7)

#     ax[1].set_ylim([0, .45])

#     plt.show()
    
##########################################################################################################
# ##########################################################################################################
# def draw_something():
#     x = np.linspace(0, 10, 100)
#     y = np.sin(x)
#     plt.plot(x, y)

##########################################################################################################
##########################################################################################################
def anova_plot():
    def gen_3samples(loc1, loc2, loc3, scale1, scale2, scale3, size1, size2, size3):
        F_statistic, p_value = [], []
        for i in range(1000):
            a = sp.stats.norm.rvs(loc1, scale1, size1)
            b = sp.stats.norm.rvs(loc2, scale2, size3)
            c = sp.stats.norm.rvs(loc3, scale3, size3)
            F, p = sp.stats.f_oneway(a,b,c)
            F_statistic.append(F)
            p_value.append(p)
        return F_statistic, p_value

    fig, ax = plt.subplots(nrows = 6, ncols = 3, figsize = (17, 34))

    mu1, mu2, mu3, sig1, sig2, sig3, size1, size2, size3 = [3, 6, 9, 6, 6, 6, 10, 20, 30]
    params = [mu1, mu2, mu3, sig1, sig2, sig3, size1, size2, size3]
    F_statistic, p_value = gen_3samples(*params)
    n, bins, patches = ax[0,0].hist(F_statistic, bins = 50)
    F_critical = sp.stats.f.ppf(.95, 2, size1+size2+size3-3)
    textstr = '\n'.join((
        '$\mu_1, \mu_2, \mu_3 = {}, {}, {}$'.format(mu1, mu2, mu3),
        '$\sigma_1, \sigma_2, \sigma_3 = {}, {}, {}$'.format(sig1, sig2, sig3),
        '$n_1, n_2, n_3 = {}, {}, {}$'.format(size1, size2, size3),
        r'$F_c = {:.4f}$'.format(F_critical)))
    props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
    ax[0,0].text(max(bins)/2, max(n)/2, textstr, fontsize=10,
            verticalalignment='top', bbox=props)
    ax[0,0].set_title('Simulation 1')
    ax[0,0].vlines(F_critical, 0, max(n)*1.1, color = 'r')
    ax[1,0].hist(p_value,bins = 50)

    mu1, mu2, mu3, sig1, sig2, sig3, size1, size2, size3 = [3, 3.1, 2.9, 6, 6, 6, 10, 20, 30]
    params = [mu1, mu2, mu3, sig1, sig2, sig3, size1, size2, size3]
    F_statistic, p_value = gen_3samples(*params)
    n, bins, patches = ax[0,1].hist(F_statistic, bins = 50)
    F_critical = sp.stats.f.ppf(.95, 2, size1+size2+size3-3)
    textstr = '\n'.join((
        '$\mu_1, \mu_2, \mu_3 = {}, {}, {}$'.format(mu1, mu2, mu3),
        '$\sigma_1, \sigma_2, \sigma_3 = {}, {}, {}$'.format(sig1, sig2, sig3),
        '$n_1, n_2, n_3 = {}, {}, {}$'.format(size1, size2, size3),
        r'$F_c = {:.4f}$'.format(F_critical)))
    props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
    ax[0,1].text(max(bins)/2, max(n)/2, textstr, fontsize=10,
            verticalalignment='top', bbox=props)
    ax[0,1].vlines(F_critical, 0, max(n)*1.1, color = 'r')
    ax[1,1].hist(p_value,bins = 50)
    ax[0,1].set_title('Simulation 2')


    mu1, mu2, mu3, sig1, sig2, sig3, size1, size2, size3 = [3, 3.1, 2.9, 6, 12, 18, 10, 20, 30]
    params = [mu1, mu2, mu3, sig1, sig2, sig3, size1, size2, size3]
    F_statistic, p_value = gen_3samples(*params)
    n, bins, patches = ax[0,2].hist(F_statistic, bins = 50)
    F_critical = sp.stats.f.ppf(.95, 2, size1+size2+size3-3)
    textstr = '\n'.join((
        '$\mu_1, \mu_2, \mu_3 = {}, {}, {}$'.format(mu1, mu2, mu3),
        '$\sigma_1, \sigma_2, \sigma_3 = {}, {}, {}$'.format(sig1, sig2, sig3),
        '$n_1, n_2, n_3 = {}, {}, {}$'.format(size1, size2, size3),
        r'$F_c = {:.4f}$'.format(F_critical)))
    props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
    ax[0,2].text(max(bins)/2, max(n)/2, textstr, fontsize=10,
            verticalalignment='top', bbox=props)
    ax[0,2].vlines(F_critical, 0, max(n)*1.1, color = 'r')
    ax[1,2].hist(p_value,bins = 50)
    ax[0,2].set_title('Simulation 3')


    mu1, mu2, mu3, sig1, sig2, sig3, size1, size2, size3 =  [3, 6, 9, 10, 10, 10, 10, 20, 30]
    params = [mu1, mu2, mu3, sig1, sig2, sig3, size1, size2, size3]
    F_statistic, p_value = gen_3samples(*params)
    n, bins, patches = ax[2,0].hist(F_statistic, bins = 50)
    F_critical = sp.stats.f.ppf(.95, 2, size1+size2+size3-3)
    textstr = '\n'.join((
        '$\mu_1, \mu_2, \mu_3 = {}, {}, {}$'.format(mu1, mu2, mu3),
        '$\sigma_1, \sigma_2, \sigma_3 = {}, {}, {}$'.format(sig1, sig2, sig3),
        '$n_1, n_2, n_3 = {}, {}, {}$'.format(size1, size2, size3),
        r'$F_c = {:.4f}$'.format(F_critical)))
    props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
    ax[2,0].text(max(bins)/2, max(n)/2, textstr, fontsize=10,
            verticalalignment='top', bbox=props)
    ax[2,0].vlines(F_critical, 0, max(n)*1.1, color = 'r')
    ax[3,0].hist(p_value,bins = 50)
    ax[2,0].set_title('Simulation 4')


    mu1, mu2, mu3, sig1, sig2, sig3, size1, size2, size3 =  [3, 5, 6, 10, 10, 10, 10, 10, 10]
    params = [mu1, mu2, mu3, sig1, sig2, sig3, size1, size2, size3]
    F_statistic, p_value = gen_3samples(*params)
    n, bins, patches = ax[2,1].hist(F_statistic, bins = 50)
    F_critical = sp.stats.f.ppf(.95, 2, size1+size2+size3-3)
    textstr = '\n'.join((
        '$\mu_1, \mu_2, \mu_3 = {}, {}, {}$'.format(mu1, mu2, mu3),
        '$\sigma_1, \sigma_2, \sigma_3 = {}, {}, {}$'.format(sig1, sig2, sig3),
        '$n_1, n_2, n_3 = {}, {}, {}$'.format(size1, size2, size3),
        r'$F_c = {:.4f}$'.format(F_critical)))
    props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
    ax[2,1].text(max(bins)/2, max(n)/2, textstr, fontsize=10,
            verticalalignment='top', bbox=props)
    ax[2,1].vlines(F_critical, 0, max(n)*1.1, color = 'r')
    ax[3,1].hist(p_value,bins = 50)
    ax[2,1].set_title('Simulation 5')


    mu1, mu2, mu3, sig1, sig2, sig3, size1, size2, size3 = [3, 5, 6, 10, 10, 10, 5000, 5000, 5000]
    params = [mu1, mu2, mu3, sig1, sig2, sig3, size1, size2, size3]
    F_statistic, p_value = gen_3samples(*params)
    n, bins, patches = ax[2,2].hist(F_statistic, bins = 50)
    F_critical = sp.stats.f.ppf(.95, 2, size1+size2+size3-3)
    textstr = '\n'.join((
        '$\mu_1, \mu_2, \mu_3 = {}, {}, {}$'.format(mu1, mu2, mu3),
        '$\sigma_1, \sigma_2, \sigma_3 = {}, {}, {}$'.format(sig1, sig2, sig3),
        '$n_1, n_2, n_3 = {}, {}, {}$'.format(size1, size2, size3),
        r'$F_c = {:.4f}$'.format(F_critical)))
    props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
    ax[2,2].text(max(bins)/2, max(n)/2, textstr, fontsize=10,
            verticalalignment='top', bbox=props)
    ax[2,2].vlines(F_critical, 0, max(n)*1.1, color = 'r')
    ax[3,2].hist(p_value,bins = 50)
    ax[2,2].set_title('Simulation 6')


    mu1, mu2, mu3, sig1, sig2, sig3, size1, size2, size3 = [3, 3, 3, 100, 100, 100, 10, 10, 10]
    params = [mu1, mu2, mu3, sig1, sig2, sig3, size1, size2, size3]
    F_statistic, p_value = gen_3samples(*params)
    n, bins, patches = ax[4,0].hist(F_statistic, bins = 50)
    F_critical = sp.stats.f.ppf(.95, 2, size1+size2+size3-3)
    textstr = '\n'.join((
        '$\mu_1, \mu_2, \mu_3 = {}, {}, {}$'.format(mu1, mu2, mu3),
        '$\sigma_1, \sigma_2, \sigma_3 = {}, {}, {}$'.format(sig1, sig2, sig3),
        '$n_1, n_2, n_3 = {}, {}, {}$'.format(size1, size2, size3),
        r'$F_c = {:.4f}$'.format(F_critical)))
    props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
    ax[4,0].text(max(bins)/2, max(n)/2, textstr, fontsize=10,
            verticalalignment='top', bbox=props)
    ax[4,0].vlines(F_critical, 0, max(n)*1.1, color = 'r')
    ax[5,0].hist(p_value,bins = 50)
    ax[4,0].set_title('Simulation 7')


    mu1, mu2, mu3, sig1, sig2, sig3, size1, size2, size3 = [3, 3, 3, 1, 1, 2, 10, 20, 30]
    params = [mu1, mu2, mu3, sig1, sig2, sig3, size1, size2, size3]
    F_statistic, p_value = gen_3samples(*params)
    n, bins, patches = ax[4,1].hist(F_statistic, bins = 50)
    F_critical = sp.stats.f.ppf(.95, 2, size1+size2+size3-3)
    textstr = '\n'.join((
        '$\mu_1, \mu_2, \mu_3 = {}, {}, {}$'.format(mu1, mu2, mu3),
        '$\sigma_1, \sigma_2, \sigma_3 = {}, {}, {}$'.format(sig1, sig2, sig3),
        '$n_1, n_2, n_3 = {}, {}, {}$'.format(size1, size2, size3),
        r'$F_c = {:.4f}$'.format(F_critical)))
    props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
    ax[4,1].text(max(bins)/2, max(n)/2, textstr, fontsize=10,
            verticalalignment='top', bbox=props)
    ax[4,1].vlines(F_critical, 0, max(n)*1.1, color = 'r')
    ax[5,1].hist(p_value,bins = 50)
    ax[4,1].set_title('Simulation 8')



    params = [3, 3.1, 2.9, .01, .01, .01, 10, 20, 30]
    mu1, mu2, mu3, sig1, sig2, sig3, size1, size2, size3 = [3, 3, 3, 1, 1, 2, 10, 20, 30]
    params = [mu1, mu2, mu3, sig1, sig2, sig3, size1, size2, size3]
    F_statistic, p_value = gen_3samples(*params)
    n, bins, patches = ax[4,2].hist(F_statistic, bins = 50)
    F_critical = sp.stats.f.ppf(.95, 2, size1+size2+size3-3)
    textstr = '\n'.join((
        '$\mu_1, \mu_2, \mu_3 = {}, {}, {}$'.format(mu1, mu2, mu3),
        '$\sigma_1, \sigma_2, \sigma_3 = {}, {}, {}$'.format(sig1, sig2, sig3),
        '$n_1, n_2, n_3 = {}, {}, {}$'.format(size1, size2, size3),
        r'$F_c = {:.4f}$'.format(F_critical)))
    props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
    ax[4,2].text(max(bins)/2, max(n)/2, textstr, fontsize=10,
            verticalalignment='top', bbox=props)
    ax[4,2].vlines(F_critical, 0, max(n)*1.1, color = 'r')
    ax[5,2].hist(p_value,bins = 50)
    ax[4,2].set_title('Simulation 9')

    #######################||Rectangle||##########################
    ##############################################################

    rect = plt.Rectangle(
        # (lower-left corner), width, height
        (0.10, 0.633), 0.2645, 0.258, fill=False, color="k", lw=2, 
        zorder=1000, transform=fig.transFigure, figure=fig
    )
    fig.patches.extend([rect])


    rect = plt.Rectangle(
        # (lower-left corner), width, height
        (0.3755, 0.633), 0.2645, 0.258, fill=False, color="k", lw=2, 
        zorder=1000, transform=fig.transFigure, figure=fig
    )
    fig.patches.extend([rect])


    rect = plt.Rectangle(
        # (lower-left corner), width, height
        (0.650, 0.633), 0.2645, 0.258, fill=False, color="k", lw=2, 
        zorder=1000, transform=fig.transFigure, figure=fig
    )
    fig.patches.extend([rect])


    rect = plt.Rectangle(
        # (lower-left corner), width, height
        (0.10, 0.37), 0.2645, 0.258, fill=False, color="k", lw=2, 
        zorder=1000, transform=fig.transFigure, figure=fig
    )
    fig.patches.extend([rect])


    rect = plt.Rectangle(
        # (lower-left corner), width, height
        (0.3755, 0.37), 0.2645, 0.258, fill=False, color="k", lw=2, 
        zorder=1000, transform=fig.transFigure, figure=fig
    )
    fig.patches.extend([rect])

    rect = plt.Rectangle(
        # (lower-left corner), width, height
        (0.650, 0.37), 0.2645, 0.258, fill=False, color="k", lw=2, 
        zorder=1000, transform=fig.transFigure, figure=fig
    )
    fig.patches.extend([rect])

    rect = plt.Rectangle(
        # (lower-left corner), width, height
        (0.650, 0.108), 0.2645, 0.258, fill=False, color="k", lw=2, 
        zorder=1000, transform=fig.transFigure, figure=fig
    )
    fig.patches.extend([rect])

    rect = plt.Rectangle(
        # (lower-left corner), width, height
        (0.3755, 0.108), 0.2645, 0.258, fill=False, color="k", lw=2, 
        zorder=1000, transform=fig.transFigure, figure=fig
    )
    fig.patches.extend([rect])

    rect = plt.Rectangle(
        # (lower-left corner), width, height
        (0.1, 0.108), 0.2645, 0.258, fill=False, color="k", lw=2, 
        zorder=1000, transform=fig.transFigure, figure=fig
    )
    fig.patches.extend([rect])
    ####################################################################

    plt.show()
###############################################################################################
###############################################################################################
def two_tail_rej_region_demo():
    data = pd.read_csv('500_Person_Gender_Height_Weight_Index.csv')
    
    df = len(data[data['Gender']=='Male'])-1
    t_975 = sp.stats.t.ppf(.975, df=df)
    t_025 = sp.stats.t.ppf(.025, df=df)
    
    male_mean = data[data['Gender']=='Male']['Height'].mean()
    male_std = data[data['Gender']=='Male']['Height'].std(ddof=1)
    male_std_error = male_std/np.sqrt(len(data[data['Gender']=='Male']))
    male_null = 172
    
    x = np.linspace(male_null-5, male_null+5, 200)
    
    
    y_t = sp.stats.t.pdf(x, df = df, loc = male_null)

    fig, ax = plt.subplots(2, 1, figsize = (18,8))

    ax[0].plot(x, y_t, color = 'tomato', lw = 3)

    rejection_lower = male_null - t_975*male_std_error
    x_rej_lower = np.linspace(rejection_lower-3, rejection_lower, 30)
    y_rej_lower = sp.stats.t.pdf(x_rej_lower, df = df, loc = male_null)
    ax[0].fill_between(x_rej_lower, y_rej_lower, color = 'tomato', alpha = .7)

    rejection_upper = male_null + t_975*male_std_error
    x_rej_upper = np.linspace(rejection_upper, rejection_upper+3, 30)
    y_rej_upper = sp.stats.t.pdf(x_rej_upper, df = df, loc = male_null)
    ax[0].fill_between(x_rej_upper, y_rej_upper, color = 'tomato', alpha = .7)
    ax[0].set_ylim([0, .45])
    ax[0].set_title('Rejection Region of Original Unit (cm)')

    x = np.linspace(-5, 5, 200)
    y_t = sp.stats.t.pdf(x, df = df, loc = 0)

    ax[1].plot(x, y_t, color = 'tomato', lw = 3)

    x_rej_lower = np.linspace(t_025-3, t_025, 30)
    y_rej_lower = sp.stats.t.pdf(x_rej_lower, df = df)
    ax[1].fill_between(x_rej_lower, y_rej_lower, color = 'tomato', alpha = .7)

    x_rej_lower = np.linspace(t_975+3, t_975, 30)
    y_rej_lower = sp.stats.t.pdf(x_rej_lower, df = df)
    ax[1].fill_between(x_rej_lower, y_rej_lower, color = 'tomato', alpha = .7)
    ax[1].set_ylim([0, .45])
    ax[1].set_title('Rejection Region of t-statistic')

    plt.show()
###############################################################################################
###############################################################################################
def one_tail_rej_region_demo():
    fig, ax = plt.subplots(2, 1, figsize = (18,8))
    x = np.linspace(-5, 5, 200)
    y_t = sp.stats.t.pdf(x, df = len(x), loc = 0)
    ax[0].plot(x, y_t, color = 'tomato')

    ax[0].annotate('$H_0: \mu = \mu_0$\n$H_1: \mu < \mu_0$', (-5, .35), size = 16)
    t_05 = sp.stats.t.ppf(.05, df = len(x))
    x_rej_lower = np.linspace(t_05, t_05-3, 30)
    y_rej_lower = sp.stats.t.pdf(x_rej_lower, df = len(x))
    ax[0].fill_between(x_rej_lower, y_rej_lower, color = 'tomato', alpha = .7)
    ax[0].set_ylim([0, .45])

    x = np.linspace(-5, 5, 200)
    y_t = sp.stats.t.pdf(x, df = len(x), loc = 0)
    ax[1].plot(x, y_t, color = 'tomato')

    ax[1].annotate('$H_0: \mu = \mu_0$\n$H_1: \mu > \mu_0$', (4, .35), size = 16)
    t_95 = sp.stats.t.ppf(.95, df = len(x))
    x_rej_lower = np.linspace(t_95, t_95+3, 30)
    y_rej_lower = sp.stats.t.pdf(x_rej_lower, df = len(x))
    ax[1].fill_between(x_rej_lower, y_rej_lower, color = 'tomato', alpha = .7)
    ax[1].set_ylim([0, .45])

    plt.show()