Introduction to Data Science¶
1MS041, 2023¶
©2023 Raazesh Sainudiin, Benny Avelin. Attribution 4.0 International (CC BY 4.0)
Estimation of mean¶
- Load data
- Make assumptions
- Estimate mean
- Write down confidence interval
In [66]:
%%bash
ls data
In [67]:
%%bash
head -n 2 data/NYPowerBall.csv
In [68]:
import csv
with open("data/NYPowerBall.csv",mode='r') as f:
reader = csv.reader(f)
header = next(reader)
#data = [i for i in reader]
data = list(reader)
In [69]:
data[:2]
Out[69]:
In [70]:
list_list_list = [line[1].split(' ') for line in data]
In [71]:
list_list_list[:2]
Out[71]:
In [72]:
flattened_list = sum(list_list_list,[])
In [73]:
import numpy as np
list_list_list_arr = np.array(list_list_list)
In [75]:
list_list_list_arr.flatten().shape
Out[75]:
In [77]:
int_list = [int(i) for i in flattened_list]
In [78]:
int_list[:2]
Out[78]:
In [79]:
int_arr = np.array(int_list)
In [80]:
int_arr.shape
Out[80]:
In [81]:
np.max(int_arr)
Out[81]:
In [82]:
np.min(int_arr)
Out[82]:
In [84]:
len(np.unique(int_arr))
Out[84]:
Assumption about data¶
- IID (Independent and identically distributed), is this really true?
- The n:th draw of the m:th trial is IID but not the n and n-1 for instance.
- Discrete distribution with all numbers between 1 and 69.
What to do¶
Compute confidence interval of what?
- The mean of the first draw
- Hoeffdings inequality, $a = 1$ and $b=69$
In [85]:
reshaped_int_arr = int_arr.reshape(-1,6)
In [87]:
first_draws = reshaped_int_arr[:,0]
In [94]:
list_list_list[:10]
Out[94]:
In [91]:
np.unique(first_draws)
Out[91]:
In [100]:
# Double check
np.unique([int(i[5]) for i in list_list_list])
Out[100]:
In [101]:
power_ball = reshaped_int_arr[:,-1]
In [102]:
power_ball_mean = np.mean(power_ball)
power_ball_mean
Out[102]:
In [103]:
39/2
Out[103]:
Hoeffdings inequality $P(X \in [a,b]) = 1$ $$ P(|\overline{X}_n - E[X]| > \epsilon) \leq 2 e^{-\frac{2 n \epsilon^2}{(a-b)^2}} $$
$$
\alpha = 2 e^{-\frac{2 n \epsilon^2}{(a-b)^2}}
$$$$
\epsilon = \sqrt{\frac{(a-b)^2}{2n}\ln(2/\alpha)}
$$
In [104]:
def compute_epsilon(n,a,b,alpha):
return np.sqrt((a-b)**2/(2*n) * np.log(2/alpha))
In [105]:
alpha = 0.05
a = 1
b = 39
n = len(power_ball)
epsilon = compute_epsilon(n,a,b,alpha)
epsilon
Out[105]:
In [106]:
(power_ball_mean-epsilon,power_ball_mean+epsilon)
Out[106]:
In [107]:
(39+1)/2
Out[107]:
In [108]:
from Utils import discrete_histogram
discrete_histogram(power_ball)
Likelihood of parameter¶
- Load data
- Work with dates
- Make assumptions
- Write down the Risk
- Split into two parts
- Minimize the risk numerically on train
- Test the risk on test
- Make confidence intervals
In [109]:
%%bash
ls data
In [111]:
%%bash
head -n 2 data/earthquakes.csv
In [115]:
import csv
with open("data/earthquakes.csv",mode='r') as f:
reader = csv.reader(f,skipinitialspace=True)
header = next(reader)
#data = [i for i in reader]
data = list(reader)
In [128]:
print(header[2],data[0][2])
In [124]:
format_string = "%Y-%m-%dT%H:%M:%S.%fZ"
In [127]:
from datetime import datetime
datetime.strptime(data[0][2],format_string)
Out[127]:
In [129]:
origin_time = [datetime.strptime(line[2],format_string) for line in data]
In [132]:
(origin_time[0]-origin_time[1]).total_seconds()
Out[132]:
In [133]:
origin_time_arr = np.array(origin_time)
In [135]:
sorted_origin_time = np.sort(origin_time_arr)
In [136]:
sorted_origin_time[1:]-sorted_origin_time[:-1]
Out[136]:
In [138]:
time_between = np.diff(sorted_origin_time)
In [139]:
time_between_seconds = [time.total_seconds() for time in time_between]
In [143]:
import matplotlib.pyplot as plt
_=plt.hist(time_between_seconds,bins=200)
We construct the model space $\mathcal{M} = \{f_\lambda(x) = \lambda e^{-\lambda x}: \lambda > 0\}$
Now define the log-loss $$ L(f_\lambda,x) = -\ln(f_\lambda(x)) $$
The risk is the expected loss $$ R(f_\lambda) = E[L(f_\lambda,X)] $$
We instead minimize the empirical risk, given i.i.d. Data $X = \{X_1,\ldots,X_n\}$ $$ \hat R_n(f_\lambda) = \frac{1}{n} \sum_{i=1}^n L(f_\lambda,X_i) $$
The goal is to minimize $\hat R_n(f_\lambda)$ w.r.t. $\lambda$.
$$ \ln(f_\lambda) = \ln(\lambda) - \lambda x $$
$$
\hat R_n(f_\lambda) = \frac{1}{n} \sum_{i=1}^n (\lambda X_i - \ln(\lambda))
$$
$$
\frac{d}{d\lambda}
\hat R_n(f_\lambda) = \frac{1}{n} \sum_{i=1}^n (X_i - 1/\lambda)
$$
$$
\frac{d}{d\lambda}
\hat R_n(f_\lambda) = 0
$$
gives $$ 1/\hat \lambda = \frac{1}{n} \sum_{i=1}^n X_i $$
In [146]:
import matplotlib.pyplot as plt
_=plt.hist(time_between_seconds,bins=200,density=True)
hat_lambda = 1/np.mean(time_between_seconds)
x_plot = np.linspace(0,20000,100)
plt.plot(x_plot,hat_lambda*np.exp(-hat_lambda*x_plot))
Out[146]:
Lets say I want to predict the time to next earthquake, lets say that we use our estimated $\hat \lambda$ to predict $1/\hat \lambda$.
Train, testing split. But for us, using the iid assumption we can just randomly split our dataset into two parts.
In [147]:
n_train = int(len(time_between_seconds)/2)
train_set = time_between_seconds[:n_train]
test_set = time_between_seconds[n_train:]
In [148]:
hat_lambda_train = 1/np.mean(train_set)
guess = 1/hat_lambda_train
In [149]:
guess
Out[149]:
In [152]:
_=plt.hist(np.abs(guess-test_set),bins=100)
In [153]:
average_error = np.mean(np.abs(guess-test_set))
print(average_error)