# Linear Regression



## The problem

Linear regression assumes that the predictor function is:

$$\hat{y} = \mathbf{w}^\top \cdot \phi(\mathbf{x}) $$

With the basis function $\phi$ equal to the identity function: 

$$\hat{y} = \mathbf{w}^\top \cdot \mathbf{x} $$

defined by parameters $\mathbf{w}$.

Given the data $D = \{(\mathbf{x_i}, y_i)\}_{i=1}^m$, we can compute the performance of an hypothesis with the Mean Square Error:

$$MSE_D = \frac{1}{m} \sum_i (\hat{y} - y)  $$

## Closed-form solution

$$\min MSE_D \implies \nabla_{\mathbf{w}} MSE_D = 0$$

$$ \implies $$

$$ \mathbf{w} = (\mathbf{X}^T \mathbf{X})^{-1}  \mathbf{X}^T \mathbf{y} $$


## Probabilistic interpretation (conditional log-likelihood maximization)

If we put the problem if linear regression in term of maximum likelihood estimation, that is, assuming the data are distributed as:

- $p(y | x) = \mathcal{N}(y; \hat{y}(\mathbf{w},\mathbf{x}), \sigma^2)$

We have that:

- MSE minimization $\equiv$ log-likelihood maximization

Indeed:
$$
\begin{align}
\mathbf{w}_{ML} &= \arg\max_\mathbf{w} p(Y | \mathbf{X}; \mathbf{w}) \\
&= \arg\max_\mathbf{w} \prod_i p(y_i | \mathbf{x_i}; \mathbf{w}) \\
&= \arg\max_\mathbf{w} \sum_i\log p(y_i | \mathbf{x_i}; \mathbf{w}) \\
&= \arg\max_\mathbf{w} - ... - ... - \sum_i \frac{||\hat{y}_i - y_i||^2}{2\sigma^2}
\end{align}
$$

## Bayesian LR

Adding a prior over the parameters and doing MAP estimation is equivalent to common regularization technique (e.g. Lasso/Ridge regression), when the prior is distributed as Gaussian or Laplacian, respectively.



## Univariate Linear Regression

[Gradient Descent based](#Gradient-Descent-based)  
[Closed form](#Closed-form)  
[Using Sklearn](#Using-sklearn)  
[Using Tensorflow](#Using-sklearn)

### Gradient Descent based
 


```python
import numpy as np
import matplotlib.pyplot as plt


def generate_linear_data(w, N):
    """generate data"""
    x = np.random.uniform(low=0, high=20, size=(N, 1))
    # one column of ones for the bias
    x = np.hstack([np.ones((N, 1)), x])
    # linear data with uniform error
    y = x.dot(w) + (2 * np.random.random(N) -1)
    y = y.reshape((N, 1))
    return x, y


x, y = generate_linear_data(np.array([0, 1]), 20)
print(x.shape, y.shape)
plt.scatter(x[:, 1], y)
plt.show()
```


```python
# our hypothesis
def hypothesis(w):
    """In this case, our hypothesis is a line"""
    return lambda x: x.dot(w)


def cost(w, x, y):
    # Mean squared error
    N = len(x)
    h = hypothesis(w)
    hx = h(x).reshape(y.shape)
    # print(x.shape, y.shape)
    # print("h(x)", h(x).shape)
    # print("np.square(h(x)-y)", (h_x -y).shape)
    # print(np.sum(np.square(h(x)-y)))
    return 1/(2*N) * np.sum(np.square(hx-y))


def cost_derivative(w, x, y):
    N = len(x)
    h = hypothesis(w)
    hx = h(x).reshape(y.shape)
    s = (hx-y)*x
    # print(hx.shape)
    # print(s.shape)
    return 1/N * np.sum(s, axis=0)
```


```python
# find parameters that minimize the cost function: Gradient Descent


def gradient_descent(x, y, alpha=0.01, num_steps=1000):
    """Given data, return two parameters, a and b, for the best fit"""
    w = np.random.random(2)
    
    for i in range(num_steps):
        c = cost(w, x, y)
        new_w = w - alpha * cost_derivative(w, x, y)
        if i % 1000 == 0: 
            print("Iteration:\t", i, "\tcost:\t", c, "\tweights:\t", w)
        w = new_w
    return w
```


```python
w = gradient_descent(x,y, alpha=0.01, num_steps=10000)
print("Final weights:", w)
```

    Iteration:	 0 	cost:	 35.4853701003202 	weights:	 [0.34663164 0.18260219]
    Iteration:	 1000 	cost:	 0.16418959105036107 	weights:	 [-0.01697571  0.98774108]
    Iteration:	 2000 	cost:	 0.16411069414476226 	weights:	 [-0.03883738  0.98943947]
    Iteration:	 3000 	cost:	 0.1641104840020508 	weights:	 [-0.03996565  0.98952712]
    Iteration:	 4000 	cost:	 0.16411048344233348 	weights:	 [-0.04002388  0.98953165]
    Iteration:	 5000 	cost:	 0.1641104834408426 	weights:	 [-0.04002688  0.98953188]
    Iteration:	 6000 	cost:	 0.16411048344083853 	weights:	 [-0.04002704  0.98953189]
    Iteration:	 7000 	cost:	 0.1641104834408387 	weights:	 [-0.04002704  0.98953189]
    Iteration:	 8000 	cost:	 0.16411048344083878 	weights:	 [-0.04002704  0.98953189]
    Iteration:	 9000 	cost:	 0.16411048344083856 	weights:	 [-0.04002704  0.98953189]
    Final weights: [-0.04002704  0.98953189]



```python
x_arr = np.arange(20).reshape((20, 1))
x_arr = np.hstack([np.ones((20, 1)), x_arr])
y_arr = hypothesis(w)(x_arr)

plt.scatter(x[:, 1], y)
plt.plot(x_arr[:, 1], y_arr, c="r")
plt.show()
```

### Closed-form


```python
w_closed_form = np.linalg.inv(x.T.dot(x)).dot(x.T).dot(y).reshape(w.size)
print("Closed-form weights:", w_closed_form)
y_closed_form = hypothesis(w_closed_form)(x_arr)

plt.scatter(x[:, 1], y)
plt.plot(x_arr[:, 1], y_arr, c="r")
plt.plot(x_arr[:, 1], y_closed_form, c="b")
plt.show()
```

### Using Sklearn


```python
from sklearn.linear_model import LinearRegression

model = LinearRegression()
model.fit(x, y)
y_sklearn = model.predict(x_arr)
print(model.coef_)
plt.scatter(x[:, 1], y)
plt.plot(x_arr[:, 1], y_arr, c="r")
plt.plot(x_arr[:, 1], y_closed_form, c="b")
plt.plot(x_arr[:, 1], y_sklearn, c="y")
plt.show()
```

### Using Tensorflow


```python
import tensorflow as tf

N = len(x)
X = tf.placeholder(tf.float32)
Y = tf.placeholder(tf.float32)

W = tf.Variable(tf.random_uniform((2, )), name="weights", expected_shape=(2, ))

Y_Pred = tf.reshape(tf.tensordot(X, W, 1), (N, 1))
error = tf.square(Y_Pred - Y)
f_error = tf.reduce_sum(error)/len(x)  # + tf.norm(W)
optimizer = tf.train.GradientDescentOptimizer(0.001).minimize(f_error)

epochs = 1000
with tf.Session() as sess:
    sess.run(tf.global_variables_initializer())
    print("Initial weights: ", W.eval())
    
    loss_expected = 0
    for i in range(epochs):
        
        sess.run(optimizer, {X: x, Y: y})
        
        loss_per_epoch = sess.run(f_error, {X: x, Y: y})
        
        # Allow the training to quit if we've reached a minimum
        if np.abs(loss_expected - loss_per_epoch) < 0.0000001:
            break
        loss_expected = loss_per_epoch
        
        if i%100==0:
            print(loss_expected, "weights: ", W.eval(session=sess))

    print("final loss: ", loss_expected, "weights: ", W.eval(session=sess))
    y_tensorflow = Y_Pred.eval(feed_dict={X: x_arr}, session=sess)
    plt.figure()
    plt.scatter(x[:, 1], y)
    plt.plot(x_arr[:, 1], y_arr, c="r")
    plt.plot(x_arr[:, 1], y_closed_form, c="b")
    plt.plot(x_arr[:, 1], y_sklearn, c="y")
    plt.plot(x_arr[:, 1], y_tensorflow, c="g")
    plt.show()
    

```

## References

- [Chapter 5 of deeplearningbook.org](http://www.deeplearningbook.org/contents/ml.html)


```python

```

bf373ff632f91e1aa0f4ba34411a495a75c52c63
ipynb
Jupyter Notebook
notebooks/linear-regression.ipynb
MarcoFavorito/ml-from-scratch
9295ac5da7fa7328b5a51288eef67b5ca84ab61b
MIT
2019-12-24T03:20:44.000Z
2021-07-28T21:40:45.000Z
notebooks/linear-regression.ipynb
MarcoFavorito/ml-from-scratch
9295ac5da7fa7328b5a51288eef67b5ca84ab61b
MIT
notebooks/linear-regression.ipynb
MarcoFavorito/ml-from-scratch
9295ac5da7fa7328b5a51288eef67b5ca84ab61b
MIT
2020-08-09T03:15:02.000Z
2020-08-09T03:15:02.000Z
Qwen/Qwen-72B
1. YES
2. YES