text: 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()