# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 2
#     language: python
#     name: python2
# ---

# ## LINEAR REGRESSION - PART 1 (ONE VARIABLE)
#
# One variable linear regression is the simplest form of machine learning we can think of. The purpose of this practice is to learn the basics of different components about linear regression. 
# #### DATASET http://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/

# #### Wisconsin Breast Cancer Database
# Samples arrive periodically as Dr. Wolberg reports his clinical cases.
#
# 1. Number of Instances: 699 (as of 15 July 1992)
#
# 2. Number of Attributes: 10 plus the class attribute
#
# 3. Attribute Information: (class attribute has been moved to last column)
#    Attribute                     Domain
#    1. Sample code number            id number
#    2. Clump Thickness               1 - 10
#    3. Uniformity of Cell Size       1 - 10
#    4. Uniformity of Cell Shape      1 - 10
#    5. Marginal Adhesion             1 - 10
#    6. Single Epithelial Cell Size   1 - 10
#    7. Bare Nuclei                   1 - 10
#    8. Bland Chromatin               1 - 10
#    9. Normal Nucleoli               1 - 10
#   10. Mitoses                       1 - 10
#   11. Class:                        (2 for benign, 4 for malignant)
#
# 4. Missing attribute values: 16
#
#    There are 16 instances in Groups 1 to 6 that contain a single missing 
#    (i.e., unavailable) attribute value, now denoted by "?".  
#
# 5. Class distribution:
#  
#    Benign: 458 (65.5%)
#    Malignant: 241 (34.5%)

# +
import numpy as np
import pandas as pd

# Read input
df = pd.read_csv("data/breast-cancer-wisconsin.data.csv")
df.head(3)

# +
# Preprocess the data so we can use it for regression. Specially the class values.
# Access function/method docstrings in jupyter via '?'. Ex: preprocessing.LabelEncoder?
from sklearn import preprocessing

encoder = preprocessing.LabelEncoder()

for col in df.columns:
    df[col] = encoder.fit_transform(df[col])
    
df.head(5)


# -

# #### Predict 
# This method takes x, w, b and returns the value of y. <br/>
# \begin{equation}
# y_i = w*x_i+b
# \end{equation}

def predict(x, w, b):
    return x*w + b


# #### Cost_function/loss_function/risk_function:
# The MSE is a measure of the quality of an estimator—it is always non-negative, and values closer to zero are better.
# We are using mean squared error. 
# \begin{equation}
# MSE = \frac{1}{n} \sum_{i=1}^n (y_i - (wx_i+b))^2
# \end{equation}

def cost_function(x, y, w, b):
    n = len(x)
    total_error = 0.0
    for i in range(n):
        total_error += (y[i] - (w*x[i] + b))**2
    return total_error / n


# #### Gradient descent
# To minimize MSE we use Gradient Descent to calculate the gradient of our cost function.
# There are two parameters (coefficients) in our cost function we can control: weight $w$ and bias $b$. Since we need to consider the impact each one has on the final prediction, we use partial derivatives.
# Recall the cost function:
# \begin{equation}
# f(w,b)= \frac{1}{n} \sum_{i=1}^{n} (y_i - (wx_i+b))^2
# \end{equation}
# First lets find the partial derivative for $w$:
# \begin{equation}
# \frac{\partial f}{\partial w} = \frac{1}{n} \sum_{i=1}^n -2x_i(y_i - (wx_i+b))
# \end{equation}
# And the partial derivative for $b$:
# \begin{equation}
# \frac{\partial f}{\partial b} = \frac{1}{n} \sum_{i=1}^n -2(y_i - (wx_i+b))
# \end{equation}
#

def update_weights(x, y, w, b, learning_rate):
    weight_deriv = 0
    bias_deriv = 0
    n = len(x)

    for i in range(n):
        # Calculate partial derivatives
        # -2x(y - (mx + b))
        weight_deriv += -2*x[i] * (y[i] - (w*x[i] + b))

        # -2(y - (mx + b))
        bias_deriv += -2*(y[i] - (w*x[i] + b))

    # We subtract because the derivatives point in direction of steepest ascent
    w -= (weight_deriv / n) * learning_rate
    b -= (bias_deriv / n) * learning_rate

    return w, b


# #### Train loop
#
# Train the model. Things to cover:
# 1. Parameter vs hyperparameter
# 2. Learning rate: possibly l1 and l2 regularization
# 3. Number of iterrations also popularly known as epochs

def train_model(x, y, w, b, learning_rate, epochs):
    cost_history = []

    for i in range(epochs):
        w,b = update_weights(x, y, w, b, learning_rate)

        #Calculate cost for auditing purposes
        cost = cost_function(x, y, w, b)
        cost_history.append(cost)

        # Log Progress
        if (i+1) % 20 == 0:
            print("Epochs: ", str(i+1), " cost: ",  str(cost))

    return w, b, cost_history


# #### Run the training
# Pick a parameter! Possibly mitoses as it's the worst one. Make them focus on the weight and bias values we get!
#

x = df['mitoses']
y = df['class']
w, b, cost_history = train_model(x, y, 0, 0, 0.02, 500)
print(w, b)


# #### Predictive separation index (PSI) as score
# We can use Predictive Separation Index (PSI), to use as the strength of a predictor. The equation is:
# \begin{equation}
# PSI ( x ) = [ \textrm{mean } ( \hat{ p } \textrm{ given } x ) \textrm{ when } y = 1 ] - [ \textrm{mean } ( \hat{ p } \textrm{ given } x ) \textrm{ when } y = 0 ] \, .
# \end{equation}
# We want PSI(x) very close to 1 as the first term should be close to 1 and the second term should be close to 0.

def get_score(x, y, w, b):
    preds_0 = []
    preds_1 = []
    for i in range(len(x)):
        p = predict(x[i], w, b)
        if y[i] == 0:
            preds_0.append(p)
        else:
            preds_1.append(p)
    if len(preds_0) != 0:
        score = (sum(preds_1) / len(preds_1) - sum(preds_0) / len(preds_0))
    else:
        score = (sum(preds_1) /  len(preds_1) - 0)
    return preds_0, preds_1, score


preds_0, preds_1, score = get_score(x, y, w, b)
print("PSI: ", score)

# ### The p distribution plot
#
# A better way to look at our prediction is to look at the distribution of $\hat{p}$ we get for each of the classes. For $y_i = 0$ the $\hat{p}$ should be close to 0 and $\hat{p}$ should be close to 1 when $y_i = 1$. This is a better way of seeing the strength of a predictor.

# +
# %matplotlib inline  
# The above is required to display matplotlib in jupyter
import matplotlib.pyplot as plt


n, bins, patches = plt.hist(preds_0, bins=100, normed=1, cumulative=0)
plt.title('Predictive distribution for class y=0')
plt.show()
n, bins, patches = plt.hist(preds_1, bins=100, normed=1, cumulative=0)
plt.title('Predictive distribution for class y=1')
plt.show()
# -

# ## Use a library to do the same thing! - SKLEARN
# If we use sklearn to do the same thing, we will get the same exact values! Voilà! Things to make sure:
# 1) Sklearn linear regression parameters: fit_intercept, normalize, copy_X, n_jobs
#     a) fit: Fit linear model. [Our train loop]
#     b) predict: Predict using the linear model. [Our predict method]
#     c) score: Returns the coefficient of determination $R^2$ of the prediction. <br/>
#     
# ##### The R squared score:
# R squared measures the proportion of overall variation of y variable that have been explained by the prediction. The equation for R squared is:
#  \begin{align*}
#  R^2 = 1 - \frac{\sum_{i=1}^n (y_i-\hat{y_i})^2}{\sum_{i=1}^n (y_i-\bar{y})^2}
#  \end{align*}
#  Here, $y_i$ is the label, $\hat{y_i}$ is the prediction and $\bar{y}$ is the average label. The error is normalized by the deviance and we get how well the regression no matter the deviation in the dataset. For a dichotomous classification this $R^2$ is same as predictive separation index but for continuous outcome it will be different.
#

x = df[['mitoses']]
y = df['class']
from sklearn import datasets, linear_model
lm_model = linear_model.LinearRegression(fit_intercept=True, normalize=True, copy_X=True, n_jobs=1)
lm_model.fit(x, y)
print("sklearn, w, b, score: ", round(lm_model.coef_[0], 6), round(lm_model.intercept_, 6), round(lm_model.score(x,y), 6))
print("manualr, w, b, score: ", round(w, 6) , round(b, 6), round(score, 6))
# So, we have learned the same co-efficents (or, almost the same) through a library.

# #### Variable selection
# Now loop through each of the feature/variable present in the dataset (except id and class) and fit the model to get a score. Sort features based on their score/R-squared-value the top variables are best predictors. <br/>
#
# A better approach would be: <br/>
# 1) Make sure all the features are in same scale. <br/>
# 2) Use better algorithm like step wise criterion function (ex. bayesian information criterion) or xgBoost which is supposedly one of the best available algorithm for feature selection. <br/>
#
# xgBoost: http://xgboost.readthedocs.io/en/latest/python/python_intro.html
#
# But this works pretty well too! We'll see!

# +
features = list(df.keys())
features.remove('id')
features.remove('class')

feature_scores = []

for feature in features:
    x = df[[feature]]
    y = df['class']
    lm_model = linear_model.LinearRegression(fit_intercept=True, normalize=True, copy_X=True, n_jobs=1)
    lm_model.fit(x, y)
    feature_scores.append([feature, lm_model.score(x,y)])

feature_scores.sort(key=lambda x: x[1], reverse=True)
for feature_score in feature_scores:
    print(feature_score)
# -

# #### Pick the top two predictors
#
# Pick the top two predictors and fit the model. The score will be better this time.

x = df[['uniformity-of-cell-shape','uniformity-of-cell-size']]
y = df['class']
lm_model = linear_model.LinearRegression(fit_intercept=True, normalize=True, copy_X=True, n_jobs=1)
lm_model.fit(x, y)
print("Score: ", lm_model.score(x,y))

# #### Pick a cutoff
#
# Now look at the distribution again and pick a cutoff. Everything below the cutoff will be 0 and above will be 1. 

# +
input_0 = df.loc[df['class'] == 0]
y_0 = input_0['class']
x_0 = input_0[['uniformity-of-cell-shape','uniformity-of-cell-size']]

input_1 = df.loc[df['class'] == 1]
y_1 = input_1['class']
x_1 = input_1[['uniformity-of-cell-shape','uniformity-of-cell-size']]
preds_0 = lm_model.predict(x_0)
preds_1 = lm_model.predict(x_1)
n, bins, patches = plt.hist(preds_0, bins=50, normed=1, cumulative=0)
plt.show()
n, bins, patches = plt.hist(preds_1, bins=50, normed=1, cumulative=0)
plt.show()
# -

# #### Accuracy
# Looks like 0.2 is a good cutoff? Let's pick that and look at the accuracy of the model

y_pred = lm_model.predict(x)
y_pred = [1 if p > 0.2 else 0 for p in y_pred]
from sklearn.metrics import accuracy_score
accuracy_score(y, y_pred)

# ## Splitting the dataset
# But we can't trust these results! We trained and tested on the same data. This may overfit the data and give us an illusion of better fit. Let's split the dataset in 80-20 for training-testing and fit and test the model again.

# +
# split into train and test data
from sklearn.model_selection import train_test_split

(train,test) = train_test_split(df, test_size=0.2)
train_output = train['class']
train_input = train[['uniformity-of-cell-shape','uniformity-of-cell-size']]

test_output = test['class']
test_input = test[['uniformity-of-cell-shape','uniformity-of-cell-size']]
# -

lm_model.fit(train_input, train_output)
print("Score: ", lm_model.score(train_input, train_output))

# +
input_0 = train.loc[train['class'] == 0]
y_0 = input_0['class']
x_0 = input_0[['uniformity-of-cell-shape','uniformity-of-cell-size']]

input_1 = train.loc[train['class'] == 1]
y_1 = input_1['class']
x_1 = input_1[['uniformity-of-cell-shape','uniformity-of-cell-size']]
preds_0 = lm_model.predict(x_0)
preds_1 = lm_model.predict(x_1)
n, bins, patches = plt.hist(preds_0, bins=50, normed=1, cumulative=0)
plt.show()
n, bins, patches = plt.hist(preds_1, bins=50, normed=1, cumulative=0)
plt.show()
# -

y_pred = lm_model.predict(test_input)
y_pred = [1 if p > 0.5 else 0 for p in y_pred]
from sklearn.metrics import accuracy_score
accuracy_score(test_output, y_pred)

# ## Break the accuracy myth
# We see that we are getting decent accuracy. But accuracy in a model doesn't mean the model is useful. There are many other statistical tests that can help us understand the effectiveness of a model. Sensitivity, specificity are widely used. We are going to look at true positive rate and false positive rates as simplest statistical significance values of the model.

from sklearn.metrics import confusion_matrix
print(confusion_matrix(test_output, y_pred))
tn, fp, fn, tp = confusion_matrix(test_output, y_pred).ravel()
fpr = fp / (fp+tn)
fnr = fn / (tp+fn)
print("False positive rate: (predicting malignant while benign)", fpr)
print("False negative rate: (predicting benign while malignant)", fnr)

# ## Which case is the worst?
# Think about this model in action. It has a very low false positive rate but a high false negative rate. False negative results in real world is very dangerous in our case. On the other hand, false positive results will cause a lot of chaos. But, if we have to make a decision, which is the best one?

# +
y_pred = lm_model.predict(test_input)
y_pred = [1 if p > 0.1 else 0 for p in y_pred]

print("Accuracy", accuracy_score(test_output, y_pred))

print(confusion_matrix(test_output, y_pred))
tn, fp, fn, tp = confusion_matrix(test_output, y_pred).ravel()
fpr = fp / (fp+tn)
fnr = fn / (tp+fn)
print("False positive rate: (predicting malignant while benign)", fpr)
print("False negative rate: (predicting benign while malignant)", fnr)

StarcoderdataJupyterScriptsDedupFiltered