import numpy as np

import scipy.stats as stats

import matplotlib.pyplot as plt

# set the random seed:

np.random.seed(123456)

# set sample size and MC simulations:

r = 10000

n = 100

# initialize arrays to later store results:

CIlower = np.empty(r)

CIupper = np.empty(r)

pvalue1 = np.empty(r)

pvalue2 = np.empty(r)

# repeat r times:

for j in range(r):

# draw a sample:

sample = stats.norm.rvs(10, 2, size=n)

sample_mean = np.mean(sample)

sample_sd = np.std(sample, ddof=1)

# test the (correct) null hypothesis mu=10:

testres1 = stats.ttest_1samp(sample, popmean=10)

pvalue1[j] = testres1.pvalue

cv = stats.t.ppf(0.975, df=n - 1)

CIlower[j] = sample_mean - cv * sample_sd / np.sqrt(n)

CIupper[j] = sample_mean + cv * sample_sd / np.sqrt(n)

# test the (incorrect) null hypothesis mu=9.5 & store the p value:

testres2 = stats.ttest_1samp(sample, popmean=9.5)

pvalue2[j] = testres2.pvalue

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

## correct H0 ##

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

plt.figure(figsize=(3, 5)) # set figure ratio

plt.ylim(0, 101)

plt.xlim(9, 11)

for j in range(1, 101):

if 10 > CIlower[j] and 10 < CIupper[j]:

plt.plot([CIlower[j], CIupper[j]], [j, j], linestyle='-', color='grey')

else:

plt.plot([CIlower[j], CIupper[j]], [j, j], linestyle='-', color='black')

plt.axvline(10, linestyle='--', color='black', linewidth=0.5)

plt.ylabel('Sample No.')

plt.savefig('PyGraphs/Simulation-Inference-Figure1.pdf')

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

## incorrect H0 ##

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

plt.figure(figsize=(3, 5)) # set figure ratio

plt.ylim(0, 101)

plt.xlim(9, 11)

for j in range(1, 101):

if 9.5 > CIlower[j] and 9.5 < CIupper[j]:

plt.plot([CIlower[j], CIupper[j]], [j, j], linestyle='-', color='grey')

else:

plt.plot([CIlower[j], CIupper[j]], [j, j], linestyle='-', color='black')

plt.axvline(9.5, linestyle='--', color='black', linewidth=0.5)

plt.ylabel('Sample No.')

plt.savefig('PyGraphs/Simulation-Inference-Figure2.pdf')