<jupyter_start><jupyter_text>4. (Chapter 4) Training ModelsIn this chapter we will begin our study of Machine Learning algorithms in earnest.  Upto this point you have gotten a taste of using some ML models as black boxes.You don't really know internally what they are doing or how they work.  But still, learningthe API of how to create data sets and set up models for training, and knowning ingeneral what supervisied learning is and how you train a model to fit a dataset issufficient to start working with ML.However, at some point you will be working with a dataset and you will get someweird or unexpected results.  Also there are myriads of parameters and hyper-parameters,not to mention issues with formatting and cleaning your data, and tweakingtraining procedures, etc. that will all always be a bit mysterious and/or you will onlylearn by wrote checklists, if you don't understand some of the fundamentals of theinternal implemenations of the various ML algorithms you use.Understanding whats under the hood will help you quickly understand what types of MLmodels will be promising or not for the data you have.  It will give you an understandingof the hyper-parameters the model uses, and what you might want to try for these parametersto get good results.  Understanding whats under the hood will also allow you to debugand spot issues you have with training and debugging your models.  It will allow youto perform error analysis, and better understand why your models may be working, and whythey may be failing, so that you can tweak things to improve their performance.In short, while you can be a push-button/check-list using data analyst, if you eventuallywant to improve your skills and become a better data scientist, you need to startlearning the internal working of ML algorithms.  Our goal in this class is to beginyour journey along this path.  We won't become experts on any of the ML methods we will look at, but we hope to introduce you to the fundamental concepts and innerworkings, so that you can beging to understand more deeply how the techniques work.In this chapter, and units of our course, we will begin by looking at **linear regression**(which we have used already), followed by **logistic regression**.  These units are veryimportant, as these models have basic concepts that are often similar across many typesof ML algorithms.  We will look at the idea of the **cost function** (or **fitness function**),and how this can be used with techniques like **gradient descent** to fit the parametersof a model in general to make predictions about a set of data.  Linear regressionis a technique for performing regression modeling (predict real valued outputs), whilelogistic regression is useful for classification (despite its name, which is a bitmisleading at best).  By comparing similarities and differences of cost functions andfitting models between regression and classification tasks, we will illustrate some generalprinciples of ML that apply to many types of supervisied learning models, whetherregression or classification tasks.  4.6 Logistic RegressionThe linear regression examples we have used up to this point have all been applied to regression problems.  Recallthat by regression we mean we are fitting a model or forming a hypothesis where the output can be a real valued numberwith in theory an infinite amount of precision.But besides regression problems, another common type of supervisied learning problem we normally need to work with areclassification problems, where the ouput is a set of discrete categories.  In the simplest case, the labeled outputsare binary, like yes/no, true/false or spam/ham classifications.We need to use a slightly modified approch to solve a classification problem.  Some regression algorithms can be used forclassification as well (and vice versa).  Logistic regression (despite being arguably misnamed) is commonlyused to estimate the probability an instance belongs to a particular class.  Estimating ProbabilitiesWe will start by considering a binary classification first, and we will use labels of 0 and 1 to represent the false/trueclasses respectively.  A logistic regression starts initially in the same way as a linear regression, by computing the weightedsum of the input features (blus a bias term).  But instead of outputting this result directly like for linear regression, it outpus the **logistic** (also called **logit**, also called the **sigmoid** function) of this result.\begin{equation}\hat{p} = h_{\theta}(x) = \sigma \big( x^T \theta \big)\end{equation}This equation has the same weighted sum calculation as for linear regression ( $x^T \theta$ ).  But thefunction $\sigma(\cdot)$ here denotes the logistic function.  The standard form of the Logistic/Sigmoid function is:\begin{equation}\sigma(z) = \frac{1}{1 + e^{-z}}\end{equation}<jupyter_code>def sigmoid(z):
    """The basic sigmoid or logistic function.  Called sigmoid because of its S-shaped
    output.
    """
    return 1.0 / (1.0 + np.exp(-z))
z = np.linspace(-10.0, 10.0, 250)
plt.plot(z, sigmoid(z), label=r'sigmoid $\sigma(z) = \frac{1}{1 + e^{-z}}$')
plt.xlabel('z')
plt.legend();<jupyter_output><empty_output><jupyter_text>Notice that the result from the weighted sum ($z = x^T \theta$), no matter if it ends up as a very large negativeresult or a very large positive result, it will always be squashed to some value between 0 and 1.  Results below-5 or above 5 get effectively logited to 0 and 1 respectively, and values between will vary somewhere in between 0 and 1.If we think of the output of the sigmoid function as a probability estimate, then we can interpret this outputas the probability $\hat{p}$ that the result is the positive (1) outcome for a binary classification.  Once the logisticregression model has estimated the probability $\hat{p} = h_{\theta}(x)$ that an instance $x$ belongs to thepositive class, we can use a threshold to make our final prediction:\begin{equation}\hat{y} =\begin{cases}0 & \text{if} \;\; \hat{p} < 0.5 \\1 & \text{if} \;\; \hat{p} \ge 0.5\end{cases}\end{equation}Notice that, using this threshold, when $z < 0$ then we will predict the negative case or 0, and when $z \ge 0$ wepredict the positive case 1.  Training and Cost FunctionGiven this modification of the weighted sum, we need to also modify the way we train our models.The objective of training for a binary classifier is to set the paramter vector $\theta$ so that themodel estimates high probabilities for positive instances ($y = 1$) and low probabilities for negative instances ($y = 0$).  This idea is captured by a cost function like the following\begin{equation}c(\theta) =\begin{cases}-\log(\hat{p}) & \text{if} \;\; y = 1 \\-\log(1 - \hat{p}) & \text{if} \;\; y = 0\end{cases}\end{equation}To understand this cost function, think of the two cases of $y = 1$ and $y = 0$ separately.  For the positive casewhen we want the model to predict $y = 1$ the cost looks like the following (remember that $\hat{p}$ isthe output from the sigmod, so it will always be a value between 0 and 1):<jupyter_code>p = np.linspace(0.001, 1, 200) # log not defined at 0, so start at very small value
plt.plot(p, -np.log(p), label=r'positive case where y = 1')
plt.xlabel(r'$\hat{p}$')
plt.ylabel(r'cost $c(\theta)$')
plt.legend();<jupyter_output><empty_output><jupyter_text>Again remember this is the cost function for all training instances where the correct answer is $y = 1$.When the prediction probability $\hat{p}$ is 1.0, the cost of our cost function approaches 0.This is what we want, the model is predicting well when it says there is a high probability itis the $y = 1$ case when that is true.  But when $y = 1$ and the model is predicting a smallprobability it is the positive case, we want the cost penality to be high for that prediction.Likewise you will get a similar though mirrored cost function for the $y = 0$ training cases:<jupyter_code>p = np.linspace(0, 0.9999, 200) # log not defined at 0, so start at very small value
plt.plot(p, -np.log(1.0 - p), label=r'negative case where y = 0')
plt.xlabel(r'$\hat{p}$')
plt.ylabel(r'cost $c(\theta)$')
plt.legend();<jupyter_output><empty_output><jupyter_text>Make sure you see that the cost functions make sense.  For the 2 cases, they give low costs when the model ismaking a good prediction, and the cost rapidly increases when the model is making a bad prediction.Finally, to use this cost function, as we did before for the MSE, we have to average this cost function over thewhole training set of training inputs.  We can use a little bit of mathematical trickery to write this as asingle expression, even though our cost function has been expressed as a piecewise function so far.\begin{equation}J(\theta) = -\frac{1}{m} \sum_{i=0}^{m-1} \Big[ y^{(i)} \log \big( \hat{p}^{(i)} \big) + \big( 1 - y^{(i)} \big) \log \big( 1 - \hat{p}^{(i)} \big)   \Big]\end{equation} As before we have changed to 0 based indexing here, which is different from the text, to reflect our Python implementation.Notice that inside of the summation, when $y = 1$ the left hand term will remain, but the right hand term will get zeroed out.  And likewise when $y = 0$ the left part gets zeroed out. There is no known closed form solution to compute the value of $\theta$ that minimizes this cost function(no equivalent of the normal equation).  So we cannot "solve" this analytically like for linear regression.But this cost function is convex, like the linear regression MSE cost function.  So gradient descent (or anyoptimization method) is guaranteed to find the global minimum (if the learning rate is not too large and youwait long enough). The partial derivatives of the cost function with regards to the $j^{th}$ model parameter $\theta_j$ is given by:\begin{equation}\frac{\partial}{\partial \theta_j}  J(\theta) = \frac{1}{m} \sum_{i=0}^{m-1} \Big( \sigma \big( \theta^T x^{(i)} \big) - y^{(i)}  \Big) x_j^{(i)}\end{equation} It is beyond the scope of this course to understand how the partial derivative of thecost function is derived here.  But you should go back and compare this equationto the derivative of the MSE cost function that was given for linear regression.  They arevery similar (on purpose, this is one of the reasons for picking the logithere).  In fact the only difference is that we have to transform theweighted sum term here using the sigmoid function, which we do not do forthe MSE derivative.  This means that we can use basically the same implementationof gradient descent, with only a very small modification, for bothlinear and logistic regression.Once we can calculate the gradients, we can use gradient descent in exactly thesame ways as discussed before for linear regression.  We can use batch gradient descent, stochastic ormini-batches.  Decision BoundariesThe best way to visualize a logistic regression model is to determine andplot the decision boundary that the model produces.  We will illustratethe concept of a decision boundary for classification problems using the well known iris dataset.  This dataset has only 150 instances, with 4separate features (sepal length, sepal width, petal length, petal width).  Thisis a classification problem, there are 3 classes of irises in thedataset, Setosa, Versicolor and Virginica.  We will later talk about how togeneralize logistic regression to problems with more than 2 classes.  But here, we will collapse the Setosa and Versicolor cases into one, and simplyclassify either Virginica or not Virginica.Start by getting the iris dataset from `sklearn`.<jupyter_code>from sklearn import datasets
iris = datasets.load_iris()

# this is a dictionary, here are the items of the dataset we can use
list(iris.keys())
# extract the input features, we are only using 1 feture to begin with, petal width
X = iris["data"][:, 3:]

# make the labels, 1 if virginica and 0 if not-virginica
y = (iris["target"] == 2).astype(np.int)<jupyter_output><empty_output><jupyter_text>Use `sklearn` to train a `LogisticRegression` model.<jupyter_code># create and train a linear model using Logistic Regression
from sklearn.linear_model import LogisticRegression
log_reg = LogisticRegression()
log_reg.fit(X, y)<jupyter_output><empty_output><jupyter_text>Here since we only use a single feature, the decision boundary will basically occurat some specific petal width.  This will be the point where theestimated probability of the learned hypothesis is at 0.5 (our thresholdvalue).Lets visualize this boundary, and plot some of the instances in ourtraining set as well.<jupyter_code>X_new = np.linspace(0, 3, 1000).reshape(-1, 1)
y_proba = log_reg.predict_proba(X_new)
plt.plot(X_new, y_proba[:,1], 'g-', label="Virginica");

# rest of figure shown in textbook
plt.plot(X_new, 1.0 - y_proba[:,1], 'b--', label='not Virginica')

# determine the decision boundary value and visualize it
decision_index = min(np.where(y_proba[:,1] > 0.5)[0])
decision_value = X_new[decision_index]
print("Decision boundary value", decision_value)
plt.plot([decision_value, decision_value], [0, 1], 'k:')

# plot Virginica and not Virginica points to visualize how well the decision boundary works
not_virginica_width = X[y == 0]
not_virginica_y = np.zeros(not_virginica_width.shape)
plt.plot(not_virginica_width, not_virginica_y, 'bs')

virginica_width = X[y == 1]
virginica_y = np.ones(virginica_width.shape)
plt.plot(virginica_width, virginica_y, 'g^')

# add labels and such
plt.annotate(r'Decision Boundary',
             xy=(decision_value, 0.1), xycoords='data',
             xytext=(-70, 0), textcoords='offset points', fontsize=18)
plt.arrow(decision_value, 0.07, -0.33, 0, color='blue', linewidth=1.0, width=0.01)
plt.arrow(decision_value, 0.93, 0.33, 0, color='green', linewidth=1.0, width=0.01)

plt.xlabel('Petal Width (cm)')
plt.ylabel('Probability')
plt.legend();<jupyter_output>Decision boundary value [1.66066066]<jupyter_text>The petal widths of Virginica range from 1.4 cm to 2.5 cm, while the other iris flowers havesmaller petal lengths, ranging from 0.1 cm to 1.8 cm.  Notice that there is a bit of overlap,we cant get a perfect prediction with this attribute alone (petal width) of the irisinstances.  There is a decision boundary at around 1.66 cm that the classifier learns here.  This will be the location of the learned parameter models where the probability for an input of thepetal width features results in a probability value of 0.5.  The black dotted line representsthe decision boundary, the location where the model predicted probability is 0.5.  A little smallerthan this decision boundary, and the model will predict the iris is not a Virginica.  A little bigger and it predicts it is Virginica.The 0.5 threshold is used as the decision boundary.  Thus if we ask the logistic regressor to give us its best prediction, we will get a definitivie answer. For example:<jupyter_code>log_reg.predict([[1.7], [1.5]])<jupyter_output><empty_output><jupyter_text>The logistic regressor we used previously with a single feature can be used with 2 or manymore features.  But the logistic regressor we are currently using creates a linear modelfor the data.  This means that the decision boundary will always be a straight line (or hyperplanein 3 or more dimensions).  For example, if we train the Virginica/Not Virginica classifier on two features, petal width and petal length, we can see that the decision boundary that resultsis a line in the 2 dimensions.<jupyter_code># extract the input features, use 2 features, petal length and petal width this time
X = iris["data"][:, 2:]

# make the labels, 1 if virginica and 0 if not-virginica
y = (iris["target"] == 2).astype(np.int)
log_reg = LogisticRegression()
log_reg.fit(X, y)
# create a mesh so we can use a contour plot to visualize the decision boundary
x0, x1 = np.meshgrid(
        np.linspace(2.9, 7, 500).reshape(-1, 1),
        np.linspace(0.8, 2.7, 200).reshape(-1, 1),
    )
X_new = np.c_[x0.ravel(), x1.ravel()]

# calculate probability on the mesh for the contour plot
y_proba = log_reg.predict_proba(X_new)

# plot the actual points, green triangles for Virginica, blue squares are not Virginica
plt.figure(figsize=(12, 12))
plt.plot(X[y==0, 0], X[y==0, 1], "bs")
plt.plot(X[y==1, 0], X[y==1, 1], "g^")

# create contour plot
zz = y_proba[:, 1].reshape(x0.shape)
contour = plt.contour(x0, x1, zz, cmap=plt.cm.brg)
plt.clabel(contour, inline=1, fontsize=12)



# plot the decision boundary as a dashed thick black line
left_right = np.array([2.9, 7])
boundary = -(log_reg.coef_[0][0] * left_right + log_reg.intercept_[0]) / log_reg.coef_[0][1]
plt.plot(left_right, boundary, "k--", linewidth=3)

# add in some labels and what not
plt.text(3.5, 1.5, "Not Virginica", fontsize=14, color="b", ha="center")
plt.text(6.5, 2.3, "Virginica", fontsize=14, color="g", ha="center")
plt.xlabel("Petal length", fontsize=14)
plt.ylabel("Petal width", fontsize=14)
plt.axis([2.9, 7, 0.8, 2.7]);<jupyter_output><empty_output><jupyter_text>Just like the linear regression models we discussed above, Logistic Regression models can be regularized using $\ell_1$ or $\ell_2$ penalities.  Scikit-learn adds an $\ell_2$ penality by default.  Softmax RegressionThe logistic regression model can be generalized to model multiple classes directly.  One approach is to train multiplelogistic classifiers on all of the combinations of class/not class and choose the one that gives the highest probability when we want to make predictions (discussed previously).Another approach supports multi-classes directly, called the **softmax regression** or **multinomial logistic regression**.  The ideais similar.  The softmax regression first computes a score $s_k(x)$ for each class $k$, then it estimates the probability of eachclass by applying the **softmax function** (also called the **normalized exponential**).  The equation to compute the initial score foreach class should look familiar, as it is just like the equation for linear regression prediciton.\begin{equation}s_k(x) = x^T \theta^{(k)}\end{equation}This means that each class $k$ has its own dedicated parameter vector $\theta^{(k)}$.  So we are essentially again proposing separate classifiersfor each class.  All these vectors are typically stored as rows in a **parameter matrix** $\Theta$.Once you have computed the score of every class for the instance $x$, you can estimate the probability $\hat{p}_k$ that the instance belongs to class$k$ by running the scores through the softmax function.  It computes the exponential of every score, then normalizes them (dividing by the sum of all the exponentials).  The scores are generally called logits or log-odds.\begin{equation}\hat{p}_k = \sigma(s(x))_k = \frac{\exp \big(s_k(x) \big)}{\sum_{j=1}^{K} \exp \big( s_j(x) \big)}\end{equation}Where - $K$ is the number of classes, if we have to we can renumber using 0 based indexing for implementation in Python- $s(x)$ is a vector containing the scores of each class for the instance $x$.- $\sigma \big( s(x) \big)_k $ is the estimated probability that the instance $x$ belongs to class $k$ given the scores of each class for that instance.Just like the logistic regression classifier, the softmax regression classifier predicts the class with the highest estimated probability:\begin{equation}\hat{y} = \underset{k}{\text{argmax}} \; \sigma(s(x))_k = \underset{k}{\text{argmax}} \; s_k(x) = \underset{k}{\text{argmax}} \; \Big( \big(\theta^{(k)} \big)^T x \Big) \end{equation} An example of training using the softmax regression classifier should make this more clear.  For training, we will minimize thedefined probability function $\hat{p}_k$ as our objective function.  The following equation is known as the **cross entropy**cost function.\begin{equation}J(\Theta) = -\frac{1}{m} \sum_{i=0}^{m-1} \sum_{k=0}^{K-1} y_k^{(i)} \log \big( \hat{p}_k^{(i)} \big)\end{equation}Here $y_k^{(i)}$ is the target probability that the $i^{th}$ instance belongs to class $k$.  It will be either equal to 1 or 0 dependingon whether the instance $i$ belongs to the class or not.Notice that when there are just two classes ($K=2$) this cost function is equivalent to the logistic regression's cost function.Given this cost function, we need the gradient vector of the cost function with regards to $\theta^{(k)}$.  The derivation is the sameas for the logistic regression.  We end up with the following formula for calculating the gradients\begin{equation}\bigtriangledown_{\theta^{(k)}} J(\Theta) = \frac{1}{m} \sum_{i=0}^{m-1} \big( \hat{p}_k^{(i)} - y_k^{(i)}  \big) x^{(i)}\end{equation}Given this expression for the gradient, we can use gradient descent to perform the softmax regression.  Of course, this will all bedone for us using scikit-learn.  Lets use scikit-learn's `LogisticRegression` class to classify all 3 iris flower classes.Scikit-learn `LogisticRegression` uses one-versus-all by default when you train it on more than two classes, but you can set the`multi_class` hyperparameter to `"multinomial"` to switch it to softmax regression instead.  You must also specify a solver(gradient descent algorith) that supports softmax regression, such as the `"lbgs"` solver.  It also applies $\ell_2$ regularizationby default, which you can control using the hyperparameter $C$.<jupyter_code>X = iris["data"][:, (2,3)] # still only using petal length and petal width
y = iris["target"] # but use all of the targets
softmax_reg = LogisticRegression(multi_class="multinomial", solver="lbfgs", C=10)
softmax_reg.fit(X, y)<jupyter_output><empty_output><jupyter_text>So once the mode is trained, we can ask it to predict.  Like with one-versus-all, asking to `predict` will give the predicted class label.And we can use the `predict_proba` function to see the individual softmax probability estimates for each class the model is trained with.<jupyter_code># prediction for an iris with 5cm petal length and 2cm petal width
softmax_reg.predict([[5,2]])
# probabilities, these are in the order defined by y=0,1,2
# which is Setosa, Versicolor and Virginica respecitvely
softmax_reg.predict_proba([[5,2]])<jupyter_output><empty_output><jupyter_text>Softmax gives an almost 94.3% probability prediction that this iris should be in the Virginica class, with a 5.7% probability that it is the versicolor.We can also plot the resulting decision boundary reached by this softmax regressor.  Notice that the decision boundaries between any 2classes are still linear.<jupyter_code># create a meshgrid so we can make a contour plot of the decision boundaries
x0, x1 = np.meshgrid(
        np.linspace(0, 8, 500).reshape(-1, 1),
        np.linspace(0, 3.5, 200).reshape(-1, 1),
    )
X_new = np.c_[x0.ravel(), x1.ravel()]


# make predictions over the grid for the contour plot
y_proba = softmax_reg.predict_proba(X_new)
y_predict = softmax_reg.predict(X_new)

# plot the contour map
plt.figure(figsize=(12, 12))
zz1 = y_proba[:, 1].reshape(x0.shape)
zz = y_predict.reshape(x0.shape)
from matplotlib.colors import ListedColormap
custom_cmap = ListedColormap(['#fafab0','#9898ff','#a0faa0'])

plt.contourf(x0, x1, zz, cmap=custom_cmap)
contour = plt.contour(x0, x1, zz1, cmap=plt.cm.brg)

# Plot the actual data points, color and shape coded
plt.plot(X[y==2, 0], X[y==2, 1], "g^", label="Iris-Virginica")
plt.plot(X[y==1, 0], X[y==1, 1], "bs", label="Iris-Versicolor")
plt.plot(X[y==0, 0], X[y==0, 1], "yo", label="Iris-Setosa")

# add labels and plot information
plt.clabel(contour, inline=1, fontsize=12)
plt.xlabel("Petal length", fontsize=14)
plt.ylabel("Petal width", fontsize=14)
plt.legend(loc="center left", fontsize=14)
plt.axis([0, 7, 0, 3.5]);<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09103594080338266, 'fraction_numerical': 0.015644820295983086, 'mean_word_length': 4.909795102448776, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 45, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '239389', 'license': 'CC-BY-3.0', 'n_tokens_mistral': 7073, 'n_tokens_neox': 6629, 'n_words': 3177, 'path': 'lectures/homl/homl-04-5-logistic-regression.ipynb', 'repo_name': 'tgrasty/CSCI574-Machine-Learning', 'size': '23650'}