# Import libraries necessary for this project
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
import seaborn as sns
import scipy.stats as stats

#Display for notebooks
%matplotlib inline

#Load the Boston housing dataset
data = pd.read_csv ( 'boston_data_1.csv' )
#Note that 'medv' been multiplicatively scaled to account for 35 years of market inflation
prices = data [ 'medv' ]
features = data.drop ( 'medv', axis = 1 )

#Obtain general info on dataframe
data.info()

#Look at first few rows of dataframe
data.head()

#Histographs for each variable
data.hist ( sharex = False, sharey = False, xlabelsize = 1, ylabelsize = 1 )
plt.show()

# ********************* *************** #
#   Learn more on the variable "crim"
# ************************************* #
# what are the basic statistics
data["crim"].describe()

def eda_numeric_plt(variable, y_max=-99, y_label=''):
    fig = plt.figure(figsize=(12,4))
    ax = fig.add_subplot(1,2,1)
    ax.hist(data[variable],50, color='lightblue', edgecolor='k')
    ax.set_xlabel(variable, fontsize=14)
    ax.set_ylabel("Counts", fontsize=14)
    ax.tick_params(labelsize=14)
    ax.set_title('Distribution of {}'.format(variable), fontsize=14)

    ax = fig.add_subplot(1,2,2)
    ax.boxplot(data[variable], widths=0.3)
    # if no value is provided, we skip next line
    if y_max != -99:
        ax.set_ylim([-1, y_max])
    ax.set_ylabel(y_label, fontsize=14)
    ax.tick_params(labelsize=14)
    ax.set_title('Boxplot', fontsize=16)
    plt.rcParams['axes.linewidth'] = 1
    plt.setp(ax.spines.values(), color='k')


    plt.show()
    plt.close()

eda_numeric_plt(variable='crim', y_max=20, y_label='Per Capita Crime Rate')

ia = np.where(data.crim>20)[0]
data.iloc[ia]

# ********************* *************** #
#   "zn": proportion of residential 
# land zoned for lots over 25,000 sq.ft.
# ************************************* #
data["zn"].describe()

eda_numeric_plt(variable='zn',y_max=50, y_label='Proportion of residential')

# ********************* *************** #
#   "indus": proportion of non-retail 
# . business acres per town
# ************************************* #
data["indus"].describe()

eda_numeric_plt(variable='indus',y_label='Proportion of non-retail business acres per town')

# ********************* *************** #
#   "chas": Charles River dummy variable (= 1 if tract bounds river; 0 otherwise).
# ************************************* #
# find out how many 0s and 1s
ia=np.where(data['chas']== 1)
print("counts of \'chas==1\' are: ", len(ia[0]))

# ********************* *************** #
#   "nox": nitrogen oxides concentration (parts per 10 million).
# ************************************* #
data['nox'].describe()

eda_numeric_plt(variable='nox',y_label='nitrogen oxides concentration')

# ********************* *************** #
#   "rm": average number of rooms per dwelling.
# ************************************* #
data['rm'].describe()

eda_numeric_plt(variable='rm',y_label='average number of rooms per dwelling.')

# ********************* *************** #
#   "age": proportion of owner-occupied units built prior to 1940.
# ************************************* #
data['age'].describe()

eda_numeric_plt(variable='age',y_label='proportion of owner-occupied units built prior to 1940')

# ********************* *************** #
#   "dis": weighted mean of distances to five Boston employment centres.
# ************************************* # 
data['dis'].describe()

eda_numeric_plt(variable='dis',y_label='weighted mean of distances')

# ********************* *************** #
#   "rad": index of accessibility to radial highways.
# ************************************* # 
data['rad'].describe()

eda_numeric_plt(variable='rad',y_label='index of accessibility to radial highways',y_max=35)

# ********************* *************** #
#   "tax": full-value property-tax rate per $10,000.
# ************************************* # 
data['tax'].describe()

eda_numeric_plt(variable='tax',y_label='full-value property-tax rate per $10,000')

# ********************* *************** #
#   "ptratio": pupil-teacher ratio by town.
# ************************************* # 
data['ptratio'].describe()

eda_numeric_plt(variable='ptratio',y_label='pupil-teacher ratio by town')

# ********************* *************** #
#   "black": 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town.
# ************************************* # 
data['black'].describe()

eda_numeric_plt(variable='black',y_label='1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town',y_max=500)

# ********************* *************** #
#   "lstat": lower status of the population (percent).
# ************************************* # 
data["lstat"].describe()

eda_numeric_plt(variable='lstat',y_label='lower status of the population (percent)')

data["medv"].describe()

eda_numeric_plt(variable='medv',y_label='Median Value of owner occupied homes')

def correlations(data, y, xs):
    rs = []
    rhos = []
    for x in xs:
        r = stats.pearsonr(data[y], data[x])[0]
        rs.append(r)
        rho = stats.spearmanr(data[y], data[x])[0]
        rhos.append(rho)
    return pd.DataFrame({"feature": xs, "r": rs, "rho": rhos})

variables = list(data.columns) 

for i, var_i in enumerate(variables):
    if i< len(variables)-1:
        print("correlations with %s" % var_i)
        print(correlations(data, var_i, variables[i+1:]))
    print("\n\n")

#Heat map
plt.figure ( figsize = ( 20, 10 ) )
sns.heatmap ( data.corr().abs(),  annot = True )

import warnings
warnings.filterwarnings('ignore')
from sklearn import preprocessing
#Let's scale the columns before plotting them against medv
min_max_scaler = preprocessing.MinMaxScaler()
column_sels = [ 'lsat', 'indus', 'nox', 'ptratio', 'rm', 'tax', 'dis', 'age' ]
x = data.loc[:,column_sels]
y = data[ 'medv' ]
x = pd.DataFrame ( data = min_max_scaler.fit_transform(x), columns = column_sels )
fig, axs = plt.subplots ( ncols = 4, nrows = 2, figsize = ( 20, 10 ) )
index = 0
axs = axs.flatten()
for i, k in enumerate ( column_sels ):
    sns.regplot ( y = y, x = x [ k ], ax = axs[ i ] )
plt.tight_layout ( pad = 0.4, w_pad = 0.5, h_pad = 5.0 )

# data transformation

# divide rad into two categories < 15 & > 15
data['rad_lt_15'] = data.rad.apply(lambda x: 1 if x<15 else 0)

# convert lstat to a log scale
data['lstat_log'] = np.log(data.lstat)

data.head()

import models

# build the second linear model
model = " medv ~ rm + rad_lt_15 + lstat_log"
result = models.bootstrap_linear_regression(model, data=data)
models.describe_bootstrap_lr(result)

import random
import patsy

def chunk(xs, n):
    k, m = divmod(len(xs), n)
    return [xs[i * k + min(i, m):(i + 1) * k + min(i + 1, m)] for i in range(n)]

def cross_validation(algorithm, formula, data, evaluate, fold_count=10, repetitions=3):
    indices = list(range(len( data)))
    metrics = []
    for _ in range(repetitions):
        random.shuffle(indices)
        folds = chunk(indices, fold_count)
        for fold in folds:
            test_data = data.iloc[fold]
            train_indices = [idx not in fold for idx in indices]
            train_data = data.iloc[train_indices]
            result = algorithm(formula, data=train_data)
            model = result["model"]
            y, X = patsy.dmatrices(formula, test_data, return_type="matrix")
            # y = np.ravel( y) # might need for logistic regression
            results = models.summarize(formula, X, y, model)
            metric = evaluate(results)
            metrics.append(metric)
    return metrics


lr_result_1 = cross_validation(models.linear_regression, model, data, lambda r: (r["sigma"], r["r_squared"]))
print(r"95% CI for R^2:", stats.mstats.mquantiles([r[1] for r in lr_result_1], [0.025, 0.975]))

# data transformation
data["log_dis"] = data["dis"].apply(np.log)
data["log_crim"] = data["crim"].apply(np.log)
data["sqrt_lstat"] = data["lstat"].apply(np.sqrt)

model2 = "medv ~ rm + sqrt_lstat + ptratio + log_dis + log_crim" 
result2 = models.bootstrap_linear_regression(model2, data=data)
models.describe_bootstrap_lr(result2)

lr_result_2 = cross_validation(models.linear_regression, model2, data, lambda r: (r["sigma"], r["r_squared"]))
print(r"95% CI for R^2:", stats.mstats.mquantiles([r[1] for r in lr_result_2], [0.025, 0.975]))

model3 = "medv ~ black + tax + dis + lstat_log + ptratio + rm"
results = models.linear_regression(model3, data)
models.simple_describe_lr(results)

lr_result_3 = cross_validation(models.linear_regression, model3, data, lambda r: (r["sigma"], r["r_squared"]))
print(r"95% CI for R^2:", stats.mstats.mquantiles([r[1] for r in lr_result_3], [0.025, 0.975]))

# build a cross validation curve for decision tree 
def decisiontree_cross_validation(formula, builder, data, fold_count=10, repetitions=3):
    indices = list(range(len( data)))
    metrics = {"train": [], "test": []}
    for _ in range(repetitions):
        random.shuffle(indices)
        folds = chunk(indices, fold_count)
        for fold in folds:
            test_data = data.iloc[fold]
            train_indices = [idx not in fold for idx in indices]
            train_data = data.iloc[train_indices]
            # y, X for training data
            y, X = patsy.dmatrices(formula, train_data, return_type="matrix")
            model = builder.fit(X, y)
            y_hat = model.predict(X)
            training_r_squared = (stats.pearsonr(y[:,0], y_hat)[0])**2
            metrics["train"].append(training_r_squared)
            # y, X for training data
            y, X = patsy.dmatrices(formula, test_data, return_type="matrix")
            y_hat = model.predict(X)
            test_r_squared = (stats.pearsonr(y[:,0], y_hat)[0])**2
            metrics["test"].append(test_r_squared)
    return metrics

from sklearn.tree import DecisionTreeRegressor

# test a range of tree depths and create a validation curve
test_curve = []
train_curve = []
for k in range(1, 21):
    builder = DecisionTreeRegressor(max_depth=k)
    results = decisiontree_cross_validation(model, builder, data)
    test_curve.append(stats.mstats.mquantiles(results['test'], [0.025, 0.5, 0.975]))
    train_curve.append(stats.mstats.mquantiles(results['train'], [0.025, 0.5, 0.975]))

# visualize the results
xs = list(range(1, 21))
figure = plt.figure(figsize=(10, 6)) # first element is width, second is height.
axes = figure.add_subplot(1, 1, 1)

test_lower, test_mid, test_upper = zip(*test_curve)
train_lower, train_mid, train_upper = zip(*train_curve)

axes.fill_between(xs, test_lower, test_upper, alpha=0.25, color="firebrick", label="test")
axes.plot(xs, test_mid, color="firebrick")

axes.fill_between(xs, train_lower, train_upper, alpha=0.25, color="cornflowerblue", label="train")
axes.plot(xs, train_mid, color="cornflowerblue")

axes.set_xticks(xs)
axes.set_xticklabels([str(x) for x in xs])

axes.legend()
axes.set_title("Validation Curves for Decision Tree")
axes.set_xlabel("Depth of Decision Tree")
axes.set_ylabel("$R^2$")

plt.show()
plt.close()

from collections import defaultdict

def data_collection():
    result = dict()
    result[ "train"] = defaultdict( list)
    result[ "test"] = defaultdict( list)
    return result

def decisiontree_learning_curves(builder, formula, data, fold_count=10, repetitions=3, increment=1):
    indices = list(range(len( data)))
    results = data_collection()
    for _ in range(repetitions):
        random.shuffle(indices)
        folds = chunk(indices, fold_count)
        for fold in folds:
            test_data = data.iloc[ fold]
            train_indices = [idx for idx in indices if idx not in fold]
            train_data = data.iloc[train_indices]
            for i in list(range(increment, 100, increment)) + [100]: # ensures 100% is always picked.
                # the indices are already shuffled so we only need to take ever increasing chunks
                train_chunk_size = int( np.ceil((i/100)*len( train_indices)))
                train_data_chunk = data.iloc[train_indices[0:train_chunk_size]]
                # we calculate the model
                y, X = patsy.dmatrices(formula, train_data_chunk, return_type="matrix")
                model = builder.fit(X, y)
                y_hat = model.predict(X)
                training_r_squared = (stats.pearsonr(y[:,0], y_hat)[0])**2
                results["train"][i].append(training_r_squared)
                
                # y, X for test data
                y, X = patsy.dmatrices(formula, test_data, return_type="matrix")
                y_hat = model.predict(X)
                test_r_squared = (stats.pearsonr(y[:,0], y_hat)[0])**2
                results["test"][i].append(test_r_squared)
    # process results
    statistics = {}
    for k, v in results["train"].items():
        statistics[ k] = (np.mean(v), np.std(v))
    results["train"] = statistics
    statistics = {}
    for k, v in results["test"].items():
        statistics[ k] = (np.mean(v), np.std(v))
    results["test"] = statistics
    return results

# plot a learning curve
builder = DecisionTreeRegressor(max_depth=13)
result = decisiontree_learning_curves(builder, model, data)
def results_to_curves( curve, results):
    all_statistics = results[ curve]
    keys = list( all_statistics.keys())
    keys.sort()
    mean = []
    upper = []
    lower = []
    for k in keys:
        m, s = all_statistics[ k]
        mean.append( m)
        upper.append( m + 2 * s)
        lower.append( m - 2 * s)
    return keys, lower, mean, upper

def plot_learning_curves( results, metric, zoom=False):
    figure = plt.figure(figsize=(10,6))

    axes = figure.add_subplot(1, 1, 1)

    xs, train_lower, train_mean, train_upper = results_to_curves( "train", results)
    _, test_lower, test_mean, test_upper = results_to_curves( "test", results)

    axes.plot( xs, train_mean, color="steelblue")
    axes.fill_between( xs, train_upper, train_lower, color="steelblue", alpha=0.25, label="train")
    axes.plot( xs, test_mean, color="firebrick")
    axes.fill_between( xs, test_upper, test_lower, color="firebrick", alpha=0.25, label="test")
    axes.legend()
    axes.set_xlabel( "training set (%)")
    axes.set_ylabel( metric)
    axes.set_title("Learning Curves")

    if zoom:
        y_lower = int( 0.9 * np.amin([train_lower[-1], test_lower[-1]]))
        y_upper = int( 1.1 * np.amax([train_upper[-1], test_upper[-1]]))
        axes.set_ylim((y_lower, y_upper))

    plt.show()
    plt.close()

plot_learning_curves(result, "$R^2$")

DT_results = decisiontree_cross_validation(model, builder, data)
print(r"95% CI for R^2:", stats.mstats.mquantiles(DT_results['test'], [0.025, 0.975]))

from sklearn.preprocessing import StandardScaler
import sklearn.neighbors as neighbors
import sklearn.linear_model as linear
import random
import patsy

columns = ['medv', 'black' , 'tax' , 'dis' , 'lstat_log' , 'ptratio' , 'rm'] 
scaler = StandardScaler()
data[columns] = scaler.fit_transform(data[columns])

def chunk(xs, n):
    k, m = divmod(len(xs), n)
    return [xs[i * k + min(i, m):(i + 1) * k + min(i + 1, m)] for i in range(n)]

def knn_cross_validation(formula, builder, data, fold_count=10, repetitions=3):
    indices = list(range(len( data)))
    metrics = {"train": [], "test": []}
    for _ in range(repetitions):
        random.shuffle(indices)
        folds = chunk(indices, fold_count)
        for fold in folds:
            test_data = data.iloc[fold]
            train_indices = [idx not in fold for idx in indices]
            train_data = data.iloc[train_indices]
            # y, X for training data
            y, X = patsy.dmatrices(formula, train_data, return_type="matrix")
            model = builder.fit(X, y)
            y_hat = model.predict(X)
            training_r_squared = (stats.pearsonr(y, y_hat)[0][0])**2
            metrics["train"].append(training_r_squared)
            # y, X for training data
            y, X = patsy.dmatrices(formula, test_data, return_type="matrix")
            y_hat = model.predict(X)
            test_r_squared = (stats.pearsonr(y, y_hat)[0][0])**2
            metrics["test"].append(test_r_squared)
    return metrics

def data_collection(test_curve, train_curve):
    xs = list(range(1, 21))
    figure = plt.figure(figsize=(10, 6)) # first element is width, second is height.
    axes = figure.add_subplot(1, 1, 1)

    test_lower, test_mid, test_upper = zip(*test_curve)
    train_lower, train_mid, train_upper = zip(*train_curve)

    axes.fill_between(xs, test_lower, test_upper, alpha=0.25, color="firebrick", label="test")
    axes.plot(xs, test_mid, color="firebrick")

    axes.fill_between(xs, train_lower, train_upper, alpha=0.25, color="cornflowerblue", label="train")
    axes.plot(xs, train_mid, color="cornflowerblue")

    axes.set_xticks(xs)
    axes.set_xticklabels([str(x) for x in xs])

    axes.legend()
    axes.set_title("Validation Curves for kNN")
    axes.set_xlabel("k")
    axes.set_ylabel("$R^2$")

    plt.show()
    plt.close()
from collections import defaultdict
def create_data_collection():
    result = dict()
    result["train"] = defaultdict(list)
    result["test"]  = defaultdict(list)
    return result
def knn_learning_curves(builder, formula, data, fold_count=10, repetitions=3, increment=1):
    indices = list(range(len( data)))
    results = create_data_collection()
    for _ in range(repetitions):
        random.shuffle(indices)
        folds = chunk(indices, fold_count)
        for fold in folds:
            test_data = data.iloc[ fold]
            train_indices = [idx for idx in indices if idx not in fold]
            train_data = data.iloc[train_indices]
            for i in list(range(increment, 100, increment)) + [100]: # ensures 100% is always picked.
                # the indices are already shuffled so we only need to take ever increasing chunks
                train_chunk_size = int( np.ceil((i/100)*len( train_indices)))
                train_data_chunk = data.iloc[train_indices[0:train_chunk_size]]
                # we calculate the model
                y, X = patsy.dmatrices(formula, train_data_chunk, return_type="matrix")
                model = builder.fit(X, y)
                y_hat = model.predict(X)
                training_r_squared = (stats.pearsonr(y, y_hat)[0][0])**2
                results["train"][i].append(training_r_squared)
                
                # y, X for test data
                y, X = patsy.dmatrices(formula, test_data, return_type="matrix")
                y_hat = model.predict(X)
                test_r_squared = (stats.pearsonr(y, y_hat)[0][0])**2
                results["test"][i].append(test_r_squared)
    # process results
    # Rely on the CLT...
    statistics = {}
    for k, v in results["train"].items():
        statistics[ k] = (np.mean(v), np.std(v))
    results["train"] = statistics
    statistics = {}
    for k, v in results["test"].items():
        statistics[ k] = (np.mean(v), np.std(v))
    results["test"] = statistics
    return results
#
def results_to_curves( curve, results):
    all_statistics = results[ curve]
    keys = list( all_statistics.keys())
    keys.sort()
    mean = []
    upper = []
    lower = []
    for k in keys:
        m, s = all_statistics[ k]
        mean.append( m)
        upper.append( m + 2 * s)
        lower.append( m - 2 * s)
    return keys, lower, mean, upper
def plot_learning_curves( results, metric, zoom=False):
    figure = plt.figure(figsize=(10,6))

    axes = figure.add_subplot(1, 1, 1)

    xs, train_lower, train_mean, train_upper = results_to_curves( "train", results)
    _, test_lower, test_mean, test_upper = results_to_curves( "test", results)

    axes.plot( xs, train_mean, color="steelblue")
    axes.fill_between( xs, train_upper, train_lower, color="steelblue", alpha=0.25, label="train")
    axes.plot( xs, test_mean, color="firebrick")
    axes.fill_between( xs, test_upper, test_lower, color="firebrick", alpha=0.25, label="test")
    axes.legend()
    axes.set_xlabel( "training set (%)")
    axes.set_ylabel( metric)
    axes.set_title("Learning Curves")

    if zoom:
        y_lower = int( 0.9 * np.amin([train_lower[-1], test_lower[-1]]))
        y_upper = int( 1.1 * np.amax([train_upper[-1], test_upper[-1]]))
        axes.set_ylim((y_lower, y_upper))

    plt.show()
    plt.close()

formula = "medv ~ black + tax + dis + lstat_log + ptratio + rm"

test_curve = []
train_curve = []
for k in range(1, 21):
    builder = neighbors.KNeighborsRegressor(k)
    results = knn_cross_validation(formula, builder, data)
    test_curve.append(stats.mstats.mquantiles(results['test'], [0.025, 0.5, 0.975]))
    train_curve.append(stats.mstats.mquantiles(results['train'], [0.025, 0.5, 0.975]))
data_collection(test_curve, train_curve)

result = knn_learning_curves(neighbors.KNeighborsRegressor(1), formula, data)
plot_learning_curves(result, "$R^2$")

builder = neighbors.KNeighborsRegressor(1)
knn_results = knn_cross_validation(formula, builder, data)
print(r"95% CI for R^2:", stats.mstats.mquantiles(knn_results['test'], [0.025, 0.975]))

def bootstrap_sample( data, f, n=100):
    result = []
    m = len( data)
    for _ in range( n):
        sample = np.random.choice( data, len(data), replace=True)
        r = f( sample)
        result.append( r)
    return np.array( result)

lr1_r2 = [r[1] for r in lr_result_1]
lr2_r2 = [r[1] for r in lr_result_2]
lr3_r2 = [r[1] for r in lr_result_3]
dt_r2 = DT_results['test']
knn_bootstrap = bootstrap_sample(knn_results['test'], np.mean,30)

sns.set(style="whitegrid")
# visualize result
figure = plt.figure(figsize=(8, 8)) # first element is width, second is height.

axes = figure.add_subplot(2, 3, 1)
axes.hist(lr1_r2, density=True)
axes.set_ylabel( "Density")
axes.set_title( "Linear Model 1")

axes = figure.add_subplot(2,3, 2)
axes.hist(lr2_r2, density=True)
axes.set_ylabel( "Density")
axes.set_title( "Linear Model 2")

axes = figure.add_subplot(2,3, 3)
axes.hist( lr3_r2, density=True)
axes.set_ylabel( "Density")
axes.set_xlabel( "$R^2$")
axes.set_title( "Linear Model 3")

axes = figure.add_subplot(2,3, 4)
axes.hist( dt_r2, density=True)
axes.set_ylabel( "Density")
axes.set_xlabel( "$R^2$")
axes.set_title( "Decision Tree")

axes = figure.add_subplot(2,3, 5)
axes.hist( knn_bootstrap, density=True)
axes.set_ylabel( "Density")
axes.set_xlabel( "$R^2$")
axes.set_title( "k-nearest neighbors")


plt.show()
plt.close()

diff1 = np.array(lr2_r2) - np.array(lr1_r2)
print("P(LR2 > LR1)", np.mean(diff1 > 0))

diff2a = np.array(lr3_r2) - np.array(lr1_r2)
print("P(LR3 > LR1)", np.mean(diff2a > 0))
diff2b = np.array(lr3_r2) - np.array(lr2_r2)
print("P(LR3 > LR2)", np.mean(diff2b > 0))

diff3a = np.array(dt_r2) - np.array(lr1_r2)
print("P(DT > LR1)", np.mean(diff3a > 0))
diff3b = np.array(dt_r2) - np.array(lr2_r2)
print("P(DT > LR2)", np.mean(diff3b > 0))
diff3c = np.array(dt_r2) - np.array(lr3_r2)
print("P(DT > LR3)", np.mean(diff3c > 0))

diff3a = np.array(knn_bootstrap) - np.array(lr1_r2)
print("P(KNN > LR1)", np.mean(diff3a > 0))
diff3b = np.array(knn_bootstrap) - np.array(lr2_r2)
print("P(KNN > LR2)", np.mean(diff3b > 0))
diff3c = np.array(knn_bootstrap) - np.array(lr3_r2)
print("P(KNN > LR3)", np.mean(diff3c > 0))
diff4c = np.array(knn_bootstrap) - np.array(dt_r2)
print("P(KNN > dt_r2)", np.mean(diff4c > 0))