<jupyter_start><jupyter_text>Before you turn this problem in, make sure everything runs as expected. First, **restart the kernel** (in the menubar, select Kernel$\rightarrow$Restart) and then **run all cells** (in the menubar, select Cell$\rightarrow$Run All).Make sure you fill in any place that says `YOUR CODE HERE` or "YOUR ANSWER HERE", as well as your name and collaborators below:<jupyter_code>NAME = ""
ID = ""<jupyter_output><empty_output><jupyter_text>---  Lab 03: Logistic RegressionThus far, the problems we've encountered have been *regression* problems, in which the target $y \in \mathbb{R}$.Today we'll start experimenting with *classification* problems, beginning with *binary* classification problems, in which the target $y \in \{ 0, 1 \}$. BackgroundThe simplest approach to classification, applicable when the input feature vector $\mathbf{x} \in \mathbb{R}^n$, is a simple generalization of what wedo in linear regression. Recall that in linear regression, we assume that the target is drawn from a Gaussian distribution whose mean is a linear functionof $\mathbf{x}$:$$ y \sim {\cal N}(\theta^\top \mathbf{x}, \sigma^2) $$In logistic regression, similarly, we'll assume that the target is drawn from a Bernoulli distribution with parameter $p$ being the probability ofclass 1:$$ y \sim \text{Bernoulli}(p) $$That's fine, but how do we model the parameter $p$? How is it related to $\mathbf{x}$?In linear regression, we assume that the mean of the Gaussian is $\theta^\top \mathbf{x}$, i.e., a linear function of $\mathbf{x}$.In logistic regression, we'll assume that $p$ is a "squashed" linear function of $\mathbf{x}$, i.e.,$$ p = \text{sigmoid}(\theta^\top \mathbf{x}) = g(\theta^\top \mathbf{x}) = \frac{1}{1+e^{-\theta^\top \mathbf{x}}}. $$Later, when we introduce generalized linear models, we'll see why $p$ should take this form. For now, though, we can simply note that the selection makessense. Since $p$ is a discrete probability, $p$ is bounded by $0 \le p \le 1$. The sigmoid function $g(\cdot)$ conveniently obeys these bounds:<jupyter_code>import numpy as np
import matplotlib.pyplot as plt

def f(z):
    return 1 / (1 + np.exp( -z ))

z = np.arange(-5, 5, 0.1)
plt.plot(z, f(z), 'b-')
plt.xlabel('z')
plt.ylabel('g(z)')
plt.title('Logistic sigmoid function')
plt.show()<jupyter_output><empty_output><jupyter_text>We see that the sigmoid approaches 0 as its input approaches $-\infty$ and approaches 1 as its input approaches $+\infty$. If its input is 0, its value is 0.5.Again, this choice of function may seem strange at this point, but bear with it! We'll derive this function from a more general principle, the generalizedlinear model, later.OK then, we now understand that for logistic regression, the assumptions are:1. The *data* are pairs $(\textbf{x}, y) \in \mathbb{R}^n \times \{ 0, 1 \}$.1. The *hypothesis function* is $h_\theta(\textbf{x}) = \frac{1}{1+e^{-\theta^\top \mathbf{x}}}$.What else do we need... ? A cost function and an algorithm for minimizing that cost function! Cost function for logistic regressionYou can refer to the lecture notes to see the derivation, but for this lab, let's just skip to the chase.With the hypothesis $h_\theta(\mathbf{x})$ chosen as above, the log likelihood function $\ell(\theta)$ can be derived as$$ \ell(\theta) = \log {\cal L}(\theta) =  \sum_{i=1}^{m}y^{(i)}\log(h_{\theta}(\mathbf{x}^{(i)})) + (1 - y^{(i)})\log(1 - (h_{\theta}(\mathbf{x}^{(i)})) .$$Negating the log likelihood function to obtain a loss function, we have$$ J(\theta) = - \sum_{i=1}^m y^{(i)}\log h_\theta(\mathbf{x}^{(i)}) + (1-y^{(i)})\log(1-h_\theta(\textbf{x}^{(i)})) .$$There is no closed-form solution to this problem like there is in linear regression, so we have to use gradient descent to find $\theta$ minimizing $J(\theta)$.Luckily, the function *is* convex in $\theta$ so there is just a single global minimum, and gradient descent is guaranteed to get us there eventually if we takethe right step size.The *stochastic* gradient of $J$, for a single observed pair $(\mathbf{x}, y)$, turns out to be (see lecture notes)$$\nabla_J(\theta) = (h_\theta(\mathbf{x}) - y)\mathbf{x} .$$Give some thought as to whether following this gradient to increase the loss $J$ would make a worse classifier, and vice versa!Finally, we obtain the update rule for the $j$th iteration selecting training pattern $i$:$$ \theta^{(j+1)} \leftarrow \theta^{(j)} + \alpha(y^{(i)} - h_\theta(\textbf{x}^{(i)}))\textbf{x}^{(i)} .$$ Note that we can perform *batch gradient descent* simply by summing the single-pair gradient over the entire training set before taking a step,or *mini-batch gradient descent* by summing over a small subset of the data. Example dataset 1: student admissions dataThis example is from <NAME>'s machine learning course on Coursera.The data contain students' scores for two standardized tests and an admission decision (0 or 1).<jupyter_code># Load student admissions data. The data file does not contain headers,
# so we use hard coded indices for exam 1, exam2, and the admission decision.

data = np.loadtxt('ex2data1.txt',delimiter = ',')
exam1_data = data[:,0]
exam2_data = data[:,1]
X = np.array([exam1_data, exam2_data]).T
y = data[:,2]

# Output some sample data

print('Exam scores', X[0:5,:])
print('-----------------------------')
print('Admission decision', y[0:5])<jupyter_output><empty_output><jupyter_text>Let's plot the data...<jupyter_code># Plot the data

idx_0 = np.where(y == 0)
idx_1 = np.where(y == 1)

fig1 = plt.figure(figsize=(5, 5)) 
ax = plt.axes()
ax.set_aspect(aspect = 'equal', adjustable = 'box')
plt.title('Distribution')
plt.xlabel('Exam 1')
plt.ylabel('Exam 2')
plt.grid(axis='both', alpha=.25)
ax.scatter(exam1_data[idx_0], exam2_data[idx_0], s=50, c='r', marker='*', label='Not Admitted')
ax.scatter(exam1_data[idx_1], exam2_data[idx_1], s=50, c='b', marker='o', label='Admitted')
plt.show()<jupyter_output><empty_output><jupyter_text>Let's see if we can find good values for $\theta$ without normalizing the data.We will definitely want to split the data into train and test, however...<jupyter_code>import random

# As usual, we fix the seed to eliminate random differences between different runs

random.seed(12)

# Partion data into training and test datasets

m, n = X.shape
XX = np.insert(X, 0, 1, axis=1)
y = y.reshape(m, 1)
idx = np.arange(0, m)
random.shuffle(idx)
percent_train = .6
m_train = int(m * percent_train)
train_idx = idx[0:m_train]
test_idx = idx[m_train:]
X_train = XX[train_idx,:];
X_test = XX[test_idx,:];

y_train = y[train_idx];
y_test = y[test_idx];<jupyter_output><empty_output><jupyter_text>Important functions needed laterLet's put all of our important functions here...<jupyter_code>def sigmoid(z):   
    return 1 / (1 + np.exp(-z))

def h(X, theta):
    return sigmoid(X @ theta)

def grad_j(X, y, y_pred):
    return X.T @ (y - y_pred) / X.shape[0]
    
def j(theta, X, y):    
    y_pred = h(X, theta)
    error = (-y * np.log(y_pred)) - ((1 - y) * np.log(1 - y_pred))
    cost = sum(error) / X.shape[0]
    grad = grad_j(X, y, y_pred)
    return cost[0], grad<jupyter_output><empty_output><jupyter_text>Initialize thetaIn any iterative algorithm, we need an initial guess. Here we'll just use zeros for all parameters.<jupyter_code># Initialize our parameters, and use them to make some predictions

theta_initial = np.zeros((n+1, 1))

print('Initial theta:', theta_initial)
print('Initial predictions:', h(XX, theta_initial)[0:5,:])
print('Targets:', y[0:5,:])<jupyter_output><empty_output><jupyter_text>Training functionHere's a function to do batch training for `num_iters` iterations.<jupyter_code>def train(X, y, theta_initial, alpha, num_iters):
    theta = theta_initial
    j_history = []
    for i in range(num_iters):
        cost, grad = j(theta, X, y)
        theta = theta + alpha * grad
        j_history.append(cost)
    return theta, j_history<jupyter_output><empty_output><jupyter_text>Do the trainingHere we run the training function for a million batches!<jupyter_code># Train for 1000000 iterations on full training set

alpha = .0005
num_iters = 1000000
theta, j_history = train(X_train, y_train, theta_initial, alpha, num_iters)

print("Theta optimized:", theta)
print("Cost with optimized theta:", j_history[-1])<jupyter_output><empty_output><jupyter_text>Plot the loss curveNext let's plot the loss curve (loss as a function of iteration).<jupyter_code>plt.plot(j_history)
plt.xlabel("Iteration")
plt.ylabel("$J(\theta)$")
plt.title("Training cost over time with batch gradient descent (no normalization)")
plt.show()<jupyter_output><empty_output><jupyter_text>In-lab exercise from Example 1 (Total 35 points)That took a long time, right?We'll see if we can do better. We will try the following:1. Try increasing the learning rate $\alpha$ and starting with a better initial $\theta$. How much does it help?   - Try at least 2 learning rate $\alpha$ with 2 difference $\theta$ (4 experiments)   - Do not forget to plot the loss curve to compare your results2. Better yet, try *normalizing the data* and see if the training converges better. How did it go?    - Be sure to plot loss curves to compare the results with unnormalized and normalized data.3. Discuss the effects of normalization, learning rate, and initial $\theta$ in your report.Do this work in the following steps.  Exercise 1.1 (5 points)Fill in two different values for $\alpha$ and $\theta$.Use variable names `alpha1`, `alpha2`, `theta_initial1`, and `theta_initial2`.<jupyter_code># grade task: change 'None' value to number(s) or function
# YOUR CODE HERE
raise NotImplementedError()
# declare your alphas
# alpha1 = None
# alpha2 = None

# initialize thetas as you want
# theta_initial1 = None
# theta_initial2 = None

# define your num iterations
# num_iters = None
alpha_list = [alpha1, alpha2]
print('alpha 1:', alpha1)
print('alpha 2:', alpha2)

theta_initial_list = [theta_initial1, theta_initial2]
print('theta 1:', theta_initial_list[0])
print('theta 2:', theta_initial_list[1])

print('Use num iterations:', num_iters)

# Test function: Do not remove
assert alpha_list[0] is not None and alpha_list[1] is not None, "Alpha has not been filled"
chk1 = isinstance(alpha_list[0], (int, float))
chk2 = isinstance(alpha_list[1], (int, float))
assert chk1 and chk2, "Alpha must be number"
assert theta_initial_list[0] is not None and theta_initial_list[1] is not None, "initialized theta has not been filled"
chk1 = isinstance(theta_initial_list[0], (list,np.ndarray))
chk2 = isinstance(theta_initial_list[1], (list,np.ndarray))
assert chk1 and chk2, "Theta must be list"
chk1 = ((n+1, 1) == theta_initial_list[0].shape)
chk2 = ((n+1, 1) == theta_initial_list[1].shape)
assert chk1 and chk2, "Theta size are incorrect"
assert num_iters is not None and isinstance(num_iters, int), "num_iters must be integer"
print("success!")
# End Test function<jupyter_output><empty_output><jupyter_text>Exercise 1.2 (5 points)Fill in the code required to train your model on a particular $\alpha$ and $\theta$:<jupyter_code># grade task: change 'None, None' value to number(s) or function
j_history_list = []
theta_list = []
for alpha in alpha_list:
    for theta_initial in theta_initial_list:
        # YOUR CODE HERE
        raise NotImplementedError()
        # theta_i, j_history_i = None, None
        j_history_list.append(j_history_i)
        theta_list.append(theta_i)
# Test function: Do not remove
assert theta_list[0] is not None and j_history_list[0] is not None, "No values in theta_list or j_history_list"
chk1 = isinstance(theta_list[0], (list,np.ndarray))
chk2 = isinstance(j_history_list[0][0], (int, float))
assert chk1 and chk2, "Wrong type in theta_list or j_history_list"
print("success!")
# End Test function<jupyter_output><empty_output><jupyter_text>Exercise 1.3 (10 points)Write code to plot loss curves for each of the sequences in `j_history_list` from the previous exercise:<jupyter_code># YOUR CODE HERE
raise NotImplementedError()<jupyter_output><empty_output><jupyter_text>Exercise 1.4 (10 points)- Repeat your training, but **normalize** your data before training- Compare the results between normalized data and unnormalized data<jupyter_code># code here<jupyter_output><empty_output><jupyter_text>Exercise 1.5 (5 points)Discuss the effects of normalization, learning rate, and initial $\theta$ in your report. Write your discussion here.  The logistic regression decision boundaryNote that when $\theta^\top \textbf{x} = 0$, we have $h_\theta(\textbf{x}) = 0.5$. That is, we areequally unsure as to whether $\textbf{x}$ belongs to class 0 or class 1. The contour at which$h_\theta(\textbf{x}) = 0.5$ is called the classifier's *decision boundary*.We know that in the plane, the equation $$ax+by+c=0$$ is the general form of a 2D line. In our case, we have$$\theta_0 + \theta_1 x_1 + \theta_2 x_2 = 0$$ as our decision boundary, but clearly, this is just a 2D linein the plane. So when we plot $x_1$ against $x_2$, it is easy to plot the boundary line.<jupyter_code>def boundary_points(X, theta):
    v_orthogonal = np.array([[theta[1,0]],[theta[2,0]]])
    v_ortho_length = np.sqrt(v_orthogonal.T @ v_orthogonal)
    dist_ortho = theta[0,0] / v_ortho_length
    v_orthogonal = v_orthogonal / v_ortho_length
    v_parallel = np.array([[-v_orthogonal[1,0]],[v_orthogonal[0,0]]])
    projections = X @ v_parallel
    proj_1 = min(projections)
    proj_2 = max(projections)
    point_1 = proj_1 * v_parallel - dist_ortho * v_orthogonal
    point_2 = proj_2 * v_parallel - dist_ortho * v_orthogonal
    return point_1, point_2
fig1 = plt.figure(figsize=(5,5)) 
ax = plt.axes() 
ax.set_aspect(aspect = 'equal', adjustable = 'box')
plt.title('Logistic regression boundary')
plt.xlabel('Exam 1')
plt.ylabel('Exam 2')
plt.grid(axis='both', alpha=.25)
ax.scatter(X[:,0][idx_0], X[:,1][idx_0], s=50, c='r', marker='*', label='Not Admitted')
ax.scatter(X[:,0][idx_1], X[:,1][idx_1], s=50, c='b', marker='o', label='Admitted')
point_1, point_2 = boundary_points(X, theta)
plt.plot([point_1[0,0], point_2[0,0]],[point_1[1,0], point_2[1,0]], 'g-')
plt.legend(loc=0)
plt.show()<jupyter_output><empty_output><jupyter_text>You may have to adjust the above code to make it work with normalized data. Test set performanceNow let's apply the learned classifier to the test data we reserved in the beginning:<jupyter_code>def r_squared(y, y_pred):
    return 1 - np.square(y - y_pred).sum() / np.square(y - y.mean()).sum()
y_test_pred_soft = h(X_test, theta)
y_test_pred_hard = (y_test_pred_soft > 0.5).astype(int)

test_rsq_soft = r_squared(y_test, y_test_pred_soft)
test_rsq_hard = r_squared(y_test, y_test_pred_hard)
test_acc = (y_test_pred_hard == y_test).astype(int).sum() / y_test.shape[0]

print('Got test set soft R^2 %0.4f, hard R^2 %0.4f, accuracy %0.2f' % (test_rsq_soft, test_rsq_hard, test_acc))<jupyter_output><empty_output><jupyter_text>For classification, accuracy is probably the more useful measure of goodness of fit.  Example 2: Loan prediction datasetLet's take another example dataset and see what we can do with it.This dataset is from [Kaggle](https://www.kaggle.com/altruistdelhite04/loan-prediction-problem-dataset).The data concern loan applications. It has 12 independent variables, including 5 categorical variables. The dependent variable is the decision "Yes" or "No" for extending a loan to an individual who applied.One thing we will have to do is to clean the data, by filling in missing values and converting categorical data to reals.We will use the Python libraries pandas and sklearn to help with the data cleaning and preparation. Read the data and take a look at it<jupyter_code># Import Pandas. You may need to run "pip3 install pandas" at the console if it's not already installed

import pandas as pd

# Import the data

data_train = pd.read_csv('train_LoanPrediction.csv')
data_test = pd.read_csv('test_LoanPrediction.csv')

# Start to explore the data

print('Training data shape', data_train.shape)
print('Test data shape', data_test.shape)

print('Training data:\n', data_train)

# Check for missing values in the training and test data

print('Missing values for train data:\n------------------------\n', data_train.isnull().sum())
print('Missing values for test data \n ------------------------\n', data_test.isnull().sum())<jupyter_output><empty_output><jupyter_text>Handle missing valuesWe can see from the above table that the `Married` column has 3 missing values in the training dataset and 0 missing values in the test dataset.Let's take a look at the distribution over the datasets then fill in the missing values in approximately the same ratio.You may be interested to look at the [documentation of the Pandas `fillna()` function](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.fillna.html). It's great!<jupyter_code># Compute ratio of each category value
# Divide the missing values based on ratio
# Fillin the missing values
# Print the values before and after filling the missing values for confirmation

print(data_train['Married'].value_counts())

married = data_train['Married'].value_counts()
print('Elements in Married variable', married.shape)
print('Married ratio ', married[0]/sum(married.values))

def fill_martial_status(data, yes_num_train, no_num_train):        
    data['Married'].fillna('Yes', inplace = True, limit = yes_num_train)
    data['Married'].fillna('No', inplace = True, limit = no_num_train)  

fill_martial_status(data_train, 2, 1)
print(data_train['Married'].value_counts()) 
print('Missing values for train data:\n------------------------\n', data_train.isnull().sum())<jupyter_output><empty_output><jupyter_text>Now the number of examples missing the `Married` attribute is 0.Let's complete the data processing based on examples given and logistic regression model on training dataset. Then we'll get the model's accuracy (goodness of fit) on the test dataset.Here is another example of filling in missing values for the `Dependents` (number of children and other dependents)attribute. We see that categorical values are all numeric except one value "3+"Let's create a new category value "4" for "3+" and ensure that all the data is numeric:<jupyter_code>print(data_train['Dependents'].value_counts())
dependent = data_train['Dependents'].value_counts()

print('Dependent ratio 1 ', dependent['0'] / sum(dependent.values))
print('Dependent ratio 2 ', dependent['1'] / sum(dependent.values))
print('Dependent ratio 3 ', dependent['2'] / sum(dependent.values))
print('Dependent ratio 3+ ', dependent['3+'] / sum(dependent.values))

def fill_dependent_status(num_0_train, num_1_train, num_2_train, num_3_train, num_0_test, num_1_test, num_2_test, num_3_test):        
    data_train['Dependents'].fillna('0', inplace=True, limit = num_0_train)
    data_train['Dependents'].fillna('1', inplace=True, limit = num_1_train)
    data_train['Dependents'].fillna('2', inplace=True, limit = num_2_train)
    data_train['Dependents'].fillna('3+', inplace=True, limit = num_3_train)
    data_test['Dependents'].fillna('0', inplace=True, limit = num_0_test)
    data_test['Dependents'].fillna('1', inplace=True, limit = num_1_test)
    data_test['Dependents'].fillna('2', inplace=True, limit = num_2_test)
    data_test['Dependents'].fillna('3+', inplace=True, limit = num_3_test)

fill_dependent_status(9, 2, 2, 2, 5, 2, 2, 1)

print(data_train['Dependents'].value_counts())

# Convert category value "3+" to "4"

data_train['Dependents'].replace('3+', 4, inplace = True)
data_test['Dependents'].replace('3+', 4, inplace = True)<jupyter_output><empty_output><jupyter_text>Once missing values are filled in, you'll want to convert strings to numbers.Finally, here's an example of replacing missing values for a numeric attribute. Typically, we would use the mean of the attribute over the training set.<jupyter_code>print(data_train['LoanAmount'].value_counts())

LoanAmt = data_train['LoanAmount'].value_counts()

print('mean loan amount ', np.mean(data_train["LoanAmount"]))

loan_amount_mean = np.mean(data_train["LoanAmount"])

data_train['LoanAmount'].fillna(loan_amount_mean, inplace=True, limit = 22)
data_test['LoanAmount'].fillna(loan_amount_mean, inplace=True, limit = 5)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12964516606622908, 'fraction_numerical': 0.01697675566303114, 'mean_word_length': 5.212139791538934, 'pattern_counts': {'":': 0, '<': 82, '<?xml version=': 0, '>': 83, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '552296', 'license': 'CC0-1.0', 'n_tokens_mistral': 6752, 'n_tokens_neox': 6430, 'n_words': 2230, 'path': 'Labs/03-Logistic-Regression/03-Logistic-Regression.ipynb', 'repo_name': 'dsai-asia/ml', 'size': '20266'}