Commit 413724df authored by Luke Ng's avatar Luke Ng
Browse files

Added topic Jupyter Notebook examples.

parent 43240906
Loading
Loading
Loading
Loading
+268 −0

File added.

Preview size limit exceeded, changes collapsed.

Topic02/Diabetes.csv

0 → 100644
+769 −0

File added.

Preview size limit exceeded, changes collapsed.

+125 −0
Changes for Topic02/MissingData-Correlation.ipynb: 125 added lines, 0 removed lines.
Original line number Diff line number Diff line
%% Cell type:markdown id: tags:

## Import Libraries

%% Cell type:code id: tags:

``` python
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

import seaborn as sns
```

%% Cell type:markdown id: tags:

## Load Diabetes Data

%% Cell type:code id: tags:

``` python
data = pd.read_csv("Diabetes.csv")
data.head()
```

%% Cell type:markdown id: tags:

## View Missing Values

%% Cell type:code id: tags:

``` python
# If the data uses '?' for missing values then we can replace them with a NaN
#data.replace('?', np.NaN, inplace=True)

# using isnull() function
# View each column and how many missing values it has
print(data.isnull().sum())
```

%% Cell type:code id: tags:

``` python
data.shape
```

%% Cell type:code id: tags:

``` python
374/768
```

%% Cell type:markdown id: tags:

## Deal with Missing Values

%% Cell type:code id: tags:

``` python
## Insulin has a large number of missing values .. so we can drop that column
data = data.drop('Insulin', axis=1)
```

%% Cell type:code id: tags:

``` python
### Replace missing values in each column with the mean or median of that column
#data.fillna(data.mean(), inplace=True)
data.fillna(data.median(), inplace=True)

# Or drop rows with missing values
# data = data.dropna()
```

%% Cell type:code id: tags:

``` python
print(data.isnull().sum())
```

%% Cell type:markdown id: tags:

## Scatter Matrix

%% Cell type:code id: tags:

``` python
sns.pairplot(data, hue="Outcome", diag_kind = 'hist');
```

%% Cell type:markdown id: tags:

## Correlation Matrix

%% Cell type:code id: tags:

``` python
data.corr()
```

%% Cell type:code id: tags:

``` python
plt.figure(figsize=(9,9))
sns.heatmap(np.abs(data.corr()), annot=True, cmap="viridis", fmt="0.2f");
```

%% Cell type:markdown id: tags:

## Histogram Plots

%% Cell type:code id: tags:

``` python
# https://seaborn.pydata.org/generated/seaborn.distplot.html
for i,col in enumerate(data.columns[:-1]):
    plt.figure(i)
    sns.distplot(data[col], kde=True, rug=True);
```

%% Cell type:code id: tags:

``` python
```
+170 −0
Changes for Topic04/Significance-Tests.ipynb: 170 added lines, 0 removed lines.
Original line number Diff line number Diff line
%% Cell type:markdown id: tags:

## Import Libraries

%% Cell type:code id: tags:

``` python
from random import randrange
import pandas as pd
import numpy as np
from scipy.stats import ttest_ind
from scipy.stats import ttest_1samp
from scipy.stats import chi2_contingency
from scipy.stats import f_oneway
from statsmodels.stats.power import TTestIndPower
```

%% Cell type:markdown id: tags:

## One-way One Sample t-Test (team height example from the slides in topic 4 part 2)

%% Cell type:code id: tags:

``` python
## generate some heights for 18 players in our team
## Notice that your results can be different as we're generating random numbers
team_heights = []
for i in range(18):
    team_heights.append(randrange(180,205))

print(team_heights)
print(f"Average player height in our team = {np.mean(team_heights):.5f}")
print(f"STD of player height in our team = {np.std(team_heights):.5f}")
```

%% Cell type:code id: tags:

``` python
## popmean is the mean player height in the league as was explained in the lecture
tscore, pvalue = ttest_1samp(team_heights, popmean=201)
print("t Statistic: ", tscore)
print(f'p-value for single sided test: {pvalue / 2:.5f}')
```

%% Cell type:markdown id: tags:

## One-way Two Sample t-Test (Sleep apnoea events example from the slides in topic 4 part 2)

%% Cell type:code id: tags:

``` python
## generate some numbers for people who followed the diet
## Notice that your results can be different as we're generating random numbers
diet = []
for i in range(13):
    diet.append(randrange(3,9))
print(diet)
```

%% Cell type:code id: tags:

``` python
## generate some numbers for people who did not follow the diet
## Notice that your results can be different as we're generating random numbers
nodiet = []
for i in range(14):
    nodiet.append(randrange(4,10))
print(nodiet)
```

%% Cell type:code id: tags:

``` python
## Run the t-test and extract the p-value
## We divide the resulting p-value by 2 beause this is a single-sided t-test
## Please read the documentation of the ttest_ind function
res = ttest_ind(diet, nodiet, equal_var=False)
print(f'p-value for single sided test: {res.pvalue / 2:.5f}')
```

%% Cell type:markdown id: tags:

## One-way ANOVA for multiple groups
https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.f_oneway.html

Data from: http://www.biostathandbook.com/onewayanova.html

%% Cell type:code id: tags:

``` python
tillamook = [0.0571, 0.0813, 0.0831, 0.0976, 0.0817, 0.0859, 0.0735,
             0.0659, 0.0923, 0.0836]
newport = [0.0873, 0.0662, 0.0672, 0.0819, 0.0749, 0.0649, 0.0835,
           0.0725]
petersburg = [0.0974, 0.1352, 0.0817, 0.1016, 0.0968, 0.1064, 0.105]
magadan = [0.1033, 0.0915, 0.0781, 0.0685, 0.0677, 0.0697, 0.0764,
           0.0689]
tvarminne = [0.0703, 0.1026, 0.0956, 0.0973, 0.1039, 0.1045]
res = f_oneway(tillamook, newport, petersburg, magadan, tvarminne)
print(f'p-value for one-sided ANOVA : {res.pvalue :.5f}')
```

%% Cell type:markdown id: tags:

## Chi-square test

%% Cell type:code id: tags:

``` python
a1 = [8, 6, 6, 4]
a2 = [6, 6, 6, 6]
a3 = [4, 5, 5, 7]
a4 = [5, 6, 6, 7]
a5 = [7, 7, 5, 4]
a6 = [5, 5, 7, 7]
dice = np.array([a1, a2, a3, a4, a5, a6])
df = pd.DataFrame(dice, columns=['Exp1','Exp2','Exp3','Exp4'])
df.index = ['One','Two','Three','Four','Five','Six']
#df.sum(axis=1)
```

%% Cell type:code id: tags:

``` python
df
```

%% Cell type:code id: tags:

``` python
chi2_stat, p_val, dof, ct = chi2_contingency(dice)
print("Chi2 Stat = ", chi2_stat)

print("Degrees of Freedom =", dof)

print("p-value =", p_val)

print("Contingency Table:")
print(ct)
```

%% Cell type:markdown id: tags:

## Power Analysis

* Effect Size: The minimum size of the effect that you hope to be able to detect in a statistical test, such as “a 20% improvement in click rates”.
* Statistical Power: The probability of detecting a given effect size with a given sample size (i.e. The probability of accepting the alternative hypothesis if it is true).
* Sample Size: The number of observations in the sample.
* Significance level: The statistical significance level at which the test will be conducted (e.g. alpha. which is often set to 5% or 0.05).

%% Cell type:code id: tags:

``` python
# estimate required sample size using power analysis
# parameters for power analysis
effect = 0.88
power = 0.85
alpha = 0.05

# perform power analysis
pow_analysis = TTestIndPower()
result = pow_analysis.solve_power(effect, power=power, alpha=alpha)

print(f'Sample Size: {result:.2f}')
```

%% Cell type:code id: tags:

``` python
```

Topic05/LinReg.ipynb

0 → 100644
+257 −0
Changes for Topic05/LinReg.ipynb: 257 added lines, 0 removed lines.
Original line number Diff line number Diff line
%% Cell type:markdown id: tags:

### Import Libraries

%% Cell type:code id: tags:

``` python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.stats import linregress
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
import statsmodels.api as sm
```

%% Cell type:markdown id: tags:

### Simple Data for Simple Linear Regression

%% Cell type:code id: tags:

``` python
# Independent variable
x = np.array([210,290,250,500,310,430,455,380,535,510])
# Dependent variable
y = np.array([5,7,6,13,8,11,12,10,15,14])
```

%% Cell type:markdown id: tags:

### Plot Data

%% Cell type:code id: tags:

``` python
#Plot the data along
plt.figure(figsize=(6,6), dpi=100)
plt.plot(x, y, 'o', label='original data')
plt.xlabel('Distance in KM') # X axis data label
plt.ylabel('Delivery Time in Days') # Y axis data label
plt.legend()
plt.show()
```

%% Cell type:code id: tags:

``` python
# Means of indep and dep variables
xbar = np.mean(x)
ybar = np.mean(y)
# Apply equations to find b and a
b = np.sum((x - xbar)*(y - ybar)) / np.sum((x - xbar)**2)
a = ybar - b * xbar
# Plug in the values of the dep variable into the line equation to obtain yhat
yhat = a + b * x

# Compute the error (i.e. residuals)
error = y - yhat

# Compute metrics
SE = error**2 # squared error
MSE = np.mean(SE) # mean squared error
RMSE = np.sqrt(MSE) # Root Mean Squared Error, RMSE

SST = np.sum((y - ybar)**2)
SSR = np.sum((yhat - ybar)**2)
Rsquared = SSR/SST
print('RMSE =',RMSE)
print('Rsquared =',Rsquared)
# This should give you the same result for Rsquared
#Rsquared = 1.0 - (np.var(error) / np.var(y))
```

%% Cell type:code id: tags:

``` python
#Plot the data along with the fitted line:
plt.figure(figsize=(6,6), dpi=100)
plt.plot(x, y, 'o', label='original data')
plt.plot(x, yhat, 's', label='predicted data')
plt.xlabel('Distance in KM') # X axis data label
plt.ylabel('Delivery Time in Days') # Y axis data label
plt.legend()
plt.show()
```

%% Cell type:markdown id: tags:

### Plot original points, fitted line and residuals

%% Cell type:code id: tags:

``` python
f = plt.figure(figsize=(6,6), dpi=100)
axes = f.add_subplot(111)

# first the raw data as a scatter plot
axes.plot(x, y,  'o')

# now the model as a line plot
axes.plot(x, yhat)

# now add individual line for each point
for i in range(len(x)):
    lineXdata = (x[i], x[i]) # same X
    lineYdata = (y[i], yhat[i]) # different Y
    plt.plot(lineXdata, lineYdata)

axes.set_xlabel('X Data') # X axis data label
axes.set_ylabel('Y Data') # Y axis data label

plt.show()
```

%% Cell type:markdown id: tags:

### Use existing package to fit model and confirm results

%% Cell type:code id: tags:

``` python
# function to apply simple linear regression
def func(x, b, a):
    return a + b*x

# curve fit the data
fitted_params, pcov = curve_fit(func, x, y)

model_predictions = func(x, *fitted_params)

absError = model_predictions - y

SE = np.square(absError) # squared errors
MSE = np.mean(SE) # mean squared errors
RMSE = np.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (np.var(absError) / np.var(y))
print('RMSE =', RMSE)
print('Rsquared=', Rsquared)
```

%% Cell type:markdown id: tags:

### Use another existing package to fit model and confirm results

%% Cell type:code id: tags:

``` python
slope, intercept, r_value, p_value, std_err = linregress(x, y)
print("slope: %f    intercept: %f" % (slope, intercept))
print("p-value: ", p_value/2)
#To get coefficient of determination (R-squared):
print("Rsquared: %f" % r_value**2)


#Plot the data along with the fitted line:
plt.figure(figsize=(6,6), dpi=100)
plt.plot(x, y, 'o', label='original data')
plt.plot(x, intercept + slope*x, 'r', label='fitted line')
plt.xlabel('Distance in KM') # X axis data label
plt.ylabel('Delivery Time in Days') # Y axis data label
plt.legend()
plt.show()
```

%% Cell type:markdown id: tags:

## Multiple Linear Regression

1. We would like to investigate the pasture rent structure with respect to the grass variety over
   various districts in Minnesota. The data set “pasture-data.csv” provides 67 rows of data.
2. The data columns include:
    - I, index
    - X1, the rent per arable acre (dollars)
    - X2, the number of milk cows per square mile
    - X3, the difference between pasturage and arable land
    - Y, the rental price per acre for this variety of grass (dollars)
3. We seek a model of the form:
      Y = B0 + B1 * X1 + B2 * X2 + B3 * X3

%% Cell type:code id: tags:

``` python
# load the data
data = pd.read_csv("pasture-data.csv")
```

%% Cell type:code id: tags:

``` python
data
```

%% Cell type:code id: tags:

``` python
# get x data .. only features, no index and outcome variables
X = data.iloc[:,1:-1]

# this is the outcome variable
y = data.iloc[:,-1]

# Create the train and test data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5)


# Create linear regression object
regr = LinearRegression()

# Train the model using the training sets
regr.fit(X_train, y_train)

# Make predictions using the testing set
y_pred = regr.predict(X_test)

# Compute residuals
residuals = (y_test - y_pred)
```

%% Cell type:markdown id: tags:

### Plot residuals

%% Cell type:code id: tags:

``` python
# Plot residuals
plt.scatter(range(len(residuals)),residuals)
plt.axhline(y=0, color='r', linestyle='-')
plt.show()
```

%% Cell type:code id: tags:

``` python
# Plot residuals
plt.hist(residuals)
plt.show()
```

%% Cell type:markdown id: tags:

### Obtain details of significance tests and CIs

%% Cell type:code id: tags:

``` python
mod = sm.OLS(y,X)
model = mod.fit()
model.summary2()
```

%% Cell type:code id: tags:

``` python
```
Loading