# Fundamentals of Data Science
Winter Semester 2020

## Info:

* Betti Martina (1799160)
* Ramse Gaurav (1965564)
* Ruà Giacomo (1770210)
* Sferragatta Stefania (1948081)
* Yildiz Mert (1951070)


## Prof. Fabio Galasso, Alessandro Flaborea and Luca Franco
<galasso@di.uniroma1.it>, <flaborea.1841969@studenti.uniroma1.it>, <franco.luca.17@gmail.com>

## Exercise 2

In Exercise 2, you will re-derive and implement logistic regression and optimize the parameters with Gradient Descent and with the Newton's method. Also, in this exercise you will re-derive and implement Gassian Discriminant Analysis.
We will use the files logistic_x.txt and logistic_y.txt. The first contains the feature values $x^{(i)}_1$ and $x^{(i)}_2$ for the $i$-th data sample $x^{(i)}$. The second contains the ground truth label $y^{(i)}$ for each corresponding data sample.

Note that, differently from what seen in lecture, $y^{(i)}\in \{-1,1\}$.

The completed exercise should be handed in as a single notebook file. Use Markdown to provide equations. Use the code sections to provide your scripts and the corresponding plots.
Submit it by sending an email to galasso@di.uniroma1.it, flaborea.1841969@studenti.uniroma1.it and franco.luca.17@gmail.com by Wednesday November 25th, 23:59.

## Notation

- $x^i$ is the $i^{th}$ feature vector
- $y^i$ is the expected outcome for the $i^{th}$ training example
- $m$ is the number of training examples
- $n$ is the number of features

Let's start by setting up our Python environment and importing the required libraries:


## [15 points] Question 1: Logistic Regression with Gradient Ascent

### (a) [5 points] Equations for the log likelihood, its gradient, and the gradient ascent update rule.

Write and simplify the likelihood $L(\theta)$ and log-likelihood $l(\theta)$ of the parameters $\theta$.

Recall the probabilistic interpretation of the hypothesis $h_\theta(x)= P(y=1|x;\theta)$ and that $h_\theta(x)=\frac{1}{1+\exp(-\theta^T x)}$.

Also derive the gradient $\frac{\delta l(\theta)}{\delta \theta_j}$ of $l(\theta)$ and write the gradient update equation. 

Question: in order to learn the optimal $\theta$, do you need to minimize or to maximize $l(\theta)$?

################# Do not write above this line #################

#### Likelihood:
$L(\theta)$ = $P(\vec{y}|x;\theta)$


=$\prod_{i=1}^{m}P(y^{(i)}|x^{(i)};\theta)$ 
            
            
=$\prod_{i=1}^{m}\frac{1}{\sqrt{2\pi}\sigma}\exp -\frac{(y^{(i)}-\theta^{T}x^{(i)})^{2}}{2\sigma^{2}}$


#### Log - Likelihood:

$l(\theta)$ = $log L(\theta)$ 

= $log \prod_{i=1}^{m}\frac{1}{\sqrt{2\pi}\sigma}\exp -\frac{(y^{(i)}-\theta^{T}x^{(i)})^{2}}{2\sigma^{2}}$

= $\sum_{i=1}^{m} [log \frac{1}{\sqrt{2\pi}\sigma} + log \exp -\frac{(y^{(i)}-\theta^{T}x^{(i)})^{2}}{2\sigma^{2}}] $

= $ m  log \frac{1}{\sqrt{2\pi}\sigma} + \sum_{i=1}^{m} - \frac{(y^{(i)}-\theta^{T}x^{(i)})^{2}}{2\sigma^{2}}$


#### Interpretation of the hypothesis $h_\theta(x)= P(y=1|x;\theta)$ :
##### Likelihood:

$L(\theta)$ = $P(\vec{y}|x;\theta)$  (according to the Bernoulli distribution)

= $ \prod_{i=1}^{m}P(y^{(i)}|x^{(i)};\theta)$ 

= $ \prod_{i=1}^{m} h_\theta(x^{(i)})^{y^{(i)}} (1-h_\theta(x^{(i)}))^{1-y^{(i)}}$


##### Log - likelihood:
$l(\theta)$ = $log L(\theta)$ 

= $\sum_{i=1}^{m} y^{(i)} log (h_\theta(x^{(i)}))+(1-y^{(i)}) log (1-h_\theta(x^{(i)})) $

##### Derive the gradient $\frac{\delta l(\theta)}{\delta \theta_j}$ of $l(\theta)$ :
$\theta_{j}$ = $\theta_{j} + \alpha \sum_{i=1}^{m} (y^{(i)}-h_\theta(x^{(i)}))x_{j}^{(i)} $

= $\theta_{j} + \alpha \frac{\delta l(\theta)}{\delta \theta_j}$



#### Question : minimize 
Our main goal is to maximize the value of the likelihood which is the function that measures the goodness of fit of a statistical model. Since the logarithmic function is a monotone one, in order to pick the best combination of model parameter values, we are looking for the maximum of the log likelihood function. 

################# Do not write below this line #################

### (b) [10 points] Implementation of logistic regression with Gradient Ascent

Code up the equations above to learn the logistic regression parameters. Consider as data the files logistic_x.txt and logistic_y.txt. The first contains the feature values $x^{(i)}_1$ and $x^{(i)}_2$ for the $i$-th data sample $x^{(i)}$. The second contains the ground truth label $y^{(i)}$ for each corresponding data sample.


```python
%matplotlib inline
import numpy as np # imports a fast numerical programming library
import scipy as sp # imports stats functions, amongst other things
import matplotlib as mpl # this actually imports matplotlib
import matplotlib.cm as cm # allows us easy access to colormaps
import matplotlib.pyplot as plt # sets up plotting under plt
import pandas as pd # lets us handle data as dataframes


# sets up pandas table display
pd.set_option('display.width', 500)
pd.set_option('display.max_columns', 100)
pd.set_option('display.notebook_repr_html', True)

import seaborn as sns # sets up styles and gives us more plotting options
```

Let's start by loading the data into two pandas DataFrames:


```python
df_x = pd.read_csv("./data/logistic_x.txt", sep="\ +", names=["x1","x2"], header=None, engine='python')
df_y = pd.read_csv('./data/logistic_y.txt', sep='\ +', names=["y"], header=None, engine='python')
df_y = df_y.astype(int)
df_x.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>x1</th>
      <th>x2</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>1.343250</td>
      <td>-1.331148</td>
    </tr>
    <tr>
      <th>1</th>
      <td>1.820553</td>
      <td>-0.634668</td>
    </tr>
    <tr>
      <th>2</th>
      <td>0.986321</td>
      <td>-1.888576</td>
    </tr>
    <tr>
      <th>3</th>
      <td>1.944373</td>
      <td>-1.635452</td>
    </tr>
    <tr>
      <th>4</th>
      <td>0.976734</td>
      <td>-1.353315</td>
    </tr>
  </tbody>
</table>
</div>




```python
df_y.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>y</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>-1</td>
    </tr>
    <tr>
      <th>1</th>
      <td>-1</td>
    </tr>
    <tr>
      <th>2</th>
      <td>-1</td>
    </tr>
    <tr>
      <th>3</th>
      <td>-1</td>
    </tr>
    <tr>
      <th>4</th>
      <td>-1</td>
    </tr>
  </tbody>
</table>
</div>



Let's get the NumPy arrays from the dataset, and add a column of 1's to $X$ to take into account the zero intercept


```python
x = np.hstack([np.ones((df_x.shape[0], 1)), df_x[["x1","x2"]].values])
y = df_y["y"].values
```

Adjust y so that $y^{(i)}\in [0,1]$

################# Do not write above this line #################

Your code here


```python
y= (y - y.min()) / (y.max() - y.min())
np.array(y,dtype=int)
```




    array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
           0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
           1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
           1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])



################# Do not write below this line #################


```python
[x[:5,:],x[-5:,:]] # Plot the first and last 5 lines of x, now containing features x0 (constant=1), x1 and x2
```




    [array([[ 1.        ,  1.3432504 , -1.3311479 ],
            [ 1.        ,  1.8205529 , -0.6346681 ],
            [ 1.        ,  0.98632067, -1.8885762 ],
            [ 1.        ,  1.9443734 , -1.635452  ],
            [ 1.        ,  0.97673352, -1.3533151 ]]),
     array([[ 1.        ,  4.7748538 ,  0.09941521],
            [ 1.        ,  5.8274854 , -0.69005848],
            [ 1.        ,  2.2894737 ,  1.9707602 ],
            [ 1.        ,  2.494152  ,  1.4152047 ],
            [ 1.        ,  2.0847953 ,  1.3567251 ]])]




```python
[y[:5],y[-5:]] # Plot the first and last 5 lines of y
```




    [array([0., 0., 0., 0., 0.]), array([1., 1., 1., 1., 1.])]



Define the sigmoid function "sigmoid", the function to compute the gradient of the log likelihood  "grad_l" and the gradient ascent algorithm.

################# Do not write above this line #################

Your code here


```python
def sigmoid(x,theta):
    z = np.dot(x,theta)
    g = 1/(1 + np.exp(-z))
    #the sigmoid of x changes according to the value of x. In case  x = 0, g(x)= 1/2; in case x -> +inf g(x) tends to 1; 
    # in case x -> -inf g(x) tends to 0
    return g # return the sigmoid of x
```


```python
def log_likelihood(theta,x,y):
    h_theta = sigmoid(x,theta)
    log_l= sum(y*np.log(h_theta) + (1-y)*np.log(1-h_theta))/len(y)
    return log_l # return the log likehood of theta according to data x and label y

```


```python
def grad_l(theta, x, y):
    m = len(y)
    h_theta=sigmoid(x,theta)
    G = np.dot((y-h_theta),x)/m
    
    return G # return the gradient G of the log likelihood
```


```python
def gradient_ascent(theta,x,y,G,alpha,iterations):
    m = len(y)
    log_l_history = np.zeros(iterations)
    theta_history = np.zeros((iterations,3)) 
        
    for i in range(1,iterations):
        hypothesis = sigmoid(x, theta) 
        log_l_history[i] = log_likelihood(theta,x,y)
        theta = theta_history[i-1] + (alpha*G(theta,x,y))
        theta_history[i] = theta        
        
    return theta, log_l_history, theta_history #as well as two lists containing the log likelihood's and values of theta at all iterations
```

################# Do not write below this line #################

Check your grad_l implementation:
grad_l applied to the theta_test (defined below) should provide a value for log_l_test close to the target_value (defined below); in other words the error_test should be 0, up to machine error precision.


```python
target_value = -0.4250958770469834
theta_test=np.array([-2,1,2])

log_l_test  = log_likelihood(theta_test,x,y)
error_test=np.abs(log_l_test-target_value)

print("{:f}".format(error_test))
```

    0.000000


Let's now apply the function gradient_ascent and print the final theta as well as theta_history 


```python
# Initialize theta0
theta0 = np.zeros(x.shape[1])

# Run Gradient Ascent method
n_iter=1000
theta_final, log_l_history, theta_history = gradient_ascent(theta0,x,y,grad_l,alpha=0.5,iterations=n_iter)
```

Let's plot the log likelihood over iterations


```python
fig,ax = plt.subplots(num=2)

ax.set_ylabel('l(Theta)')
ax.set_xlabel('Iterations')
_=ax.plot(range(len(log_l_history)),log_l_history,'b.')
```

Plot the data and the decision boundary:


```python
df_x_pl=df_x.copy()
df_y_pl=df_y.copy()
df_x_pl.insert(0, "y", df_y)
df_x_pl["y"] = pd.to_numeric(df_x_pl["y"], downcast='signed')
df_x_pl.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>y</th>
      <th>x1</th>
      <th>x2</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>-1</td>
      <td>1.343250</td>
      <td>-1.331148</td>
    </tr>
    <tr>
      <th>1</th>
      <td>-1</td>
      <td>1.820553</td>
      <td>-0.634668</td>
    </tr>
    <tr>
      <th>2</th>
      <td>-1</td>
      <td>0.986321</td>
      <td>-1.888576</td>
    </tr>
    <tr>
      <th>3</th>
      <td>-1</td>
      <td>1.944373</td>
      <td>-1.635452</td>
    </tr>
    <tr>
      <th>4</th>
      <td>-1</td>
      <td>0.976734</td>
      <td>-1.353315</td>
    </tr>
  </tbody>
</table>
</div>




```python
# Generate vector to plot decision boundary
x1_vec = np.linspace(df_x_pl["x1"].min(),df_x_pl["x1"].max(),2);

# Plot raw data
sns.scatterplot(x="x1", y="x2", hue="y", data=df_x_pl);

# Plot decision boundary
plt.plot(x1_vec,(-x1_vec*theta_final[1]-theta_final[0])/theta_final[2], color="red");

# Save the theta_final value for later comparisons
theta_GA = theta_final.copy()
print(theta_GA)
```

################# Do not write above this line #################

Write now your considerations. Discuss in particular:
- what effects do the learning rate $\alpha$ and the number of iterations have on the computed $\theta$
- may you still achieve the same value if you decrease the learning rate $\alpha$ and increase the number of iterations n_iter up to an infinite number of them? Why?
- may you still achieve the same value if you increase the learning rate $\alpha$ substantially and decrease the number of iterations n_iter? Why?

Substantiate your answer with experiments and plots of the $l(\theta)$ and the decision boundary.


```python
alphas=[0.001,0.005,0.01,0.05,0.2,0.4,0.7,1]
iterations=range(10,1000,10)
log_l=[]
p=[]
for alpha in alphas:
    log_l=[]
    for i in iterations:
        theta_final, log_l_history, theta_history = gradient_ascent(theta0,x,y,grad_l,alpha,int(i))
        log_l.append(np.mean(log_l_history))
    p.append(log_l)
```


```python
plt.figure(figsize=(15,10))
palette=sns.color_palette("rocket_r",8)
g=sns.lineplot(data=p,palette=palette)
plt.legend(title='alpha', loc='lower right', labels=alphas)
plt.show(g)
```

According to our analysis we observe that the behaviour of the log-likelihood of theta is strictly dependent on both the number of iterations and the value of the learning rate. As we can see in the plot above for small values of alpha the maximum is approched slowly, consequently the steepness of the curve is minimal. So for small values of the learning rate we obtain a good result only after a large number of iterations. However, the log likelihood that we reach with larger values of alpha, is never obtained with alpha smaller than 0.2. \
In the case in which we have a small number of iterations and also a small value of alpha we risk to undershoot the maximum, since more iterations would be needed in order to reach the maximum. To conclude we can state that the learning rate needed to reach the maximum is inversely proportional to the number of iterations.

################# Do not write below this line #################


## [15 points] Question 2: Logistic Regression with the Newton's method

### (a) [5 points] Equations for the Hessian of the log likelihood

Derive the equation for the Hessian matrix $H$ of the log-likelihood $l(\theta)$ of the parameters $\theta$.

Recall that: $H_{i,j} = \frac{\delta^2 l(\theta)}{\delta \theta_i \delta \theta_j}$

Also write the parameter update rule for the Newton's method, as a function of the Hessian and gradient of $l(\theta)$.

Question: in order to learn the optimal $\theta$, do you need to minimize or to maximize $\nabla_{\theta} l$?

################# Do not write above this line #################

##### Hessian:

$H_{i,j} = \frac{\delta^2 l(\theta)}{\delta \theta_i \delta \theta_j}$


$z = \frac{1}{1 + \exp(- \theta^{T}x)}$
 
$H = -\frac{1}{m}\sum_{i=1}^{m}{(z(1-z))x^{(i)}x^{(i)^{T}}}$

##### Parameter Update Rule for Newton's Method $\frac{\delta l(\theta)}{\delta \theta_j}$ of $l(\theta)$ :
$\theta_{j + 1}$ = $\theta_{j} - H^{-1} G $


##### Question : Do we minimizing or maximizing $\nabla_{\theta} l$?

As we have calculate and store the $\nabla_{\theta} l$ values for each iteration, in order to learn the optimal value of $\theta$, we need to minimize the $\nabla_{\theta} l \,$ function. 

################# Do not write below this line #################

### (b) [10 points] Implementation of logistic regression with the Newton's method

Code up the equations above to learn the logistic regression parameters. Consider as data the files logistic_x.txt and logistic_y.txt. The first contains the feature values  $x^{(i)}_1$ and $x^{(i)}_2$  for the  𝑖 -th data sample  $x^{(i)}$ . The second contains the ground truth label  $y^{(i)}$  for each corresponding data sample.

Define the Hessian function hess_l and the Newton's method function "newton":

################# Do not write above this line #################

Your code here


```python
def hess_l(theta, x, y):
    z = sigmoid(x, theta)
    diag = np.diag(z * (1 - z))
    hess = -(x.T.dot(diag.dot(x))/len(y))
    return hess

def newton(theta0, x, y, G, H, eps):
   
    theta_list = []
    delta_l =  np.linalg.norm(grad_l(theta0,x, y))
    l = log_likelihood(theta0, x,y)
    log_l_history = []
    theta_history = []
    G_history = []
    
    while delta_l>eps:
        H = hess_l(theta0,x,y)
        G = grad_l(theta0,x,y)
        H_inv = np.linalg.inv(H) 
        theta_history.append(theta0)
        theta0 = theta0 - H_inv.dot(G)
        log_l_history.append(log_likelihood(theta0,x,y))
        delta_l = np.linalg.norm((grad_l(theta0,x,y))) 
   
    theta_history= np.array(theta_history)
    log_l_history= np.array(log_l_history)
    # return the optimized theta parameters,
    # as well as two lists containing the log likelihood's and values of theta at all iterations
    return theta0, theta_history, log_l_history
```

################# Do not write below this line #################

Check your hess_l implementation: hess_l applied to the theta_test (defined below) should provide a value for hess_l_test close to the target_value (defined below); in other words the error_test should be 0, up to machine error precision.


```python
target_value = np.array([[-0.14357024, -0.43836134,  0.18767787], \
                         [-0.43836134, -1.72295801,  0.61817001], \
                         [ 0.18767787,  0.61817001, -0.51247901]])
theta_test=np.array([-2,1,2])

hess_l_test  = hess_l(theta_test,x,y)
error_test=np.sum(np.abs(hess_l_test-target_value),axis=(0,1))

print("{:f}".format(error_test))
```

    2.879550


Now run Newton's method to perform logistic regression:


```python
# Initialize theta0
theta0 = np.zeros(x.shape[1])


# Run Newton's method
theta_final, theta_history, log_l_history = newton(theta0,x,y,grad_l,hess_l,1e-6)
print(theta_final)
```

    [-2.62050954  0.76037096  1.17194549]


Let's plot the log likelihood over iterations


```python
fig,ax = plt.subplots(num=2)

ax.set_ylabel('l(Theta)')
ax.set_xlabel('Iterations')
_=ax.plot(range(len(log_l_history)),log_l_history,'b.')
```

Plot the data and our decision boundary:


```python
df_x_pl=df_x.copy()
df_y_pl=df_y.copy()
df_x_pl.insert(0, "y", df_y)
df_x_pl["y"] = pd.to_numeric(df_x_pl["y"],downcast='signed')
df_x_pl.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>y</th>
      <th>x1</th>
      <th>x2</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>-1</td>
      <td>1.343250</td>
      <td>-1.331148</td>
    </tr>
    <tr>
      <th>1</th>
      <td>-1</td>
      <td>1.820553</td>
      <td>-0.634668</td>
    </tr>
    <tr>
      <th>2</th>
      <td>-1</td>
      <td>0.986321</td>
      <td>-1.888576</td>
    </tr>
    <tr>
      <th>3</th>
      <td>-1</td>
      <td>1.944373</td>
      <td>-1.635452</td>
    </tr>
    <tr>
      <th>4</th>
      <td>-1</td>
      <td>0.976734</td>
      <td>-1.353315</td>
    </tr>
  </tbody>
</table>
</div>




```python
# Generate vector to plot decision boundary
x1_vec = np.linspace(df_x_pl["x1"].min(),df_x_pl["x1"].max(),2);

# Plot raw data
sns.scatterplot(x="x1", y="x2", hue="y", data=df_x_pl);

# Plot decision boundary
plt.plot(x1_vec,(-x1_vec*theta_final[1]-theta_final[0])/theta_final[2], color="red");
```

################# Do not write above this line #################

Write now your considerations. Discuss in particular:
- the different convergence speed of the Newton's method, compared to Gradient Ascent
- the incurred cost for the different convergence speed

In comparing the efficiency of Newton's method and the gradient ascent to find the optimum theta in the logistic regression model, the runtime complexity for each iteration in Newton's method is higher than the runtime complexity for each iteration in gradient descent. The memory complexity for each iteration in gradient ascent is $O(p)$, while the memory complexity for each iteration in Newton's method is $O(p^2)$. It can be seen, that Newton's method requires more number of operations and more amount of memory than gradient ascent at each iteration. However, Newton's method requires much less iterations to converge to the minimum because it uses the knowledge of the second derivative at each iteration, while gradient ascent only uses the knowledge of the first derivative at each iteration. This is a known trade-off between complexity and the speed of convergence of Newton's method and gradient ascent. Newton's method is very fast to converge yet requires a high amount of complexity for computing the inverse of the Hessian matrix.

**Note**: As we have dig a bit more to understand the logic of the methods which can be used in logistic regression we have found out that, according to Thomas (2003), the fastest algorithms turn out to be conjugate gradient ascent and *quasi-Newton* algorithms, which far outstrip Iterative Scaling and its variants. It inspired us to use these methods as a future work.

################# Do not write below this line #################


## [7 points] Question 3: Logistic Regression with non linear boundaries

### (a) Polynomial features for logistic regression

Define new features, e.g. of 2nd and 3rd degrees, and learn a logistic regression classifier by using the new features, by using the Newton's optimization algorithm (or the gradient ascent one, if you like).

In particular, we would consider a polynomial boundary with equation:

$f(x_1, x_2) = c_0 + c_1 x_1 + c_2 x_2 + c_3 x_1^2 + c_4 x_2^2 + c_5 x_1 x_2 + c_6 x_1^3 + c_7 x_2^3 + c_8 x_1^2 x_2 + c_9 x_1 x_2^2$

We would therefore compute 7 new features: 3 new ones for the quadratic terms and 4 new ones for the cubic terms.

Create new arrays by stacking x (from the previous questions) and the new 7 features (in the order x1x1, x2x2, x1x2, x1x1x1, x2x2x2, x1x1x2, x1x2x2). In particular create x_new_quad by additionally stacking with x the quadratic features, and x_new_cubic by additionally stacking with x the quadratic and the cubic features.


```python
# First extract features x1 and x2 from x and reshape them to x1 vector arrays
x1 = x[:,1]
x2 = x[:,2]
x1 = x1.reshape(x1.shape[0], 1)
x2 = x2.reshape(x2.shape[0], 1)
print(x[:5,:]) # For visualization of the first 5 values
print(x1[:5,:]) # For visualization of the first 5 values
print(x2[:5,:]) # For visualization of the first 5 values
```

    [[ 1.          1.3432504  -1.3311479 ]
     [ 1.          1.8205529  -0.6346681 ]
     [ 1.          0.98632067 -1.8885762 ]
     [ 1.          1.9443734  -1.635452  ]
     [ 1.          0.97673352 -1.3533151 ]]
    [[1.3432504 ]
     [1.8205529 ]
     [0.98632067]
     [1.9443734 ]
     [0.97673352]]
    [[-1.3311479]
     [-0.6346681]
     [-1.8885762]
     [-1.635452 ]
     [-1.3533151]]


################# Do not write above this line #################

Your code here


```python
# Define features: x1x1, x2x2, x1x2, x1x1x1, x2x2x2, x1x1x2, x1x2x2

# first intialize the variables ...
x1x1 = [];  x2x2 = []; x1x2= []; x1x1x1= []; x2x2x2= []; x1x1x2= []; x1x2x2= [];


# ... then define them 
for i in range(0, len(x1)):
    x1x1.append(np.square(x1[i]))
    x2x2.append(np.square(x2[i]))
    x1x2.append(x1[i]*x2[i])
    
    x1x1x1.append(np.power(x1[i], 3))
    x2x2x2.append(np.power(x2[i], 3))
    x1x1x2.append(x1[i]*x1[i]*x2[i])
    x1x2x2.append(x1[i]*x2[i]*x2[i])

# Stack the features together together with x, to define
# x_new_quad [x, x1x1, x2x2, x1x2]
# and x_new_cubic [x, x1x1, x2x2, x1x2, x1x1x1, x2x2x2, x1x1x2, x1x2x2]

x_new_quad = np.hstack ((x, x1x1, x2x2, x1x2))

x_new_cubic = np.hstack((x, x1x1, x2x2, x1x2, x1x1x1, x2x2x2, x1x1x2, x1x2x2))
```

################# Do not write below this line #################

Now use the Newton's optimization algorithm to learn theta by maximizing the log-likelihood, both for the case of x_new_quad and x_new_cubic.

Note: you may alternatively use gradient ascent, as in Question 1, if you like.


```python
# Initialize theta0, in case of quadratic features
theta0_quad = np.zeros(x_new_quad.shape[1])

# Run Newton's method, in case of quadratic features
theta_final_quad, theta_history_quad, log_l_history_quad = newton(theta0_quad,x_new_quad,y,grad_l,hess_l,1e-6)

# Initialize theta0, in case of quadratic and cubic features
theta0_cubic = np.zeros(x_new_cubic.shape[1])

# Run Newton's method, in case of quadratic and cubic features
theta_final_cubic, theta_history_cubic, log_l_history_cubic = newton(theta0_cubic,x_new_cubic,y,grad_l,hess_l,1e-6)

# check and compare with previous results
print(theta_final_quad)
print(theta_final_cubic)
```

    [-2.01513398  0.17762534  1.91642971  0.07051926  0.18229119 -0.1295329 ]
    [-15.80125013  11.84659121   8.30909256  -3.03962928   1.86524086
      -2.79099534   0.25914767   0.14781713   0.24456441  -0.30583232]


    <ipython-input-9-0d66dff45168>:3: RuntimeWarning: divide by zero encountered in log
      log_l= sum(y*np.log(h_theta) + (1-y)*np.log(1-h_theta))/len(y)
    <ipython-input-9-0d66dff45168>:3: RuntimeWarning: invalid value encountered in multiply
      log_l= sum(y*np.log(h_theta) + (1-y)*np.log(1-h_theta))/len(y)



```python
# Plot the log likelihood values in the optimization iterations, in one of the two cases.
fig,ax = plt.subplots(num=2)

ax.set_ylabel('l(Theta)')
ax.set_xlabel('Iterations')
_=ax.plot(range(len(log_l_history_cubic)),log_l_history_cubic,'b.')
```

### (b) Plot the computed non-linear boundary and discuss the questions



First, define a boundary_function to compute the boundary equation for the input feature vectors $x_1$ and $x_2$, according to estimated parameters theta, both in the case of quadratic (theta_final_quad) and of quadratic and cubic features (theta_final_cubic). Refer for the equation to the introductory part of Question 3.

################# Do not write above this line #################

Your code here


```python
def boundary_function(x1_vec, x2_vec, theta_final):
    
    x1_vec, x2_vec = np.meshgrid(x1_vec,x2_vec)
    
    if len(theta_final) == 6:
        # boundary function value for features up to quadratic
        c_0, c_1, c_2, c_3, c_4, c_5 = theta_final
        
        f = c_0 + c_1*x1_vec + c_2*x2_vec + c_3*np.square(x1_vec) +\
        c_4*np.square(x2_vec)+c_5*np.multiply(x1_vec, x2_vec)
    
    elif len(theta_final) == 10:
        # boundary function value for features up to cubic
        c_0, c_1, c_2, c_3, c_4, c_5, c_6, c_7, c_8, c_9 = theta_final
        
        f = c_0 + c_1*x1_vec + c_2*x2_vec + c_3*np.square(x1_vec) +\
        c_4*np.square(x2_vec)+c_5*np.multiply(x1_vec, x2_vec) +\
        c_6*np.power(x1_vec,3) + c_7*np.power(x2_vec,3) +\
        c_8*np.multiply(np.square(x1_vec), x2_vec) +\
        c_9*np.multiply(x1_vec, np.square(x2_vec))
     
    else:
        raise("Number of Parameters is not correct")
        
    return x1_vec, x2_vec, f
```

################# Do not write below this line #################

Now plot the decision boundaries corresponding to the theta_final_quad and theta_final_cubic solutions.


```python
x1_vec = np.linspace(df_x_pl["x1"].min()-1,df_x_pl["x1"].max()+1,200);
x2_vec = np.linspace(df_x_pl["x2"].min()-1,df_x_pl["x2"].max()+1,200);

x1_vec, x2_vec, f = boundary_function(x1_vec, x2_vec, theta_final_quad)

sns.scatterplot(x="x1", y="x2", hue="y", data=df_x_pl);

plt.contour(x1_vec, x2_vec, f, colors="red", levels=[0])
```


```python
x1_vec = np.linspace(df_x_pl["x1"].min()-1,df_x_pl["x1"].max()+1,200);
x2_vec = np.linspace(df_x_pl["x2"].min()-1,df_x_pl["x2"].max()+1,200);

x1_vec, x2_vec, f = boundary_function(x1_vec, x2_vec, theta_final_cubic)

sns.scatterplot(x="x1", y="x2", hue="y", data=df_x_pl);

plt.contour(x1_vec, x2_vec, f, colors="red", levels=[0])
```

################# Do not write above this line #################

Write now your considerations. Discuss in particular:
- What happens if you further increase the degree of the polynomial. How would the boundary change? Would you incur underfitting or overfitting?
- Look at the boundary of the quadratic and cubic polynomial features. Do they meet your expectations? Why? Consider in particular the classifier decision boundary in areas where there are fewer training data, i.e. where the classifier is supposed to generalize. Hint: Try changing the minimum and maximum limits in vectors x1_vec and x2_vec before plotting the boundaries.
- Look at the values of the original features, the 2nd and 3rd degree ones. Do you expect that normalizing the input features (either by setting their range to [-1,1] or by subtracting their mean and dividing by their standard deviation) would improve the classifier? Why yes or why not?

**Considerations:**

* In general when we increase the degree of the polynomial boundary, the boundary becomes more flexible. Therefore it will fit more closely the training points.
    For this reason if we increase the degree we would get a more winding red line, that will select more and more closely every training point. Consequently we might incur in overfitting. 

* For the boundary of the quadratic polynomial features we expect to plot a parabolic graph (and that's exactly what we got). The problem with this polynomial boundary (as with all quadratic boundaries) is that while we might get a good classification in areas where we have a high density of training data, as we move in areas where the data is more scattered the graph will decrease rapidly, causing an underfittitng.  

    Let's analyze this point more closely. Consider the quadratic polynomial boundary in the section of the plane $[5,8] \times  [-5, 1] $ : 


```python
x1_vec = np.linspace(5,8, 200)
x2_vec = np.linspace(-5, 1, 200)

x1_vec, x2_vec, f = boundary_function(x1_vec, x2_vec, theta_final_quad)

sns.scatterplot(x="x1", y="x2", hue="y", data=df_x_pl);
plt.contour(x1_vec, x2_vec, f, colors="red", levels=[0])
```

 From the graph above we can clearly see that the points on the right half of the plane will not be detected by the quadratic polynomial boundary. 
    One way to fix this problem is to add to the quadratic component a cubic one:


```python
x1_vec = np.linspace(5,8, 200)
x2_vec = np.linspace(-5, 1, 200)

x1_vec, x2_vec, f = boundary_function(x1_vec, x2_vec, theta_final_cubic)

sns.scatterplot(x="x1", y="x2", hue="y", data=df_x_pl);
plt.contour(x1_vec, x2_vec, f, colors="red", levels=[0])
```

* Let's explore the ranges of the training data, the features we have started with have pretty similar ranges:


```python
print(max(x1) - min(x1))
print(max(x2) - min(x2))
```

    [7.13460119]
    [7.8193692]


When we square them (or elevate them  to the third power), obviously the ranges start increasing rapidly. For example: 


```python
print(max(x2x2) - min(x2x2))
print(max(x1x1x1) - min(x1x1x1))
```

    [16.0517872]
    [457.30830535]


Therefore we end up with features at different scales. 
In general having features at similar scales is good practice, since we don't want the existence of some *privileged* axis, meaning that a unit increase along one axis counts more than the others.  

In particular in this code we are using the Newton's optimization algorithm to calculate `theta_final`: this algorithm uses the euclidean distance to calculate each iteration, meaning that the features with higher magnitude will weigh in more in the calculation of `theta_final`. 

################# Do not write below this line #################


## [10 points, extra, not compulsory] Question 4: Gaussian Discriminant Analysis



### (a) Review and implement classification with Gradient Discriminant Analysis (GDA)

Recall the assumptions for GDA

$\begin{aligned} y & \sim \operatorname{Bernoulli}(\phi) \\ x \mid y=0 & \sim \mathcal{N}\left(\mu_{0}, \Sigma\right) \\ x \mid y=1 & \sim \mathcal{N}\left(\mu_{1}, \Sigma\right) \end{aligned}$

Recall the probability distributions:

$\begin{aligned} p(y) &=\phi^{y}(1-\phi)^{1-y} \\ p(x \mid y=0) &=\frac{1}{(2 \pi)^{d / 2}|\Sigma|^{1 / 2}} \exp \left(-\frac{1}{2}\left(x-\mu_{0}\right)^{T} \Sigma^{-1}\left(x-\mu_{0}\right)\right) \\ p(x \mid y=1) &=\frac{1}{(2 \pi)^{d / 2}|\Sigma|^{1 / 2}} \exp \left(-\frac{1}{2}\left(x-\mu_{1}\right)^{T} \Sigma^{-1}\left(x-\mu_{1}\right)\right) \end{aligned}$

Overall, the parameters of the model are $\phi, \Sigma, \mu_{0}$ and $\mu_{1}$ and the log-likelihood is given by:

$\begin{aligned} \ell\left(\phi, \mu_{0}, \mu_{1}, \Sigma\right) &=\log \prod_{i=1}^{n} p\left(x^{(i)}, y^{(i)} ; \phi, \mu_{0}, \mu_{1}, \Sigma\right) \\ &=\log \prod_{i=1}^{n} p\left(x^{(i)} \mid y^{(i)} ; \mu_{0}, \mu_{1}, \Sigma\right) p\left(y^{(i)} ; \phi\right) \end{aligned}$

Parameters may be estimated by maximum likelihood estimate as follows:

$\begin{aligned} \phi &=\frac{1}{n} \sum_{i=1}^{n} 1\left\{y^{(i)}=1\right\} \\ \mu_{0} &=\frac{\sum_{i=1}^{n} 1\left\{y^{(i)}=0\right\} x^{(i)}}{\sum_{i=1}^{n} 1\left\{y^{(i)}=0\right\}} \\ \mu_{1} &=\frac{\sum_{i=1}^{n} 1\left\{y^{(i)}=1\right\} x^{(i)}}{\sum_{i=1}^{n} 1\left\{y^{(i)}=1\right\}} \\ \Sigma &=\frac{1}{n} \sum_{i=1}^{n}\left(x^{(i)}-\mu_{y^{(i)}}\right)\left(x^{(i)}-\mu_{y^{(i)}}\right)^{T} \end{aligned}$

Consider the data variables as from the previous questions, x and y.
For this question, we would select features $x^{(i)}_1$ and $x^{(i)}_2$ into the new varible x12, dropping the $x_0=1$ convention.


```python
# First extract the sole features x1 and x2 from x into an x12 array and reshape it to x2 vector arrays
x12 = x[:,1:]
print(x[:5,:]) # For visualization of the first 5 values
print(x12[:5,:]) # For visualization of the first 5 values
print(y[:5]) # For visualization of the first 5 values
```

    [[ 1.          1.3432504  -1.3311479 ]
     [ 1.          1.8205529  -0.6346681 ]
     [ 1.          0.98632067 -1.8885762 ]
     [ 1.          1.9443734  -1.635452  ]
     [ 1.          0.97673352 -1.3533151 ]]
    [[ 1.3432504  -1.3311479 ]
     [ 1.8205529  -0.6346681 ]
     [ 0.98632067 -1.8885762 ]
     [ 1.9443734  -1.635452  ]
     [ 0.97673352 -1.3533151 ]]
    [0. 0. 0. 0. 0.]


Now code up the equations above to learn the GDA parameters from the training data.


################# Do not write above this line #################


Your code here


```python
import math

def compute_phi(y):
    
    phi = 1/len(y)*sum(y)
    
    return phi


def compute_mu0(x, y):
    
    numerator=np.zeros((2,))    
    sum = list(y).count(0)
    
    for i in range (len(y)):
        if(y[i]==0):
            numerator+=x[i]
    mu0 = numerator/sum 
    
    return mu0


def compute_mu1(x, y):
    
    numerator=np.zeros((2,))
    sum = list(y).count(1)
    
    for i in range (len(y)):
        if(y[i]==1):
            numerator+=x[i]
    mu1=numerator/sum
    
    return mu1


def compute_sigma(x, y, mu0, mu1):
    mu0=mu0.reshape(2,1)
    mu1=mu1.reshape(2,1)
    sigma=np.zeros((2,2))
    for i in range (len(y)):
        if(y[i]==1):
            sigma+=(x12[i].reshape((1,2)).T-mu1).dot((x12[i].reshape((1,2)).T-mu1).T)
        else:
            sigma+=(x12[i].reshape((1,2)).T-mu0).dot((x12[i].reshape((1,2)).T-mu0).T)
    sigma = sigma/len(y)
    return sigma


# computing p(x|y) for the Bayes rule
def p_x_given_y(x,mu,sigma):
    
    d = 2
    det = np.linalg.det(sigma)
    det_s = det**0.5
    m_pi =(2*math.pi)**(d/2)
    c = m_pi*det_s
    d = 1/c
    
    a = ((x-mu).T).dot(np.linalg.inv(sigma))
    b = a.dot(x-mu)
    e = -(0.5*b)
    p_x_giv_y = d*np.exp(e)

    return p_x_giv_y
```


################# Do not write below this line #################



```python
# Support functions

# Function p(y) for applying the Bayes rule
def p_y(y,phi):
    if y==1: return phi
    else: return 1-phi
```


```python
# Now estimate the GDA parameters
phi = compute_phi(y)
mu0 = compute_mu0(x12,y)
mu1 = compute_mu1(x12,y)
sigma = compute_sigma(x12, y, mu0, mu1)

phi,mu0,mu1,sigma
```




    (0.494949494949495,
     array([ 2.93068144, -1.44476469]),
     array([4.93616668, 0.76771897]),
     array([[ 2.68174521, -0.46404053],
            [-0.46404053,  1.88367821]]))



Now that we know both the distributions of $p(y)$ and $p(x|y)$, we may estimate the posterior probability for an input $x$ via Bayes Rule:

\begin{equation}
p(y \mid x)=\frac{p(x \mid y) p(y)}{p(x)}
\end{equation}
where $p(x)$ is the same for both classes and may be disregarded.

This is used to classify each data point, by comparing $p(y=1|x) \propto p(x \mid y=1) p(y=1)$ and $p(y=0|x) \propto p(x \mid y=0) p(y=0)$.


```python
selected_example=1
print('Selected example =', x12[selected_example,:])

# compute p(y=0|x) ~ p(x|y=0)*p(y=0)  &  p(y=1|x) ~ p(x|y=1)*p(y=1)

# y=0
print('p(y=0|x) ~', p_x_given_y(x12[selected_example,:],mu0,sigma)*p_y(0,phi))

# y=1
print('p(y=1|x) ~', p_x_given_y(x12[selected_example,:],mu1,sigma)*p_y(1,phi))
```

    Selected example = [ 1.8205529 -0.6346681]
    p(y=0|x) ~ 0.026129116841520316
    p(y=1|x) ~ 0.0020617882761784934


### (b) Plot the estimated likelihood functions, the decision boundary and write your remarks



```python
def boundary_gda(x1_vec, x2_vec, mu0, mu1, sigma, phi):
    
    x1_vec, x2_vec = np.meshgrid(x1_vec,x2_vec)
    
    f=np.zeros(x1_vec.shape)
    
    for i in range(f.shape[0]):
        for j in range(f.shape[1]):
            x12_ij = np.array( (x1_vec[i,j],x2_vec[i,j]) )
            f[i,j] = p_x_given_y(x12_ij,mu1,sigma)*p_y(1,phi) - p_x_given_y(x12_ij,mu0,sigma)*p_y(0,phi)
            
    return x1_vec, x2_vec, f


x1_vec = np.linspace(df_x_pl["x1"].min()-1,df_x_pl["x1"].max()+1,200);
x2_vec = np.linspace(df_x_pl["x2"].min()-1,df_x_pl["x2"].max()+1,200);

x1_vec, x2_vec, f = boundary_gda(x1_vec, x2_vec, mu0, mu1, sigma, phi)
```


```python
from matplotlib import cm
from scipy.stats import multivariate_normal

# Define Gaussian pdfs
rv0 = multivariate_normal(mu0, sigma)
rv1 = multivariate_normal(mu1, sigma)

pairx12 = np.dstack((x1_vec, x2_vec))

# plot pdf contours
plt.contourf(x1_vec, x2_vec, rv0.pdf(pairx12),alpha=0.5,cmap=cm.Reds)
plt.contourf(x1_vec, x2_vec, rv1.pdf(pairx12),alpha=0.35,cmap=cm.Blues)

# Add the GDA classification to the plot
sns.scatterplot(x="x1", y="x2", hue="y", data=df_x_pl);
plt.contour(x1_vec, x2_vec, f, colors="red", levels=[0])

plt.ylim(df_x_pl["x2"].min()-1,df_x_pl["x2"].max()+1)
plt.xlim(df_x_pl["x1"].min()-1,df_x_pl["x1"].max()+1)
```


```python
# For comparison with the Logistic Regression boundary

# Generate vector to plot decision boundary
x1vec = np.linspace(df_x_pl["x1"].min()-1,df_x_pl["x1"].max()+1,2)

# Plot raw data
sns.scatterplot(x="x1", y="x2", hue="y", data=df_x_pl)

# Plot decision boundary
log, = plt.plot(x1vec,(-x1vec*theta_GA[1]-theta_GA[0])/theta_GA[2], color="red")
gda = plt.contour(x1_vec, x2_vec, f, colors="green", levels=[0])

plt.legend([log, gda.collections[0]],['logistic regression','GDA'])

plt.ylim(df_x_pl["x2"].min()-1,df_x_pl["x2"].max()+1)
plt.xlim(df_x_pl["x1"].min()-1,df_x_pl["x1"].max()+1)

# plt.legend()
plt.show()
```

Write now your considerations. Discuss in particular:
- Gaussian Discriminant Analysis makes stronger modelling assumptions than logistic regression. When may this be helpful and where may this be harmful?

#### Conclusion marks:

When the assumption of Gaussian Discriminant Analysis (GDA) hold, for example if $P(x|y)$ is multivariate gaussian with shared $Σ$, the GDA will perform better than any models, including logistic regression. \
On the other hand, logistic regression is more robust and less sensative to model assumptions.
In this model we have used the same covariance matrix Σ for both classes. This means that whatever my normal distribution looks like for one class - no matter how tall/fat/slanty it is - we assumed that the other class' covariance matrix looks exactly like that as well. \
This is a very strong assumption which can certainly be false when Σ varies for different classes and it might not be applicable for other cases.

Credits for material: Andrew Ng, C. Combier

0d5480d56289aa29c1cc11cef23063dc1a440a97
ipynb
Jupyter Notebook
FDS_Exercise_2_new_assigned.ipynb
MertYILDIZ19/Logistic_Regression_with_Gradient_Descent-Newton_Method_and_Gassian-Discriminant-Analysis
33e1315f3941e307d0baef3c1b9ce5654fac00d3
MIT
FDS_Exercise_2_new_assigned.ipynb
MertYILDIZ19/Logistic_Regression_with_Gradient_Descent-Newton_Method_and_Gassian-Discriminant-Analysis
33e1315f3941e307d0baef3c1b9ce5654fac00d3
MIT
FDS_Exercise_2_new_assigned.ipynb
MertYILDIZ19/Logistic_Regression_with_Gradient_Descent-Newton_Method_and_Gassian-Discriminant-Analysis
33e1315f3941e307d0baef3c1b9ce5654fac00d3
MIT
Qwen/Qwen-72B
1. YES
2. YES