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


###############################################################################################################
###############################################################################################################
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, 0.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=0.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, -0.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 - 0.6 * lrg2.x, color="#FA954D", alpha=0.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=0.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=0.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=0.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=0.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=0.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=0.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=0.9
    )  # vertical line
    ax[0, 0].plot([x[-1], x[-1]], [0, 1 / (b - a)], lw=3, color="r", alpha=0.9)
    ax[0, 0].fill_between(x, 1 / (b - a), 0, alpha=0.5, color="r")

    ax[0, 0].set_xlim([1, 9])
    ax[0, 0].set_ylim([0, 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=0.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=0.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=0.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=0.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=0.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=0.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=0.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=0.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=0.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=0.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=0.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 - 0.2, max(y_null) / 2), size=15)
    ax.annotate("Alternative", (alter_loc - 0.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, 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(0.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(0.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(0.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(0.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(0.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(0.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(0.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(0.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, 0.01, 0.01, 0.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(0.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(0.975, df=df)
    t_025 = sp.stats.t.ppf(0.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=0.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=0.7)
    ax[0].set_ylim([0, 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=0.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=0.7)
    ax[1].set_ylim([0, 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, 0.35), size=16)
    t_05 = sp.stats.t.ppf(0.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=0.7)
    ax[0].set_ylim([0, 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, 0.35), size=16)
    t_95 = sp.stats.t.ppf(0.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=0.7)
    ax[1].set_ylim([0, 0.45])

    plt.show()