```python
import numpy as np
import numpy.linalg as LA
from sklearn import datasets
from sklearn.linear_model import LogisticRegression
import matplotlib.pyplot as plt
```

# Logistic Regression

## Maximim Log-likelihood

Here we will use logistic regression to conduct a binary classification.  The logistic regression process will be formulated using Maximum Likelihood estimation.  To begin, consider a random variable $y \in\{0,1\}$.  Let the $p$ be the probability that $y=1$. Given a set of $N$ trials, the liklihood of the sequence $y_{1}, y_{2}, \dots, y_{N}$ is given by:

$\mathcal{L} = \Pi_{i}^{N}p^{y_{i}}(1-p)^{1 - y_{i}}$

Given a set of labeled training data, the goal of maximum liklihood estimation is to determine a probability distribution that best recreates the empirical distribution of the training set. 

The log-likelihood is the logarithmic transformation of the likelihood function.  As logarithms are strictly increasing functions, the resulting solution from maximizing the likelihood vs. the log-likelihood are the equivalent.  Given a dataset of cardinality $N$, the log-likelihood (normalized by $N$) is given by:

$l = \frac{1}{N}\sum_{i=1}^{N}\Big(y_{i}\log(p) + (1 - y_{i})\log(1 - p)\Big)$

## Logistic function

Logistic regression performs binary classification based on a probabilistic interpretation of the data.  Essentially, the process seeks to assign a probability to new observations.  If the probability associated with the new instance of data is greater than 0.5, then the new observation is assigned to 1 (for example).  If the probability associated with the new instance of the data is less than 0.5, then it is assigned to 0.  To map the real numerical values into probabilities (which must lie between 0 and 1), logistic regression makes use of the logistic (sigmoid) function, given by:

$\sigma(t) = \frac{1}{1 + e^{-t}}$

Note that by setting $t=0$, $\sigma(0) = 0.5$, which is the decision boundary.  We should also note that the derivative of the logistic function with respect to the parameter $t$ is:

$\frac{d}{dt}\sigma(t) = \sigma(t)(1 - \sigma(t))$

## Logistic Regression and Derivation of the Gradient

Let's assume the training data consists of $N$ observations, where observation $i$ is denoted by the pair $(y_{i},\mathbf{x}_{i})$, where $y \in \{0,1\}$ is the label for the feature vector $\mathbf{x}$.  We wish to compute a linear decision boundary that best seperates the labeled observations.  Let $\mathbf{\theta}$ denote the vector of coefficients to be estimated.  In this problem, the log likelihood can be expressed as:

$l =  \frac{1}{N}\sum_{i=1}^{N}\Big(y_{i}\log\big(\sigma(\mathbf{\theta}^{T}\mathbf{x}_{i})\big) + (1 - y_{i}) \log\big( 1 - \sigma(\mathbf{\theta}^{T}\mathbf{x}_{i})\big)\Big)$

The gradient of the objective with respect to the $j^{th}$ element of $\mathbf{\theta}$ is:
$$
\begin{aligned}
\frac{d}{d\theta^{(j)}} l &=  \frac{1}{N}\sum_{i=1}^{N}\Bigg( \frac{d}{d\theta^{(j)}} y_{i}\log\big(\sigma(\mathbf{\theta}^{T}\mathbf{x}_{i})\big) + \frac{d}{d\theta^{(j)}}(1 - y_{i}) \log\big( 1 - \sigma(\mathbf{\theta}^{T}\mathbf{x}_{i})\big)\Bigg) \\
&= \frac{1}{N}\sum_{i=1}^{N}\Bigg(\frac{y_{i}}{\sigma(\mathbf{\theta}^{T}\mathbf{x}_{i})} - \frac{1 - y_{i}}{1 - \sigma(\mathbf{\theta}^{T}\mathbf{x}_{i})} \Bigg)\frac{d}{d\theta^{(j)}}\sigma(\mathbf{\theta}^{T}\mathbf{x}_{i})\\
&= \frac{1}{N}\sum_{i=1}^{N}\Bigg(\frac{y_{i}}{\sigma(\mathbf{\theta}^{T}\mathbf{x}_{i})} - \frac{1 - y_{i}}{1 - \sigma(\mathbf{\theta}^{T}\mathbf{x}_{i})} \Bigg)\sigma(\mathbf{\theta}^{T}\mathbf{x}_{i})\Big(1 - \sigma(\mathbf{\theta}^{T}\mathbf{x}_{i})\Big)x_{i}^{(j)}\\
&=  \frac{1}{N}\sum_{i=1}^{N}\Bigg(\frac{y_{i} - \sigma(\mathbf{\theta}^{T}\mathbf{x}_{i})}{\sigma(\mathbf{\theta}^{T}\mathbf{x}_{i})\Big(1 - \sigma(\mathbf{\theta}^{T}\mathbf{x}_{i})\Big)}\Bigg)\sigma(\mathbf{\theta}^{T}\mathbf{x}_{i})\Big(1 - \sigma(\mathbf{\theta}^{T}\mathbf{x}_{i})\Big)x_{i}^{(j)}\\
&= \frac{1}{N}\sum_{i=1}^{N}\Bigg(y_{i} - \sigma(\mathbf{\theta}^{T}\mathbf{x}_{i})\Bigg)x_{i}^{(j)}
\end{aligned}
$$

where the last equation has the familiar form of the product of the prediciton error and the $j^{th}$ feature.  With the gradient of the log likelihood function, the parameter vector $\mathbf{\theta}$ can now be calculated via gradient ascent (as we're <em>maximizing</em> the log likelihood):

$$
\begin{equation}
    \mathbf{\theta}^{(j)}(k+1) = \mathbf{\theta}^{(j)}(k) + \alpha \frac{1}{N}\sum_{i=1}^{N}\Bigg( y_{i} - \sigma(\mathbf{\theta}^{T}(k)\mathbf{x}_{i}))\Bigg)x_{i}^{(j)}
\end{equation}
$$


```python
# Supporting Methods

#logistic function
def sigmoid(a):
    return 1/(1 + np.exp(-a))

#ll function
def log_likelihood(x, y, theta):
    logits = np.dot(x, theta)
    log_like = np.sum(y * logits - np.log(1 + np.exp(logits)))
    return log_like
```


```python
#Load the data
iris = datasets.load_iris()
x = iris.data[:,2:]                    #features will be petal width and petal length
y = (iris.target==2).astype(np.int).reshape(len(x),1)    #1 of iris-virginica, and 0 ow

#Prepare Data for Regression
#pad x with a vector of ones for computation of intercept
x_aug = np.concatenate( (x,np.ones((len(x),1))) , axis=1)
```


```python
#sklearn logistic regression
log_reg = LogisticRegression(penalty='none')
log_reg.fit(x,y)
log_reg.get_params()
coefs = log_reg.coef_.reshape(-1,1)
intercept = log_reg.intercept_
theta_sklearn = np.concatenate((coefs, intercept.reshape(-1,1)), axis=0)
print("sklearn coefficients:")
print(theta_sklearn)
print("sklearn log likelihood: ", log_likelihood(x_aug, y, theta_sklearn))
```

    sklearn coefficients:
    [[  5.75452053]
     [ 10.44681116]
     [-45.27248307]]
    sklearn log likelihood:  -10.281754052558687


    C:\Users\danie\Anaconda3\Lib\site-packages\sklearn\utils\validation.py:63: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().
      return f(*args, **kwargs)



```python
#Perform Logistic Regression
num_iterations = 40000
alpha = 1e-2

theta0 = np.ones((3,1))
theta = []
theta.append(theta0)

k=0
while k < num_iterations:
    
    #compute prediction error
    e = y - sigmoid(np.dot(x_aug, theta[k]))

    #compute the gradient of the log-likelihood
    grad_ll = np.dot(x_aug.T, e)

    #gradient ascent
    theta.append(theta[k] + alpha * grad_ll)

    #update iteration step
    k += 1

    if k % 4000 == 0:
        #print("iteration: ", k, " delta: ", delta)
        print("iteration: ", k, " log_likelihood:", log_likelihood(x_aug, y, theta[k]))
    
theta_final = theta[k]
print("scratch coefficients:")
print(theta_final)
```

    iteration:  4000  log_likelihood: -11.529302688612685
    iteration:  8000  log_likelihood: -10.800986140073512
    iteration:  12000  log_likelihood: -10.543197464480874
    iteration:  16000  log_likelihood: -10.425775111214602
    iteration:  20000  log_likelihood: -10.36535749825992
    iteration:  24000  log_likelihood: -10.331961877189842
    iteration:  28000  log_likelihood: -10.312622172168293
    iteration:  32000  log_likelihood: -10.301055609225223
    iteration:  36000  log_likelihood: -10.293975480174714
    iteration:  40000  log_likelihood: -10.289566343598294
    scratch coefficients:
    [[  5.5044475 ]
     [ 10.17562424]
     [-43.60242548]]



```python
#Plot the data and the decision boundary

#create feature data for plotting
x_dec_bnd = np.linspace(0,7,100).reshape(-1,1)
#classification boundary from sklearn
y_sklearn = (theta_sklearn[2] * np.ones((100,1)) + theta_sklearn[0] * x_dec_bnd) / -theta_sklearn[1]
#classification boundary from scratch
y_scratch = (theta_final[2] * np.ones((100,1)) + theta_final[0] * x_dec_bnd) / -theta_final[1]

y_1 = np.where(y==1)[0] #training data, iris-virginica
y_0 = np.where(y==0)[0] #training data, not iris-virginica
plt.plot(x[y_0,0],x[y_0,1],'bo',label="not iris-virginica")
plt.plot(x[y_1,0],x[y_1,1],'k+',label="iris-virginica")
plt.plot(x_dec_bnd,y_sklearn,'r',label="sklearn dec. boundary")
plt.plot(x_dec_bnd,y_scratch,'g',label="scratch dec. boundary")
plt.xlabel('petal length')
plt.ylabel('petal width')
plt.title('Logistic Regression Classification')
plt.xlim((3,7))
plt.ylim((0,3.5))
plt.legend()
plt.show()
```

9fafb598db9b786adf895a5ebc4b7bfb11eb67f8
ipynb
Jupyter Notebook
.ipynb_checkpoints/Logistic_Regression-checkpoint.ipynb
dbarnold0220/from_scratch
088b059377fd735326daaff6fc924123a13f8324
MIT
.ipynb_checkpoints/Logistic_Regression-checkpoint.ipynb
dbarnold0220/from_scratch
088b059377fd735326daaff6fc924123a13f8324
MIT
.ipynb_checkpoints/Logistic_Regression-checkpoint.ipynb
dbarnold0220/from_scratch
088b059377fd735326daaff6fc924123a13f8324
MIT
2021-11-01T02:23:56.000Z
2021-11-01T02:23:56.000Z
Qwen/Qwen-72B
1. YES
2. YES