-
Notifications
You must be signed in to change notification settings - Fork 0
/
simulation.py
100 lines (82 loc) · 2.73 KB
/
simulation.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
import math
import statistics
import random
import numpy as np
import matplotlib.pyplot as plt
numTrials = 200000
def printResults(data):
avgCount = statistics.mean(data)
std = statistics.stdev(data)
print(" e approx:", round(avgCount,4), " actual:", round(math.e,4))
print(" std dev:", round(std,4))
print(" 95% confidence interval: (", round(avgCount - std/np.sqrt(numTrials),4), ",", round(avgCount + std/np.sqrt(numTrials),4), ")\n" )
return std
# Initial simulation: straightforward but requires many trials (N) to decrease variance
print("\nSimulation 1: Initial/Straightforward")
data1 = np.zeros(numTrials)
n = numTrials
s = 0
tempCount = 0
while (n > 0):
rand = random.random()
s += rand
# print(s)
tempCount += 1
if (s >1):
# print(tempCount)
n -=1
data1[n] = tempCount
tempCount =0
s = 0
std1 = printResults(data1)
#Improved simulation that reduces variance by using antithetic variables
print("Simulation 2: using antithetic variables to reduce variance")
data2 = np.zeros(numTrials)
n = numTrials
s1, s2 = 0, 0
N1, N2 = 0, 0
N1Less1, N2Less1 = True, True
while (n > 0):
u1 = random.random()
s1 += u1
u2 = 1-u1
s2 += u2
if(N1Less1):
N1 += 1
if(N2Less1):
N2 += 1
if (s1 >1 and s2>1):
n -=1
data2[n] = (N1 + N2)/2.0
s1, s2 = 0, 0
N1, N2 = 0, 0
N1Less1, N2Less1 = True, True
elif(s1 > 1):
N1Less1 = False
elif(s2 > 1):
N2Less1 = False
std2 = printResults(data2)
#Calculate improved efficiency
# 2e-2 observations per simulation for the improved version compared to
# e observations per simulation for the original
obs2 = round((2*math.e -2)/math.e,4)
obs1 = round((std1/std2)**2,4)
print("Calculating improved efficiency: ")
print("", obs2, "times as many observations for Simulation 2, (2e-2), compared to Simulation 1 (e).")
print("", obs1, "times as many observations for Simulation 1 to have same variance as Simulation 2.")
print(" Simulation 2 is", obs1,"/",obs2,"=", round(obs1/obs2,4), "times as efficient as Simulation 1.\n" )
#Show histogram of the datasets for both simulations
# data = np.zeros(numTrials*2)
# data = data.reshape(numTrials,2)
# for i in range(numTrials):
# data[i][0] = data1[i]
# data[i][1] = data2[i]
#plt.hist(data, bins=np.arange(2,9,0.5), density=True, label=["a","Simple Simulation"])
plt.title("Histogram of Observations for Simulation 1 and Simulation 2")
plt.hist(data1, bins=np.arange(2,8,0.25), label="Simulation 1")
plt.hist(data2, bins=np.arange(2,8,0.25), label="Simulation 2")
plt.legend(prop={'size': 10})
plt.ylabel("# of Occurences")
plt.xlabel("# of Observations")
plt.xticks(np.arange(2,8,0.5))
plt.show()