import numpy as np import statsmodels.api as sm import scipy.stats as stats import matplotlib.pyplot as plt from numpy import random #from scipy.stats import chi2 # set the random seed: np.random.seed(123456) # set sample size: n_2 = 2 n_10 = 10 n_20 = 20 n_30 = 30 # initialize ybar to an array of length r=10000 to later store results: r = 10000 ybar_2 = np.empty(r) ybar_10 = np.empty(r) ybar_20 = np.empty(r) ybar_30 = np.empty(r) # repeat r times: for j in range(r): # draw a sample and store the sample mean in pos. j=0,1,... of ybar: sample_2 = random.chisquare(df=1, size=n_2) ybar_2[j] = np.mean(sample_2) for j in range(r): # draw a sample and store the sample mean in pos. j=0,1,... of ybar: sample_10 = random.chisquare(df=1, size=n_10) ybar_10[j] = np.mean(sample_10) for j in range(r): # draw a sample and store the sample mean in pos. j=0,1,... of ybar: sample_20 = random.chisquare(df=1, size=n_20) ybar_20[j] = np.mean(sample_20) for j in range(r): # draw a sample and store the sample mean in pos. j=0,1,... of ybar: sample_30 = random.chisquare(df=1, size=n_30) ybar_30[j] = np.mean(sample_30) mean_2 = np.mean(ybar_2) variance_2 = np.var(ybar_2) print("Mean of sampling distribution w/ df=1 & n=2 is :", mean_2) print("Variance of sampling distribution w/ df=1 & n=2 is :", variance_2) mean_10 = np.mean(ybar_10) variance_10 = np.var(ybar_10) print("Mean of sampling distribution w/ df=1 & n=10 is :", mean_10) print("Variance of sampling distribution w/ df=1 & n=10 is :", variance_10) mean_20 = np.mean(ybar_20) variance_20 = np.var(ybar_20) print("Mean of sampling distribution w/ df=1 & n=20 is :", mean_20) print("Variance of sampling distribution w/ df=1 & n=20 is :", variance_20) mean_30 = np.mean(ybar_30) variance_30 = np.var(ybar_30) print("Mean of sampling distribution w/ df=1 & n=30 is :", mean_30) print("Variance of sampling distribution w/ df=1 & n=30 is :", variance_30) # simulated density: kde_2 = sm.nonparametric.KDEUnivariate(ybar_2) kde_2.fit() kde_10 = sm.nonparametric.KDEUnivariate(ybar_10) kde_10.fit() kde_20 = sm.nonparametric.KDEUnivariate(ybar_20) kde_20.fit() kde_30 = sm.nonparametric.KDEUnivariate(ybar_30) kde_30.fit() # normal density: x_range_2 = np.linspace(0, 10) y_2 = stats.norm.pdf(x_range_2, 1, np.sqrt(2/n_2)) x_range_10 = np.linspace(0, 4) y_10 = stats.norm.pdf(x_range_10, 1, np.sqrt(2/n_10)) x_range_20 = np.linspace(0, 3) y_20 = stats.norm.pdf(x_range_20, 1, np.sqrt(2/n_20)) x_range_30 = np.linspace(0, 2) y_30 = stats.norm.pdf(x_range_30, 1, np.sqrt(2/n_30)) # create graph: fig = plt.figure(figsize=(7,7)) plt.subplot(221) _ = plt.plot(kde_2.support, kde_2.density, color='black', label='ybar_2') _ = plt.plot(x_range_2, y_2, linestyle='--', color='black', label='chi-square distribution') _ = plt.ylim(0, 0.8) _ = plt.ylabel('density') _ = plt.xlabel('ybar') _ = plt.legend() plt.subplot(222) _ = plt.plot(kde_10.support, kde_10.density, color='black', label='ybar_10') _ = plt.plot(x_range_10, y_10, linestyle='--', color='black', label='chi-square distribution') _ = plt.ylim(0, 1.0) _ = plt.ylabel('density') _ = plt.xlabel('ybar') _ = plt.legend() plt.subplot(223) _ = plt.plot(kde_20.support, kde_20.density, color='black', label='ybar_20') _ = plt.plot(x_range_20, y_20, linestyle='--', color='black', label='chi-square distribution') _ = plt.ylim(0, 1.5) _ = plt.ylabel('density') _ = plt.xlabel('ybar') _ = plt.legend() plt.subplot(224) _ = plt.plot(kde_30.support, kde_30.density, color='black', label='ybar_30') _ = plt.plot(x_range_30, y_30, linestyle='--', color='black', label='chi-square distribution') _ = plt.ylim(0, 2.0) _ = plt.ylabel('density') _ = plt.xlabel('ybar') _ = plt.legend() fig.subplots_adjust(hspace=.3, wspace=0.4) # 하위그림의 위아래 간격을 조정 #plt.savefig('C:/BOOK/PyBasics/PyStat/code/b1-ch5-5.png')