-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsolution.py
59 lines (52 loc) · 1.87 KB
/
solution.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
# Uros Bojanic 2019/0077
import math
import numpy as np
from matplotlib import pyplot as plt
def BernoullisTrials(n, p, N):
f = np.zeros(n + 1)
for simulation in range(0, N):
countSuccessfulTrials = 0
for trial in range(0, n):
if np.random.rand() < p:
countSuccessfulTrials += 1
f[countSuccessfulTrials] += 1
return f / N
def BinomialDistribution(n, p):
f = np.arange(n + 1) * 1.0
for x in f.astype(int):
f[x] = math.comb(n, x) * (p**(x)) * ((1-p)**(n-x))
return f
def PlotHistogram(a, b, N):
x = np.arange(len(a))
fig, ax = plt.subplots()
w = 0.4
bar1 = ax.bar(x-w/2, a, w, align='center', alpha=0.5, label='Simulacija')
bar2 = ax.bar(x+w/2, b, w, align='center', alpha=0.5, label='Binomna raspodela')
ax.set_xticks(x)
ax.set_title('Histogram (N={})'.format(N))
ax.set_xlabel('Broj uspeha')
ax.set_ylabel('Relativna frekvencija')
ax.legend()
plt.ylim([0, 0.3])
plt.show()
def RunSimulation(n, p, N):
# Simulate trials and binomial distribution
relativeFreqSimulation = BernoullisTrials(n=n, p=p, N=N)
relativeFreqTrue = BinomialDistribution(n=n, p=p)
np.set_printoptions(formatter={'float': lambda x: "{0:0.5f}".format(x)})
print("Simulation (n={}, p={}) for {} trials:\n{}".format(n, p, N, relativeFreqSimulation))
print("Binomial distribution (n={}, p={}) for {} trials:\n{}".format(n, p, N, relativeFreqTrue))
# Calculate and print error
absoluteError = np.abs(relativeFreqTrue - relativeFreqSimulation)
print("Absolute error:\n{}".format(absoluteError))
RMSE = math.sqrt(np.square(np.subtract(relativeFreqTrue, relativeFreqSimulation)).mean())
print("Root Mean Square Error: {}\n\n".format(RMSE))
# Plot histogram
PlotHistogram(relativeFreqSimulation, relativeFreqTrue, N)
def main():
# First simulation (100 trials)
RunSimulation(n=12, p=0.4, N=100)
# Second simulation (1000 trials)
RunSimulation(n=12, p=0.4, N=1000)
if __name__=="__main__":
main()