This is a notebook explaining the mathematics of Linear and Logistic Regression. It is adapted from by [blog](http://ataspinar.com/2016/03/28/regression-logistic-regression-and-maximum-entropy). All of the code is available on [GitHub](https://github.com/taspinar/siml).


<h1 style="text-align: center;"> 1. Regression Analysis</h1>
Regression Analysis is the field of mathematics where the goal is to find a function which best [correlates](http://guessthecorrelation.com/) with a dataset. Let's say we have a dataset containing $ n $ datapoints; 

<p style="text-align: center;">$ X = ( x^{(1)}, x^{(2)}, .., x^{(n)} ) $. </p>

For each of these (input) datapoints there is a corresponding (output) $ y^{(i)}$-value. 

Here, the $ x$-datapoints are called the [independent variables](https://en.wikipedia.org/wiki/Dependent_and_independent_variables) and $ y$ the dependent variable; 
the value of $ y^{(i)} $ depends on the value of $ x^{(i)} $, while the value of $ x^{(i)}$ may be freely chosen without any restriction imposed on it by any other variable.

The goal of Regression analysis is to find a function $ f(X) $ which can best describe the correlation between $ X $ and $ Y $. In the field of Machine Learning, this function is called the hypothesis function and is denoted as $ h_{\theta}(x) $. If we can find such a function, we can say we have successfully built a Regression model. 

<br>



If the input-data lives in a 2D-space, this boils down to finding a curve which fits through the data points. In the 3D case we have to find a plane and in higher dimensions a hyperplane.

To give an example, let's say that we are trying to find a predictive model for the success of students in a course called Machine Learning. We have a dataset $ Y $ which contains the final grade of $ n $ students. Dataset $ X $ contains the values of the independent variables. Our initial assumption is that the final grade only depends on the studying time. The variable $ x^{(i)} $ therefore indicates how many hours student $ i $ has studied. The first thing we would do is visualize this data:



If the results looks like the figure on the left, then we are out of luck. It looks like the points are distributed randomly and there is no correlation between $ Y$ and $ X$ at all. However, if it looks like the figure on the right, there is probably a strong correlation and we can start looking for the function which describes this correlation.

&nbsp;

This function could for example be:
<p style="text-align: center;">$ h_{\theta}(X) =  \theta_0+ \theta_1 \cdot x $</p>
<p style="text-align: left;">or</p>
<p style="text-align: center;">$ h_{\theta}(X) = \theta_0 + \theta_1 \cdot x^2 $</p>
where $ \theta $ are the dependent parameters of our model.

<br><br><br>

<h2> 1.2 Multivariate Regression</h2>
<br>
In evaluating the results from the previous section, we may find the results unsatisfying; the function does not correlate with the datapoints strongly enough. Our initial assumption is probably not complete. Taking only the studying time into account is not enough. 
The final grade does not only depend on the studying time, but also on how much the students have slept the night before the exam. Now the dataset contains an additional variable which represents the sleeping time. 

Our dataset is then given by $ X = ( (x_1^{(1)}, x_2^{(1)}), (x_1^{(2)}, x_2^{(2)}), .., (x_1^{(n)}, x_2^{(n)}) ) $. In this dataset $ x_1^{(i)} $ indicates how many hours student $ i $ has studied and $ x_2^{(i)} $ indicates how many hours he has slept.



This is an example of [Multivariate Regression](https://en.wikipedia.org/wiki/Linear_regression#Simple_and_multiple_regression). The function has to include both variables. For example:
<p style="text-align: center;">$ h_{\theta}(x) = \theta_0 + \theta_1 \cdot x_1 + \theta_2 \cdot x_2$</p>
or
<p style="text-align: center;">$ h_{\theta}(x) = \theta_0 + \theta_1 \cdot x_1 + \theta_2 \cdot x_2^3 $.</p>


<h2> 1.3 Linear vs Non-Linear Regression</h2>
<br>
All of the above examples are examples of linear regression. We have seen that in some cases $ y^{(i)} $ depends on a linear form of $ x^{(i)} $, but it can also depend on some power of $ x^{(i)}$, or on the log or any other form of $ x^{(i)}$. However, in all cases its dependence on the parameters $\theta$ was linear.

So, what makes linear regression linear is not that $ Y$ depends in a linear way on $ X $, but that it depends in a linear way on $ \theta$. 
$ Y $ needs to be [linear](http://www.investopedia.com/terms/l/linearrelationship.asp) with respect to the model-parameters $ \theta $. Mathematically speaking it needs to satisfy the [superposition principle](http://www.cut-the-knot.org/do_you_know/superposition.shtml).  

Examples of nonlinear regression would be:
<p style="text-align: center;">$ h_{\theta}(x) = \theta_0 + x_1^{\theta_1} $</p>
or
<p style="text-align: center;">$ h_{\theta}(x) = \theta_0 + \theta_1 / x_1 $</p>
&nbsp;

The reason why the distinction is made between linear and nonlinear regression is that nonlinear regression problems are more difficult to solve and therefore more computational intensive algorithms are needed.

Linear regression models can be written as a linear system of equations, which can be solved by finding the closed-form solution $ \theta = ( X^TX )^{-1}X^TY $ with Linear Algebra. See <a href="http://www.stat.purdue.edu/~jennings/stat514/stat512notes/topic3.pdf" target="_blank">these</a> statistics notes for more on solving linear models with linear algebra.

As discussed before, such a closed-form solution can only be found for linear regression problems. However, even when the problem is linear in nature, we need to take into account that calculating the inverse of a $ n $ by $ n $ matrix has a time-complexity of $ O(n^3) $. This means that for large datasets ($ n \gt 10.000 $ ) finding the closed-form solution will take more time than solving it iteratively (gradient descent method) as is done for nonlinear problems. So solving it iteratively is usually preferred for larger datasets, even if it is a linear problem.



<h2> 1.4 Gradient Descent</h2>
<br>
The Gradient Descent method is a general optimization technique in which we try to find the value of the parameters $ \theta $ with an iterative approach.

First, we construct a <a href="https://en.wikipedia.org/wiki/Loss_function" target="_blank">cost function</a> (also known as loss function or error function) which gives the difference between the values of the hypothesis function $h_{\theta}(x)$ (the values you expect $ Y$ to have with the current values of $ \theta$ ) and the actual values of $ Y $. The better your estimation of $ \theta $ is, the better the values of $ h_{\theta}(x) $ will approach the values of $ Y$ and the smaller the cost function will be.

Usually, the cost function is expressed as the squared error of the difference between these functions:
<p style="text-align: center;">$ J(x) = \frac{1}{2n} \sum_i^n ( h_{\theta}(x^{(i)}) - y^{(i)} )^2 $</p>
&nbsp;

At each iteration we choose new values for the parameters $ \theta$, and move towards the 'true' values of these parameters, i.e. the values which make this cost function as small as possible. The direction in which we have to move is the negative gradient direction;

<p style="text-align: center;"> $ \Delta\theta = - \alpha \frac{d}{d\theta} J(x) $.</p>

The reason for this is that  a function's value decreases fastest if we move towards the direction of the negative gradient (the <a href="https://en.wikipedia.org/wiki/Directional_derivative" target="_blank">directional derivative</a> is maximal in the direction of the gradient).

Taking all this into account, this is how gradient descent works:
<ul>
	<li>Make an initial but intelligent guess for the values of the parameters $ \theta $.</li>
	<li>Keep iterating while the value of the cost function has not met your criteria:
<ul>
	<li>With the current values of $ \theta $, calculate the gradient of the cost function J  ( $ \Delta \theta = - \alpha \frac{d}{d\theta} J(x)$ ).</li>
	<li>Update the values for the parameters $ \theta := \theta + \alpha \Delta \theta $</li>
	<li>Fill in these new values in the hypothesis function and calculate again the value of the cost function;</li>
</ul>
</li>
</ul>

Just as important as the initial guess of the parameters is the value you choose for the learning rate $ \alpha $. This learning rate determines how fast you move along the slope of the gradient. If the selected value of this learning rate is too small, it will take too many iterations before you reach your convergence criteria. If this value is too large, you might overshoot and not converge.

<h1 style="text-align: center;"> 2. Logistic Regression </h1>
Logistic Regression is similar to (linear) regression, but adapted for the purpose of classification. The difference is small; for Logistic Regression we also have to iteratively apply the gradient descent method to estimate the values of the parameter $ \theta$. And again, during the iteration, the values are estimated by taking the gradient of the cost function. And also, the cost function is given by the squared error of the difference between the hypothesis function $ h_{\theta}(x)$ and $ Y$. The major difference however, is the hypothesis function itself.

In order to understand the hypothesis function of Logistic Regression, we must first understand the idea behind classification. 

When you want to classify something, there are a limited number of classes it can belong to. And for each of these possible classes there can only be two states for $ y^{(i)} $;
either $ y^{(i)} $ belongs to the specified class and $ y = 1 $, or it does not belong to the class and $ y = 0 $. Even though the output values $ Y $ are binary, the independent variables $ X $ are still continuous. So, we need a function which has as input a large set of continuous variables $ X $ and for each of these variables produces a binary output. This function, the hypothesis function, has the following form:

$ h_{\theta} = \frac{1}{1 + \exp(-z)} = \frac{1}{1 + \exp(-\theta x)} $.

This function is also known as the <a href="https://en.wikipedia.org/wiki/Logistic_function" target="_blank">logistic function</a>, which is a part of the <a href="https://en.wikipedia.org/wiki/Sigmoid_function" target="_blank">sigmoid function</a> family. These functions are widely used in the natural sciences because they provide the simplest model for population growth. However, the reason why the logistic function is used for classification in Machine Learning is its 'S-shape'.



As you can see this function is bounded in the y-direction by 0 and 1. If the variable $ z$ is very negative, the output function will go to zero (it does not belong to the class). If the variable $ z$ is very positive, the output will be one and it does belong to the class. 
(Such a function is called an <a href="https://en.wikipedia.org/wiki/Indicator_function" target="_blank">indicator function</a>.)

The question then is, what will happen to input values which are neither very positive nor very negative, but somewhere 'in the middle'. We have to define a decision boundary, which separates the positive from the negative class. Usually this decision boundary is chosen at the middle of the logistic function, namely at $ z = 0 $ where the output value $ y$ is $ 0.5$.

\begin{equation}
y =\begin{cases}
1, \text{if $z \gt  0 $}.\\
0, \text{if $z \lt 0 $}.
\end{cases}
\end{equation}

As we can see in the formula of the logistic function, $ z = \theta \cdot x $. Meaning, the dependent parameter $ \theta$ (also known as the feature), maps the input variable $ x$  to a position on the $ z$-axis. With its $ z$-value,  we can use the logistic function to calculate the $ y$ -value. If this $ y$-value $ \gt 0.5 $ we assume it does belong in this class and vice versa.

So the feature $ \theta $ should be chosen such that it predicts the class membership correctly. It is therefore essential to know which features are useful for the classification task. Once the appropriate features are <a href="https://en.wikipedia.org/wiki/Feature_selection" target="_blank">selected</a> , gradient descent can be used to find the optimal value of these features.

How can we do gradient descent with this logistic function? Except for the hypothesis function having a different form, the gradient descent method is exactly the same. We again have a cost function, of which we have to iteratively take the gradient w.r.t. the feature $ \theta $ and update the feature value at each iteration.

This cost function is given by

&nbsp;

\begin{equation}
\begin{split}
J(x) = -\frac{1}{2n} \sum_i^n \left(  y^{(i)} log( h_{\theta}(x^{(i)})) + (1-y^{(i)})log(1-h_{\theta}(x^{(i)})) \right) \\
 = -\frac{1}{2n} \sum_i^n \left( y^{(i)} log(\frac{1}{1+exp(-\theta x)}) + (1-y^{(i)})log(1-\frac{1}{1+exp(-\theta x)}) \right)
\end{split}
\end{equation}

&nbsp;

We know that:
$ log(\frac{1}{1+exp(-\theta x)}) = log(1) - log(1+exp(-\theta x)) = - log(1+exp(-\theta x))$

and

$ log(1-\frac{1}{1+exp(-\theta x)}) = log( \frac{exp(-\theta x)}{1+exp(-\theta x)}) $

$ = log(exp(-\theta x)) - log(1+exp(-\theta x)) $

$ = -\theta x^{(i)} -  log(1+exp(-\theta x)) $

Plugging these two equations back into the cost function gives us:

\begin{equation}
\begin{split}
J(x)  = - \frac{1}{2n} \sum_i^n \left( - y^{(i)} log(1+exp(-\theta x)) - (1-y^{(i)})(\theta x^{(i)} +  log(1+exp(-\theta x))) \right) \\
 = - \frac{1}{2n} \sum_i^n \left(  y^{(i)} \theta x^{(i)} -\theta x^{(i)} -log(1+exp(-\theta x)) \right)
\end{split}
\end{equation}
&nbsp;

The gradient of the cost function with respect to $ \theta $ is given by

\begin{align}
\frac{d}{d\theta} J(x) = - \frac{1}{2n} \sum_i^n \left(  y^{(i)} x^{(i)} - x^{(i)} + x^{(i)} \frac{ exp(-\theta x)}{1+exp(-\theta x)} \right) \\
 = - \frac{1}{2n} \sum_i^n \left( x^{(i)} ( y^{(i)} - 1 +\frac{exp(-\theta x)}{1+exp(-\theta x)} ) \right) \\
 = - \frac{1}{2n} \sum_i^n  \left( x^{(i)} ( y^{(i)} - \frac{1}{1+exp(-\theta x)} ) \right) \\
 = - \frac{1}{2n} \sum_i^n  \left( x^{(i)} ( y^{(i)} - h_{\theta}(x^{(i)}) )\right)
\end{align}

&nbsp;

So the gradient of the seemingly difficult cost function, turns out to be a much simpler equation. And with this simple equation, gradient descent for Logistic Regression is again performed in the same way:

<ul>
	<li>Make an initial but intelligent guess for the values of the parameters $ \theta $.</li>
	<li>Keep iterating while the value of the cost function has not met your criteria:
<ul>
	<li>With the current values of $ \theta $, calculate the gradient of the cost function J  ( $ \Delta \theta = - \alpha \frac{d}{d\theta} J(x)$ ).</li>
	<li>Update the values for the parameters $ \theta := \theta + \alpha \Delta \theta $</li>
	<li>Fill in these new values in the hypothesis function and calculate again the value of the cost function;</li>
</ul>
</li>
</ul>



<h2> 2.2 Logistic Regression - Example 1</h2>
<br>
Now that we have looked at the theory, let's lets have a look at an example of Logistic Regression. We will start out with the example of students passing a course or not.

Let's generate some data points. There are $n = 300$ students participating in the course Machine Learning and whether a student $ i $ passes ( $ y_i = 1 $) or not ( $ y_i = 0$ ) depends on two variables;
<ul>
 	<li>$x_i^{(1)}$ : how many hours student $ i $ has studied for the exam.</li>
 	<li>$x_i^{(2)}$ : how many hours student $ i $ has slept the day before the exam.</li>
</ul>
&nbsp;

In our example, the results are pretty binary; everyone who has studied less than 4 hours fails the course, as well as everyone whose studying time + sleeping time is less than or equal to 13 hours ($x_i^{(1)} + x_i^{(2)} \leq 13$). 




```python
import random
import numpy as np

def func2(x_i):
    if x_i[1]+x_i[2] >= 13: 
        return 0
    else:
        return 1

def generate_data2(no_points):
    X = np.zeros(shape=(no_points, 3))
    Y = np.zeros(shape=no_points)
    for ii in range(no_points):
        X[ii][0] = 1
        X[ii][1] = random.random()*9+0.5
        X[ii][2] = random.random()*9+0.5
        Y[ii] = func2(X[ii])
    return X, Y

X, Y = generate_data2(300)

```

The results looks like this (the green dots indicate a pass and the red dots a fail)





```python
import math
import numpy as np
 
def z_to_y(z_i):
    return 1 if z_i > 0.5 else 0
 
def determine_correct_guesses(X, Y, theta, m):
    determined_z = [np.dot(theta, X[ii]) for ii in range(m)]
    determined_Y = [z_to_y(elem) for elem in determined_z]
    correct_guesses = [y_i - det_y_i for y_i, det_y_i in zip(Y, determined_Y)]
    return correct_guesses
 
def hypothesis(theta, x_i):
    z = np.dot(theta, x_i)
    sigmoid = 1.0 / (1.0 + math.exp(-1.0*z))
    return sigmoid
    
def gradient_descent(X, Y, theta, alpha, m, number_of_iterations):
    for iter in range(0,number_of_iterations):
        cost = (-1.0/m)*sum([Y[ii]*math.log(hypothesis(theta, X[ii]))+(1-Y[ii])*math.log(1-hypothesis(theta, X[ii])) for ii in range(m)])
        grad = (-1.0/m)*sum([X[ii]*(Y[ii]-hypothesis(theta, X[ii])) for ii in range(m)])
        theta = theta - alpha * grad
        correct_guesses = determine_correct_guesses(X, Y, theta, m).count(0)
        if iter % 100 == 0: print("iteration %s : cost %s : correct_guesses %s" % (iter, cost, correct_guesses))
    return theta
 
numIterations = 1500
alpha = 0.6
m,n = np.shape(X)
theta = np.ones(n)
theta = gradient_descent(X, Y, theta, alpha, m, numIterations)

```

    iteration 0 : cost 3.12146529928 : correct_guesses 241
    iteration 100 : cost 0.852070281571 : correct_guesses 107
    iteration 200 : cost 0.852743840457 : correct_guesses 241
    iteration 300 : cost 0.168394744623 : correct_guesses 225
    iteration 400 : cost 1.84940208932 : correct_guesses 257
    iteration 500 : cost 1.26800878408 : correct_guesses 222
    iteration 600 : cost 0.97047975629 : correct_guesses 214
    iteration 700 : cost 0.780941221589 : correct_guesses 217
    iteration 800 : cost 0.571828461512 : correct_guesses 227
    iteration 900 : cost 0.353647771742 : correct_guesses 247
    iteration 1000 : cost 0.203740088896 : correct_guesses 261
    iteration 1100 : cost 0.140443855068 : correct_guesses 268
    iteration 1200 : cost 0.109774575449 : correct_guesses 275
    iteration 1300 : cost 0.0910747686516 : correct_guesses 280
    iteration 1400 : cost 0.078624219286 : correct_guesses 285


<h2> 2.3 Logistic Regression - Example 2: Medical Data</h2>
<br>
Now  that the concept of Logistic Regression is a bit more clear, let's classify real-world data!
The <a href="https://www.umass.edu/statdata/statdata/stat-logistic.html" target="_blank">University of Massachusetts</a> provides some datasets which are ideal to perform Logistic Regression on. They are small (so my small laptop can also perform it in a reasonable amount of time) and there are various datasets with different (amount of) features.

The dataset "<a href="https://www.umass.edu/statdata/statdata/data/myopia.dat" target="_blank">myopia.dat</a>" contains the medical data of 618 subjects, and has 15 features describing the characteristics of each subject. We can read this data in as follows:



```python
datafile = '../datasets/myopia.dat'
file = open(datafile, 'r')
no_points = 618
X = np.zeros(shape=(no_points, 16))
Y = np.zeros(shape=no_points)
rownum = 0


for line in file:
    line = line.split()
    Y[rownum] = int(line[2])
    X[rownum][0] = 1
    X[rownum][1] = int(line[3])/10.0
    X[rownum][2] = int(line[4])
    X[rownum][3] = float(line[5])
    X[rownum][4] = float(line[6])/10.0
    X[rownum][5] = float(line[7])/10.0
    X[rownum][6] = float(line[8])/10.0
    X[rownum][7] = float(line[9])/10.0
    X[rownum][8] = int(line[10])/10.0
    X[rownum][9] = int(line[11])/10.0
    X[rownum][10] = int(line[12])/10.0
    X[rownum][11] = int(line[13])/10.0
    X[rownum][12] = int(line[14])/10.0
    X[rownum][13] = int(line[15])/10.0
    X[rownum][14] = int(line[16])
    X[rownum][15] = int(line[17])
    rownum+=1

numIterations = 1500
alpha = 0.5
m,n = np.shape(X)
theta = np.ones(n)
theta = gradient_descent(X, Y, theta, alpha, m, numIterations)

```

    iteration 0 : cost 11.9698579441 : correct_guesses 81
    iteration 100 : cost 0.271623687304 : correct_guesses 543
    iteration 200 : cost 0.255773830584 : correct_guesses 550
    iteration 300 : cost 0.251054357413 : correct_guesses 551
    iteration 400 : cost 0.249151329807 : correct_guesses 552
    iteration 500 : cost 0.2482092657 : correct_guesses 553
    iteration 600 : cost 0.247665724849 : correct_guesses 553
    iteration 700 : cost 0.247315253506 : correct_guesses 553
    iteration 800 : cost 0.247071379562 : correct_guesses 553
    iteration 900 : cost 0.24689286954 : correct_guesses 553
    iteration 1000 : cost 0.246757673012 : correct_guesses 553
    iteration 1100 : cost 0.246652753284 : correct_guesses 553
    iteration 1200 : cost 0.246569757092 : correct_guesses 553
    iteration 1300 : cost 0.246503005543 : correct_guesses 553
    iteration 1400 : cost 0.24644847617 : correct_guesses 553


While reading in, all values are normalized to a value around 1. This is done in order to speed up the calculations and to ensure that you never take the logarithm of a zero value (log functions don't really like that).

Once the data has been read into the $X$ and $Y$ matrices, logistic regression can be applied. This simple algorithm for logistic regression correctly classifies ~550 of the 618 subjects, giving it an accuracy of ~90%.



```python

```

4c42a8bb3363a356031a22f2334367e757639eb9
ipynb
Jupyter Notebook
notebooks/Linear Regression, Logistic Regression.ipynb
elenita1221/siml
3e75f7f09b7c72a55d3475f7b373dc23cf9ec38d
MIT
notebooks/Linear Regression, Logistic Regression.ipynb
elenita1221/siml
3e75f7f09b7c72a55d3475f7b373dc23cf9ec38d
MIT
notebooks/Linear Regression, Logistic Regression.ipynb
elenita1221/siml
3e75f7f09b7c72a55d3475f7b373dc23cf9ec38d
MIT
Qwen/Qwen-72B
1. YES
2. YES