# Linear Regression
## Unique variable

We have a dataset formed by $n$ values of $(x,y)$.

*Problem*

If we get some $x_{new}$ value could we estimate a correspondent $y_{new}$?

*Hypothesis 1*

We admit a linear relation between $x$ and $y$, i.e.,
\begin{equation}
\hat{y} = a x + b
\end{equation}
where $a$ is the linear coeficient or slope and $b$ is the independent term.

As we have any guarantee that our hypothesis holds, it's probable that our estimative $\hat{y}$ differs from the actual $y$. Therefore we have a potential associated error, which we will call residual $r$:
\begin{equation}
r = y - \hat{y} = y - (a x + b) = \sum_i^n r_i \implies r_i = y_i - ax_i -b_i
\end{equation}

Our estimative could be:
1. above the actual value, i.e., $y < \hat{y} \implies y - \hat{y} = r < 0$
2. below the actual value, i.e., $y > \hat{y} \implies y - \hat{y} = r > 0$
3. equal the actual value, i.e., $y = \hat{y} \implies y = \hat{y} = r = 0$

In order to be able to deal with all these cases, we square this residual:
\begin{equation}
r^2 = (y - \hat{y})^2 = (y - (a x + b))^2
\end{equation}

Let's improve the notation, treating each sample as $i$. And calculate the sum $S$ of the squared residuals:
\begin{equation}
S(a,b) = \sum_{i=1}^n r_i^2 = \sum_{i=1}(y_i - \hat{y}_i)^2 = \sum_{i=1}(y_i - a x_i - b)^2
\end{equation}
as our $n$ $(x,y)$ variables are known, we have a function of the coeficients $(a,b)$.

In order to minimize this sum $S$, we will partially derivate this expression towards each variable $a,b$; apply the chain rule and equalize this to zero.
\begin{align}
\frac{\partial S}{\partial a} &
= \frac{\partial S}{\partial r}\frac{\partial r}{\partial a} = 0\\
\frac{\partial S}{\partial b} & 
= \frac{\partial S}{\partial r}\frac{\partial r}{\partial b} = 0
\end{align}

Each term is calculated by means of:
\begin{align}
\frac{\partial S}{\partial r} &
= \frac{\mathrm{d} \sum_{i=1}^n  r_i^2}{\mathrm{d} r_i} = 2 \sum_{i=1}^n r_i
= 2\sum_{i=1}^n (y_i-ax_i-b)\\
\frac{\partial r_i}{\partial a} &= \frac{\partial (y_i-ax_i-b)}{\partial a} = -x_i\\
\frac{\partial r_i}{\partial b} &= \frac{\partial (y_i-ax_i-b)}{\partial b} = -1
\end{align}

We could replace as:
\begin{align}
\frac{\partial S}{\partial r}\frac{\partial r}{\partial a}
& = 2\sum_{i=1}^n(y_i-ax_i-b) (-x_i) = - 2 \sum_{i=1}^n x_i (y_i-ax_i-b)
= 0 \\
\frac{\partial S}{\partial r}\frac{\partial r}{\partial b}
& = 2 \sum_{i=1}^n (y_i-ax_i-b) (-1) 
= -2 \sum_{i=1}^n (y_i-ax_i-b)
= 0
\end{align}

If we evaluate the second expression:
\begin{align}
0 &= \sum_{i=1}^n (y_i-ax_i-b)
= \sum_{i=1}^n y_i - \sum_{i=1}^n a x_i - \sum_{i=1}^n b 
\end{align}

Dividing by $n$ samples we get mean values:
\begin{align}
0 = \frac{\sum_{i=1}^n y_i}{n} - \frac{\sum_{i=1}^n a x_i}{n} - \frac{\sum_{i=1}^n b}{n}
= \bar{y} - a \bar{x} - b
\implies b = \bar{y} - a \bar{x}
\end{align}

If we replace $b$ in the first derivative expression:
\begin{align}
0 &= \sum_{i=1}^n x_i (y_i-ax_i-b)
= \sum_{i=1}^n x_i (y_i-ax_i-(\bar{y} - a \bar{x})) \\
& = \sum_{i=1}^n x_i (y_i - ax_i - \bar{y} + a \bar{x})
= \sum_{i=1}^n x_i ( y_i - \bar{y}) - a \sum_{i=1}^n x_i ( x_i - \bar{x})
\end{align}

Isolating $a$
\begin{align}
a = \frac{\sum_{i=1}^n x_i ( y_i - \bar{y})}{\sum_{i=1}^n x_i ( x_i - \bar{x})}
\end{align}

In a nutshell, linear regression for one variable consist in admit a linear relation between $x$ and $y$, i.e., $\hat{y}=a+bx$, where the coeficients are determined by:
\begin{align}
a &= \frac{\sum_{i=1}^n x_i ( y_i - \bar{y})}{\sum_{i=1}^n x_i ( x_i - \bar{x})}\\
b &= \bar{y} - a \bar{x}
\end{align}
with the sample means:
\begin{align}
\bar{x} = \frac{\sum_{i=1}^n x_i}{n}, \qquad \bar{y} &= \frac{\sum_{i=1}^n y_i}{n}
\end{align}

## Multivariate

Hypothesis:

Admit a linear relation:
\begin{equation}
y = Xb + r
\end{equation}

Define the residual as:
\begin{equation}
r = y - Xb
\end{equation}

Compute the squared residual as:
\begin{align}
S(b) &= (y - X b)^T (y - X b)\\
&=  ( y^T - (X b)^T ) (y - X b)\\
&=  ( y^T - b^T X^T ) (y - X b)\\
&=  y^T y - b^T X^T y -  y^T X b + b^T X^T X b
\end{align}

The second and third term are equal, like we can see through:
\begin{equation}
b^T X^T y = (X b)^T y = y^T X b = X^T y b
\end{equation}

The last term could be rewritten as:
\begin{equation}
b^T X^T X b = (X^T X b)^T b = X^T X b b^T
\end{equation}

If we calculate the partial derivative and equals the result to zero i.e., (null residual):
\begin{align}
\frac{\partial S}{\partial b} 
& =  - 2 X^T y  \frac{\partial b^T}{\partial b}   
+ X^T X \frac{\partial b^Tb}{\partial b}\\
&= - 2 X^T y + 2 X^T X b = 0
\end{align}

Therefore we get the coefficient vector as:
\begin{align}
X^T y = X^T X b \implies b = (X^T X)^{-1} X^T y 
\end{align}

Notice the consistence as we recover our initial hypothesis for the case with a null residual through:
\begin{align}
b = (X^T X)^{-1} X^T y = X^{-1}X^{-T} X^T y = X^{-1} y
\end{align}

In a nutshell
\begin{equation}
y = Xb + r, \qquad b = (X^T X)^{-1} X^T y
\end{equation}

Observation:
Remember to add the independent variable as $ x^0=1 $:

\begin{equation}
\begin{Bmatrix}
y_1 \\
y_2 \\
\vdots \\
y_n
\end{Bmatrix}
=
\begin{bmatrix}
1 & x^{1}_{1} & x^{2}_{1} & \cdots & x^{m}_{1} \\
1 & x^{1}_{2} & x^{2}_{2} & \cdots & x^{m}_{2}  \\
\vdots & \vdots & \vdots & \ddots & \vdots\\
1 & x^{1}_{n} & x^{2}_{n} & \cdots & x^{m}_{n}
\end{bmatrix}
\times
\begin{Bmatrix}
b_0 \\
b_1 \\
\vdots \\
b_m
\end{Bmatrix}
+
\begin{Bmatrix}
r_1 \\
r_2 \\
\vdots \\
r_n
\end{Bmatrix}
\end{equation}

Hypothesis:
1. Fixed regressors (X not stochastic)
2. Aleatory error with null mean
3. Constant error variance (homoscedasticity)
4. No error correlation
5. b constant
6. Linear model
7. Error normal distributed


## $R^2$

$R^2$ determination coeficient
* quality of estimative
* $R^2$ of the $y$ variance explained by the $X$ variance
\begin{equation}
R^2 = 1 - \frac{SS_{res}}{SS_{tot}}, \qquad
SS_{res} = \sum_{i=1}^n r_i^2, \qquad
SS_{tot} = \sum_{i=1}^n (y_i - \bar{y})^2
\end{equation}
* $SS_{res}$: squared sum of the residuals $r$
* $SS_{tot}$: squared sum of the distance of each variable to $y_i$ to average $\bar{y}$

Estimative and residual
\begin{equation}
\hat{y} = Xb, \qquad
r_i = \hat{y}_i - y_i
\end{equation}


Estimator has a normal distribution
\begin{equation}
b\sim \mathcal{N}\left(\beta,\sigma^2(X^T X)^{-1}\right)
\end{equation}

Sample variance approximation
\begin{equation}
s^2 = \frac{SS_{res}}{n-(m+1)}
\end{equation}
* $n$: number of samples
* $m$: number of variables
* 1: $x_0$

t-Student-statistic
\begin{equation}
t_i = \frac{b_i}{s\sqrt{(X^T X)^{-1}_{ii}}}
\end{equation}

Null hypothesis ($H_0$)
* $\beta_i=0$
* the coeficients have no influence


If $|t_i| > t(1-\alpha) \implies H_0$ could be rejected with at least $(1-\alpha)$ confidence.

where $\alpha$ is the significance level




## Homemade implementation


```
# Multivariate Linear Regression
# Author: Vinícius Rios Fuck
# Date: 20/07/2020


import pandas as pd
import numpy as np
from scipy.stats import t

# Sample input
# Source: https://pt.wikipedia.org/wiki/M%C3%A9todo_dos_m%C3%ADnimos_quadrados

d = {"i": [1,2,3,4,5,6,7,8,9,10],
     "y": [122,114,86,134,146,107,68,117,71,98],
     "x1": [139,126,90,144,163,136,61,62,41,120],
     "x2": [0.115,0.120,0.105,0.090,0.100,0.120,0.105,0.080,0.100,0.115]
}
df = pd.DataFrame(data=d)

# Adjust input data

# select the y variable
y = df['y'].to_numpy()
# select the X variables
X = df.iloc[:,2:].to_numpy()
# add x_0 column of ones (independent coeficient)
X = np.c_[ np.ones(len(X)), X ]
```


```
def LR_predict(x, b):
  y_hat = x @ b
  return y_hat


def squared_sum(x):
  squared_x_sum = np.square(x).sum()
  return squared_x_sum


def residual_f(actual, predict):
  residual = actual - predict
  return residual


def R_2_score(squared_residuals_sum, squared_y_y_bar_distance_sum):
  R_2 = 1 - squared_residuals_sum / squared_y_y_bar_distance_sum
  return R_2


def adjusted_R_2_score(R_2, n, dof):
  adjusted_R_2 = 1 - (n-1)/dof * (1-R_2)
  return adjusted_R_2


# Standard Deviation
def std_sample(squared_residuals_sum, dof):
  sd_sample = (squared_residuals_sum / dof)**0.5
  return sd_sample


def t_stat(b, sd_sample, XtX_inv):
  t_student_coef = b / (sd_sample * np.diag(XtX_inv)**0.5)
  return t_student_coef

def significance(confidence=0.975):
  alpha = 1 - confidence
  return alpha


def critical_value_calc(dof, confidence=0.975):
  critical_value = t.ppf(confidence, dof)
  return critical_value


def round_3(input):
  return round(input, 3)


def exp_2(input):
  return '{:.2e}'.format(input)


def print_hypothesis(b, p_value, t_student_coef, critical_value, alpha):
  cv_out = round_3(critical_value)
  alpha_out = round_3(alpha)
  print("H0: null hypothesis: the coefficients bi are not relevant\n")

  accept_H0_critical_value = np.zeros(len(b))
  accept_H0_p_value = np.zeros(len(b))

  for i in range(len(b)): # all variables
    t_out = round_3(abs(t_student_coef[i]))
    # interpret via critical value
    accept_H0_critical_value[i] = (abs(t_student_coef[i]) <= critical_value)
    if accept_H0_critical_value[i]:
      print(f'Accept H0 that b{i} is not relevant.')
      print(f't_student = {t_out} <= critical_value = {cv_out}\n')
    else:
      print(f'Reject H0 that b{i} is not relevant.')
      print(f't_student = {t_out} > critical_value = {cv_out}\n')

    p_out = exp_2(p_value[i])
    # interpret via p-value
    accept_H0_p_value[i] = (p_value[i] > alpha)
    if accept_H0_p_value[i]:
      print(f'Accept H0 that b{i} is not relevant.')
      print(f'p_value = {p_out} > alpha = {alpha_out}\n')
    else:
      print(f'Reject H0 that b{i} is not relevant.')
      print(f'p_value = {p_out} <= alpha = {alpha_out}\n')

  return accept_H0_critical_value, accept_H0_p_value

```


```
# R^2 (quality of estimative)
def R2_SqResidSum(X, y, b):
  # y mean
  y_bar = y.mean()
  # predicted y
  y_hat = LR_predict(X, b=b)
  # residual
  residual = residual_f(y, y_hat)
  # squared residuals sum
  squared_residuals_sum = squared_sum(residual)
  # squared_y_y_bar_distance_sum
  squared_y_y_bar_distance_sum = squared_sum(y - y_bar)
  # R_2
  R_2 = R_2_score(squared_residuals_sum, squared_y_y_bar_distance_sum)
  return R_2, squared_residuals_sum
```


```
# t-Student stats (statistical relevance)
def t_student_stats(XtX_inv, b, squared_residuals_sum, R_2,
                    confidence=0.975, n=len(y)):
  # sample size
  # n = len(y)
  # number of variables # already take in account the independent coeficient
  # m = X.shape[1]
  m = len(b)
  # degree of freedom
  dof = n - m
  # adjusted_R_2
  adjusted_R_2 = adjusted_R_2_score(R_2, n, dof)
  # Standard Deviation
  sd_sample = std_sample(squared_residuals_sum, dof)
  # t_student_coef
  t_student_coef = t_stat(b, sd_sample, XtX_inv)
  alpha = significance(confidence)
  critical_value = critical_value_calc(dof, confidence)
  # calculate the p-value
  # 2: two-tailed distribution
  p_value = (1 - t.cdf(abs(t_student_coef), dof)) * 2
  (accept_H0_critical_value,
   accept_H0_p_value) = print_hypothesis(b, p_value, t_student_coef,
                                         critical_value, alpha)

  dict_H0 = {"t_student": t_student_coef,
             'p_value': p_value,
             'accept_H0_critical_value': accept_H0_critical_value,
             'accept_H0_p_value': accept_H0_p_value
                }
  df_H0 = pd.DataFrame(data=dict_H0)

  
  return adjusted_R_2, sd_sample, t_student_coef, critical_value, df_H0
```


```
# Minimum Least Squares

# XtX_inv
XtX_inv = np.linalg.inv(X.T @ X)
# coeficient vector  
b = XtX_inv @ X.T @ y

# R^2 (quality of estimative)
R_2, squared_residuals_sum = R2_SqResidSum(X, y, b)

# t-Student stats (statistical relevance)
(adjusted_R_2, sd_sample, t_student_coef, 
 critical_value, df_H0) = t_student_stats(XtX_inv, b,squared_residuals_sum, R_2,
                            confidence=0.975, n=len(y))
df_H0

```

    H0: null hypothesis: the coefficients bi are not relevant
    
    Reject H0 that b0 is not relevant.
    t_student = 5.641 > critical_value = 2.365
    
    Reject H0 that b0 is not relevant.
    p_value = 7.82e-04 <= alpha = 0.025
    
    Reject H0 that b1 is not relevant.
    t_student = 7.307 > critical_value = 2.365
    
    Reject H0 that b1 is not relevant.
    p_value = 1.62e-04 <= alpha = 0.025
    
    Reject H0 that b2 is not relevant.
    t_student = 3.874 > critical_value = 2.365
    
    Reject H0 that b2 is not relevant.
    p_value = 6.10e-03 <= alpha = 0.025
    





<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>t_student</th>
      <th>p_value</th>
      <th>accept_H0_critical_value</th>
      <th>accept_H0_p_value</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>5.641207</td>
      <td>0.000782</td>
      <td>0.0</td>
      <td>0.0</td>
    </tr>
    <tr>
      <th>1</th>
      <td>7.306917</td>
      <td>0.000162</td>
      <td>0.0</td>
      <td>0.0</td>
    </tr>
    <tr>
      <th>2</th>
      <td>-3.873954</td>
      <td>0.006100</td>
      <td>0.0</td>
      <td>0.0</td>
    </tr>
  </tbody>
</table>
</div>



## Validation with Sklearn


```
from sklearn.linear_model import LinearRegression

reg = LinearRegression().fit(X, y)

# coeficients b
b_skt = reg.coef_.copy()
# reference issue: without .copy()
# reg.score changes from 0.88728964083039 to -36.13590607656064
# because it changes reg.coef_[0]
b_skt[0] = reg.intercept_

print('b from sklearn and homemade implementation are equal?', 
      np.allclose(b, b_skt))

# R^2
R_2_skt = reg.score(X, y)
print('R^2 from sklearn and homemade implementation are equal?', 
      np.allclose(R_2, R_2_skt))

# predict
x_new = np.array([[1, 3, 5]])
skt_predict = reg.predict(x_new)
print('predict from sklearn and homemade implementation are equal?', 
      np.allclose(LR_predict(x_new, b), skt_predict))

```

    b from sklearn and homemade implementation are equal? True
    R^2 from sklearn and homemade implementation are equal? True
    predict from sklearn and homemade implementation are equal? True



```

```

# Logistic Regression

a6010ddbef141950eb94ab32e4ddc3beefcdea30
ipynb
Jupyter Notebook
GoogleColab/LinearRegression.ipynb
viniciusriosfuck/dscodenation
b59fb4417f6b348538f500123fe0ed4d048e2505
MIT
2020-11-09T18:18:41.000Z
2020-11-09T18:18:41.000Z
GoogleColab/LinearRegression.ipynb
inaborges/dscodenation
710113aeed64f3302fe207201a967cc7cdc9e8e7
MIT
2020-07-21T17:28:07.000Z
2020-08-05T22:30:35.000Z
GoogleColab/LinearRegression.ipynb
inaborges/dscodenation
710113aeed64f3302fe207201a967cc7cdc9e8e7
MIT
2020-07-21T22:40:33.000Z
2020-07-21T22:40:33.000Z
Qwen/Qwen-72B
1. YES
2. YES