import numpy as np
import time
# Stat 624 lecture on simulation studies
# Written by Robert Richardson
# Setting the number of samples to be drawn from a normal distribution
n_samp = 5
# Segment 1: Using a for loop to conduct Monte Carlo simulation
# Initializing variables to store results
# Measuring the time taken for this segment of code to run using time.time()
start_time = time.time()
B = 100000 # Setting the number of Monte Carlo draws
mus = np.zeros(B) # Creating an array to store mean values
meds = np.zeros(B) # Creating an array to store median values
# Looping over B iterations
for i in range(B):
x = np.random.normal(size=n_samp) # Generating n_samp random numbers from the standard normal distribution
mus[i] = np.mean(x) # Calculating and storing the mean of the generated numbers
meds[i] = np.median(x) # Calculating and storing the median of the generated numbers
print("For loop time: ", time.time() - start_time)
# Segment 2: Using vectorization to conduct Monte Carlo simulation
# Generating an array of random numbers and using numpy functions to calculate mean and median
# Measuring the time taken for this segment of code to run using time.time()
start_time = time.time()
draws = np.random.normal(size=(B, n_samp)) # Generating a 2D array of random numbers from the normal distribution
mus = np.mean(draws, axis=1) # Calculating the mean of each row using numpy mean function
meds = np.median(draws, axis=1) # Calculating the median of each row using numpy median function
print("Vectorization time: ", time.time() - start_time)
# Calculating and printing the estimated mean of the means and its 95% confidence interval
mean_mus = np.mean(mus)
std_mus = np.std(mus)
ci_mus = mean_mus + np.array([-1, 1]) * 1.96 * std_mus / np.sqrt(B)
print(f"The means have an estimate of {mean_mus:.5f}")
print(f"The 95% confidence interval is ({ci_mus[0]:.5f}, {ci_mus[1]:.5f}).")
# Calculating and printing the estimated mean of the medians and its 95% confidence interval
mean_meds = np.mean(meds)
std_meds = np.std(meds)
ci_meds = mean_meds + np.array([-1, 1]) * 1.96 * std_meds / np.sqrt(B)
print(f"The medians have an estimate of {mean_meds:.5f}")
print(f"The 95% confidence interval is ({ci_meds[0]:.5f}, {ci_meds[1]:.5f}).")