In this problem, we are given a dataset which has two columns 'BRT'
and 'metro'
and each row shows how many BRTs/trains have arrived at the given range of time (like an hour for example). Let's answer each question in the problem:
Drawing the histograms for BRT (y
) and metro (x
) datasets:
import numpy as np
import scipy.stats as stats
import pandas as pd
import matplotlib.pyplot as plt
green = '#57cc99'
red = '#e56b6f'
blue = '#22577a'
yellow = '#ffca3a'
bar_width = 0.5
TARBIAT_CSV_FILE_PATH = "../given-files/Tarbiat.csv"
data = pd.read_csv(TARBIAT_CSV_FILE_PATH)
fig, (ax1, ax2) = plt.subplots(1, 2)
x = data["metro"]
y = data["BRT"]
ax1.hist(x, bins='auto', alpha=0.7, label='Metro', color=red, edgecolor='black', width=bar_width)
ax1.set_title('Metro')
ax1.grid(True)
ax1.legend()
ax1.set_xlabel('Number of arrives')
ax1.set_ylabel('Freq')
ax2.hist(y, bins='auto', alpha=0.7, label='BRT', color=yellow, edgecolor='black', width=bar_width)
ax2.set_title('BRT')
ax2.grid(True)
ax2.legend()
ax2.set_xlabel('Number of arrives')
ax2.set_ylabel('Freq')
plt.show()
We know that for events that happen not so frequently in time, poisson is a good prediction distribution (
A poisson distribution can be fully described with it's parameter, the x
and y
datasets and create two poisson variables with calculated $\lambda$s:
x_values = np.arange(0, 10)
x_mean = np.mean(x)
y_mean = np.mean(y)
poisson_x = stats.poisson(x_mean)
pmf_x = poisson_x.pmf(x_values)
poisson_y = stats.poisson(y_mean)
pmf_y = poisson_y.pmf(x_values)
Now, let's draw the denisity histograms for our datasets and
fig, (ax1, ax2) = plt.subplots(1, 2)
counts1, bins1, _ = plt.hist(x, bins='auto')
counts2, bins2, _ = plt.hist(y, bins='auto')
plt.cla()
ax1.bar(x_values + bar_width, pmf_x, color=blue, label='X ~ Pois({:.2f})'.format(x_mean), alpha=0.7, width=bar_width)
ax1.bar(bins1[:-1], counts1 / len(x), width=bar_width, alpha=0.7, label='Metro', color=red, edgecolor='black')
ax1.set_title('Metro')
ax1.grid(True)
ax1.legend()
ax2.bar(x_values + bar_width, pmf_y, color=blue, label='Y ~ Pois({:.2f})'.format(y_mean), alpha=0.7, width=bar_width)
ax2.bar(bins2[:-1], counts2 / len(y), width=bar_width, alpha=0.7, label='BRT', color=yellow, edgecolor='black')
ax2.set_title('BRT')
ax2.grid(True)
ax2.legend()
plt.show()
Great! as we can see, it is a good prediction model for our dataset.
We know that if z = x + y
:
fig, ax = plt.subplots()
z = x + y
counts3, bins3, _ = plt.hist(z, bins='auto')
plt.cla()
poisson_z = stats.poisson(x_mean + y_mean)
pmf_z = poisson_z.pmf(x_values)
ax.bar(x_values + bar_width, pmf_z, color=blue, label='Z ~ Pois({:.2f})'.format(x_mean + y_mean), width=bar_width, alpha=0.7)
ax.bar(bins3[:-1], counts3 / len(z), width=bar_width, alpha=0.7, label='BRT + Metro', color=green, edgecolor='black')
ax.set_title('BRT + Metro')
ax.grid(True)
ax.legend()
plt.show()
As you can see, the pridiction model is again to good fit for our dataset.
Let's hand-calculate the pmf of
Applying the conditional probability formula:
We know that for a random variable
Applying that to our formula we'll get:
Which is the pmf of a binomial variable with
Let's draw
fig, ax = plt.subplots()
p = x_mean / (x_mean + y_mean)
n = 8
w_counts = [0] * (9)
x_values = np.arange(len(w_counts))
pmf_w = stats.binom.pmf(x_values, n, p)
for _, row in data.iterrows():
if row['metro'] + row['BRT'] == n:
w_counts[row['metro']] += 1
count = np.count_nonzero((data['metro'] + data['BRT'] == n))
w_counts = [x / count for x in w_counts]
ax.bar(x_values + bar_width, w_counts, color=red, label='Conditional probability in action', alpha=0.7, edgecolor='black', width=bar_width)
ax.bar(x_values, pmf_w, color=blue, label='W ~ Bin(8, {:.2f})'.format(p), alpha=0.7, width=bar_width, edgecolor='black')
ax.set_title('Conditional probability')
ax.grid(True)
ax.legend()
plt.show()
happy happy happy
In this problem, we must calculate answer of coupon collector problem with two methods:
In this method, we perform the test k
times and return the average of answers as the final result:
import random
def monte_carlo(n):
k = 1000
ans = 0
for _ in range(k):
seen = 0
trial = 0
while seen != n:
trial += 1
if random.random() <= (n - seen) / n:
seen += 1
ans += trial
return ans / k
In this method, we define
The sum of all $X_i$s is our answer's random variable:
It's clear that
We could easily calculate
The generating function for
To calculate
And since expectation in linear:
Let's write the code for this:
import sympy as sp
def theorical_solution(n):
n = int(n)
expectation = 0
s = sp.symbols('s')
x = sp.symbols('x')
for i in range(1, n + 1):
p = (n - i + 1) / n
gen_func = p * sp.exp(s) / (1 - (1 - p) * sp.exp(s))
deriv = sp.diff(gen_func, s)
expectation += deriv.subs({s:0})
return expectation
And now let's draw the diagrams of our two methods:
dom = 20
x = np.linspace(1, dom, dom)
mc = [monte_carlo(n) for n in x]
ts = [theorical_solution(n) for n in x]
plt.plot(x, mc, 'o', label='Monte Carlo Approach', color=blue, alpha=0.5, markersize=10)
plt.plot(x, ts, 'o', label='Theorical Approach', color=red, alpha=0.5, markersize=10)
plt.xlabel('n')
plt.ylabel('Solution for n')
plt.title('Coupon collector’s problem')
plt.legend()
plt.grid(True)
plt.show()
Let's calculate
Let's draw
line_thickness = 2
dom = 60
x = np.linspace(1, dom, dom)
mc = [monte_carlo(n) for n in x]
plt.plot(x, mc, 'o', label='Monte Carlo Approach', color=blue, alpha=0.7)
plt.plot(x, x * np.log(x), label='nlog(n)', color=green, linewidth=line_thickness)
plt.xlabel('n')
plt.ylabel('Solution for n')
plt.title('Coupon collector’s problem')
plt.legend()
plt.grid(True)
plt.show()
Here's the code for part 1:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
from time import sleep
from IPython import display
t = 1000
p = np.linspace(0,1,t)
fy = stats.beta.pdf(p, a=1, b=1)
DIGITS_CSV_FILE_PATH = "../given-files/digits.csv"
df = pd.read_csv(DIGITS_CSV_FILE_PATH)
r201 = df.iloc[200]
r202 = df.iloc[201]
df.drop(200, inplace=True)
df.drop(201, inplace=True)
We should replace all color numbers
df.loc[:, df.columns != 'label'] = df.loc[:, df.columns != 'label'].map(lambda x: 0 if x < 128 else 1)
I'll choose row number 87:
chosen_row = df.iloc[87].values[1:].reshape(28, 28)
plt.imshow(chosen_row)
plt.show()
happy happy happy
This part is the main part. Our program should update fy
) based on it's previous state (fy
), pny
) and integral
:
def update(fy: np.array, n: bool) -> np.array:
p = np.linspace(0,1,t)
pny = p if n else 1 - p
integral = 0
for i in range(t):
integral += fy[i] * pny[i] / 1000
post = fy * pny / integral
return post
plt.figure(figsize=(10,8))
for i in range(100):
n = df[df['label'] == 8].iloc[i, df.columns.get_loc('pixel404')]
fy = update(fy, n)
plt.plot(p, fy, 'r', label='1')
plt.ylim(-1, 10)
plt.xlim(0, 1)
plt.text(0.1,9,f'number of seen data : {i + 1}, p = {fy.argmax() / t :.2f}', color='purple')
display.clear_output(wait=True)
display.display(plt.gcf())
plt.clf()
sleep(0.05)
<Figure size 1000x800 with 0 Axes>
We can see that as we observe more data our