```python
%pylab inline
```

    Populating the interactive namespace from numpy and matplotlib



```python
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function
```


```python
import pandas as pd
```

# Motivation

I have downloaded prices of flats in Wrocław (the web-scraping code is in the notebook `03-scrape-house-prices`).

The scatter plot of the `price` vs the `area` of the flat reveals that there is a strong dependency between price and area. Linear regression finds a *linear* relation between them.


```python
house_prices = pd.read_csv('03-house-prices.csv',index_col=0)

house_prices.head()
```




<div>
<style>
    .dataframe thead tr:only-child th {
        text-align: right;
    }

    .dataframe thead th {
        text-align: left;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>rooms</th>
      <th>area</th>
      <th>price</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>3</td>
      <td>66.80</td>
      <td>367.801</td>
    </tr>
    <tr>
      <th>1</th>
      <td>3</td>
      <td>55.03</td>
      <td>289.458</td>
    </tr>
    <tr>
      <th>2</th>
      <td>4</td>
      <td>70.80</td>
      <td>389.000</td>
    </tr>
    <tr>
      <th>3</th>
      <td>1</td>
      <td>29.02</td>
      <td>195.695</td>
    </tr>
    <tr>
      <th>4</th>
      <td>2</td>
      <td>45.76</td>
      <td>265.754</td>
    </tr>
  </tbody>
</table>
</div>



House area and price seem to be correlated:


```python
house_prices = house_prices[:100]

_ = scatter(house_prices.area, house_prices.price, alpha=0.3)
_ = xlabel('area m2')
_ = ylabel('price kPLN')
```

## Linear regression analytical solution

We will now find the linear relation between house area and price.

Our data forms a table:


```python
house_prices.head()
```




<div>
<style>
    .dataframe thead tr:only-child th {
        text-align: right;
    }

    .dataframe thead th {
        text-align: left;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>rooms</th>
      <th>area</th>
      <th>price</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>3</td>
      <td>66.80</td>
      <td>367.801</td>
    </tr>
    <tr>
      <th>1</th>
      <td>3</td>
      <td>55.03</td>
      <td>289.458</td>
    </tr>
    <tr>
      <th>2</th>
      <td>4</td>
      <td>70.80</td>
      <td>389.000</td>
    </tr>
    <tr>
      <th>3</th>
      <td>1</td>
      <td>29.02</td>
      <td>195.695</td>
    </tr>
    <tr>
      <th>4</th>
      <td>2</td>
      <td>45.76</td>
      <td>265.754</td>
    </tr>
  </tbody>
</table>
</div>



## Overview of learning

We have gathered a dataset consisting of pairs $\{(x^{(i)}, y^{(i)}), i=1\ldots N \}$. We want to find the relation between $x$, the *independent* variable, area in our example and $y$, the *dependent* variable, price in the example.

We will define a family of *hypotheses* (functions) parameterized by a vector of parameters $\Theta$. 

$$ y \approx h(x; \Theta). $$

To *learn* the relation between areas and prices, or in general between $x$ and $y$, we will search the space of hypotheses (by changing the parameters $\Theta$) to find the one whose predictions best match the true targets.

We will need a *cost function* to measure the goodness of fit between hypotheses and the data. 

Finally, we will need a training (learning) algorithm that will perform the search for us.

## The hypothesis

We will start with a simple linear family of hypotheses:
$$ h(x;\Theta) = \Theta_0 + \Theta_1 x. $$


## The cost function

We will measure the discrepancy between predictions of a hypothesis and the true $y$ on the dataset by computing the sum of squares of errors:

$$ J(\Theta) = \sum_{i=1}^{N} (y^{(i)}-h(x^{(i)};\Theta))^2 $$

## The training algorithm

We will first find the optimal $\Theta$ by deriving the analytical solution. However we will first generalize the problem a little, by assuming that $x$ can be a vector of values, $x=[x_1,\ldots,x_m]$. We will also assume, that $x_0=1$. Therefore:
    
$$ h(x;\Theta) = \Theta{}_0 + \sum_{j=1}^m \Theta_j x_j = \sum_{j=0}^m \Theta_j x_j = \Theta^T x $$

The cost function then becomes

$$ J(\Theta) = \sum_{i=1}^N (y^{(i)} - \Theta^T x^{(i)})^2 $$

We will make even greater use of matrix calculus. If each data sample $x$ is a vector, we can stack them into a large matrix $X \in \mathbb{R}^{m\times N}$, such that the $i$-th column of $X$ is the sample $x^{(i)}$ likewise we will concatenate all targets $y^{(i)}$ into a matrix $Y\in\mathbb{R}^{1 \times N}$. Armed with this notation, we can compute the predictions on all data points as $\Theta^T X$ and compute the matrix of errors $E\in \mathbb{R}^{1 \times N}$ as:

$$ E = \Theta^T X - Y. $$

Finally, please note that the cost function becomes simply:

$$ J(\Theta) = E E^T = (\Theta^T X - Y)(\Theta^T X - Y)^T $$

We will compute the *gradient* of $J$ with respect to $\Theta$. By definition, we will assume that the gradient is also a vector of the same shape as $\Theta$. To quickly compute the gradient, we will use the matrix calculus relations described here: http://www.ee.ic.ac.uk/hp/staff/dmb/matrix/calculus.html and here: http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3274/pdf/imm3274.pdf:

$$ \frac{\partial J(\Theta)}{\partial\Theta} = \frac{\partial}{\partial\Theta} (\Theta^T X - Y)(\Theta^T X - Y)^T =
\frac{\partial}{\partial\Theta} [\Theta^T X X^T \Theta - YX^T\Theta - \Theta^T X Y^T + YY^T] = 
2 X X^T \Theta -2 X Y^T  = 2X(X^T\Theta - Y^T) = 2X E^T$$

Solving for the zero-crossing of the derivative we obtain:

$$
\begin{align}
2 X X^T \Theta &= 2 X Y^T \\
\Theta &= (X X^T) ^{-1} X Y^T 
\end{align}
$$

Please note: this is also called the Moore-Penrose pseudo-inverse of matrix $X$.

## Example of operation


```python
X = np.vstack((np.ones_like(house_prices.area), house_prices.area))
X[:,:10]
```




    Shape: (2, 10)
    [[  1.     1.     1.     1.     1.     1.     1.     1.     1.     1.  ]
     [ 66.8   55.03  70.8   29.02  45.76  38.99  55.7   76.42  37.6   59.51]]




```python
Y = np.asarray(house_prices.price)[None,:]
Y
```




    Shape: (1, 100)
    [[ 367.801  289.458  389.     195.695  265.754  236.614  317.967  409.
       227.211  351.704  259.893  675.615  413.699  302.305  314.908  389.
       270.     527.     347.     390.     439.946  390.     895.     730.
       313.5    251.     289.9    384.     417.7    421.5    285.     374.
       455.     604.7    339.     359.9    545.     239.729  333.196  320.737
       285.     365.     240.     359.     410.     429.     315.     269.
       305.     270.5    245.     299.     183.     300.     489.9    595.
       207.928  238.     349.9    185.     323.4    279.     320.     245.2
       249.     257.     310.     494.     244.57   208.8    185.13   197.
       380.     340.     361.     425.     399.     499.     465.     188.
       276.85   340.     330.     295.     369.     444.     444.     494.
       499.     504.     385.     325.     424.427  650.76   210.     279.
       306.735  349.     466.5    524.064]]



Remember, $\Theta = (X X^T)^{-1} X Y^T$.


```python
Theta = np.linalg.inv(X.dot(X.T)).dot(X).dot(Y.T)
Theta
```




    Shape: (2, 1)
    [[ 53.89879917]
     [  4.72624117]]




```python
scatter(house_prices.area, house_prices.price, alpha=0.3, label='data')

test_areas = np.linspace(xlim()[0], xlim()[1], 512)
test_prices = Theta[0] + Theta[1]*test_areas

plot(test_areas, test_prices,'g', label='lin. regression')
xlim(test_areas[0], test_areas[-1])
xlabel('area m2')
ylabel('price kPLN')
text(150, 500, "$y=%.2fx+%.2f$" % (Theta[1], Theta[0]))
legend()
```

## Interpretation of the least squares criterion

We will now work out a *probabilistic interpretation* of our model. We will assume that the data is generated according to the following process:
1. $x$ is sampled
2. $y$ is computed as $y=\Theta^T x + \epsilon$, where $\epsilon \sim \mathcal{N}(\mu=0,\sigma)$

We assume that the dependent variable $y$ is produced by linearly transforming $x$ and adding a little bit of Gaussian noise. The intuition is that there is the main linear relation that we model, and all other phenomena add up as noise, which by the central limit theorem we assume to be normal.

We assume that for each sample in our dataset the noise was sampled independently from the same noise distribution (IID). This allows us to write the following conditional distribution of $y$ given $x$:

$$ p(y|x;\Theta) = \mathcal{N}(y; \mu=\Theta^T x, \sigma) $$

Due to the IID assumption for the noise variables $\epsilon^{(i)}$, we can also assume that the targets $y$ in our data aree conditionally independent given the $x$s. This allows us to easily compute the probability of observing the collected data under our data-generating hypothesis as:

$$ p(Data; \Theta) = \prod_{i=1}^N p(y^{(i)}|x^{(i)};\Theta) = \prod_{i=1}^N \mathcal{N}(y^{(i)}; \mu=\Theta^T x^{(i)}, \sigma) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(y^{(i)}-\Theta^Tx^{(i)})^2}{2\sigma^2}\right)$$

We will use the maximum likelihood principle to choose $\Theta$. More specifically, this time we will minimize the *negative log-likelihood*, written in short nll:

$$ 
\begin{align}
nll(\Theta) &= - \log p(Data; \Theta) \\
&= -\sum_{i=1}^{N}\log \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(y^{(i)}-\Theta^Tx^{(i)})^2}{2\sigma^2}\right)  \\
&= -\sum_{i=1}^{N}\log \frac{1}{\sqrt{2\pi}\sigma} - \frac{(y^{(i)}-\Theta^Tx^{(i)})^2}{2\sigma^2} \\
&= C + \frac{1}{2\sigma^2}\sum_{i=1}^{N}(y^{(i)}-\Theta^Tx^{(i)})^2 
\end{align}
$$

Note that minimization of nll leads to maximization of the likelihood. Therefore the Gaussian noise assumption leads to the least squares cost function! 

## Validity of the Gaussian noise assumption

Statisticians like to check if the residuals (the errors made by our model on the train set) follow indeed the assumed noise distribution. In the house prices example we can see that the residuals do not strictly follow the Gaussian distribution. Maybe we should take into account another house feature, such as neighborhood, age, whether it is a new or used one... Think how you would incorporate such new information. Is it still possible to use linear regression to fit a model dependent on many variables?

In the following example the Gaussian noise assumption will be satisfied. We will download data of people's height vs weight.


```python
# url = './03-SOCR_Data_Dinov_020108_HeightsWeights.html'
url = 'http://socr.ucla.edu/docs/resources/SOCR_Data/SOCR_Data_Dinov_020108_HeightsWeights.html'
height_weight = pd.read_html(url)[0]
```


```python
height_weight.head()
```




<div>
<style>
    .dataframe thead tr:only-child th {
        text-align: right;
    }

    .dataframe thead th {
        text-align: left;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>0</th>
      <th>1</th>
      <th>2</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>Index</td>
      <td>Height(Inches)</td>
      <td>Weight(Pounds)</td>
    </tr>
    <tr>
      <th>1</th>
      <td>1</td>
      <td>65.78331</td>
      <td>112.9925</td>
    </tr>
    <tr>
      <th>2</th>
      <td>2</td>
      <td>71.51521</td>
      <td>136.4873</td>
    </tr>
    <tr>
      <th>3</th>
      <td>3</td>
      <td>69.39874</td>
      <td>153.0269</td>
    </tr>
    <tr>
      <th>4</th>
      <td>4</td>
      <td>68.2166</td>
      <td>142.3354</td>
    </tr>
  </tbody>
</table>
</div>




```python
inch_to_meter = 0.0254
height = np.asarray(height_weight.iloc[1:,1].astype(float))[None,:] * inch_to_meter #convert to meters
pound_to_kg = 0.453592
weight = np.asarray(height_weight.iloc[1:,2].astype(float))[None,:] * pound_to_kg # convert to kg
```


```python
scatter(height, weight, alpha=0.2, linewidths=0)
xlabel('height [m]')
ylabel('weight [kg]')
```


```python
X = np.vstack((np.ones_like(height), height))
X
```




    Shape: (2, 25000)
    [[ 1.          1.          1.         ...,  1.          1.          1.        ]
     [ 1.67089607  1.81648633  1.762728   ...,  1.64334317  1.71524117
       1.74949129]]




```python
Y = weight
Y
```




    Shape: (1, 25000)
    [[ 51.25249406  61.90954738  69.41177762 ...,  53.64428468  59.99579737
       56.64193813]]




```python
Theta = np.linalg.inv(X.dot(X.T)).dot(X).dot(Y.T)
Theta
```




    Shape: (2, 1)
    [[-37.45569645]
     [ 55.06457669]]




```python
scatter(height, weight, alpha=0.1)
test_heights = np.linspace(xlim()[0], xlim()[1], 512)
test_weights = Theta[0] + Theta[1]*test_heights

plot(test_heights, test_weights,'g')
xlim(test_heights[0], test_heights[-1])
xlabel('height [m]')
ylabel('weight [kg]')
```


```python
E = Theta.T.dot(X) - Y
E
```




    Shape: (1, 25000)
    [[ 3.29899449  0.65880721 -9.80360316 ..., -0.60998512 -3.00246477
       2.23736295]]




```python
scatter (X[1,:], E, alpha=0.1)
title("Scatter plot of residuals.")
```


```python
counts, bins, patches = hist(E.ravel(),100, normed=True)

mu = E.mean()
std = E.std()
e_range = np.linspace(bins[0], bins[-1], 100)
plot(e_range, 1/np.sqrt(2*np.pi)/std*np.exp(-e_range**2/2/std**2), 'r', linewidth=2)
title("Histogram of residuals\nNote the nice gaussian fit")
```


```python
# Analyze the histogram of errors along a few heights

slice_height = 1.85
selected = np.abs(X[1,:]-slice_height) < 0.05

selected_E = E[0,selected]
counts, bins, patches = hist(selected_E, 10, normed=True)

selected_mu = selected_E.mean()
selected_std = selected_E.std()
selected_e_range = np.linspace(bins[0], bins[-1], 100)
plot(selected_e_range, 1/np.sqrt(2*np.pi)/selected_std*np.exp(-(selected_mu-selected_e_range)**2/2/selected_std**2), 'r', linewidth=2)

title('Histogram of residuals for heights=%f +-1cm' 
      'Slice residual mean: %f stdev: %f'% (slice_height,selected_mu, selected_std))
```


```python
# Are the errors homoscedastic?

err_stats = []

slice_width = 0.02

for slice_height in np.arange(height.min(), height.max(), slice_width):
    selected = np.abs(X[1,:]-slice_height) < slice_width/2.0

    selected_E = E[0,selected]

    selected_mu = selected_E.mean()
    selected_std = selected_E.std()
    err_stats.append((slice_height, selected_mu, selected_std))
    
err_stats = np.asarray(err_stats)

plot(err_stats[:,0], err_stats[:,1], label='mean of residuals in slice')
plot(err_stats[:,0], err_stats[:,2], label='stdev of residuals in slice')
legend(loc='lower right')
title("mean and standard deviation of residuals vs height")
xlabel("height")
```


```python
err_stats
```




    Shape: (19, 3)
    [[  1.53107034e+00   2.14386532e+00   5.23856017e+00]
     [  1.55107034e+00  -3.69963253e+00   2.29752329e+00]
     [  1.57107034e+00   5.59137546e+00   5.63231137e+00]
     [  1.59107034e+00  -1.42154057e-01   4.97655815e+00]
     [  1.61107034e+00  -3.35868391e-01   4.36489480e+00]
     [  1.63107034e+00  -9.43881477e-03   4.65754915e+00]
     [  1.65107034e+00  -2.56189507e-02   4.70251580e+00]
     [  1.67107034e+00  -1.13242658e-01   4.61848993e+00]
     [  1.69107034e+00   7.56820063e-02   4.52490018e+00]
     [  1.71107034e+00  -3.49821388e-03   4.58388038e+00]
     [  1.73107034e+00   1.24413620e-01   4.56959502e+00]
     [  1.75107034e+00  -1.38191096e-02   4.56081512e+00]
     [  1.77107034e+00  -1.83136847e-01   4.56016427e+00]
     [  1.79107034e+00   7.44721543e-02   4.52788329e+00]
     [  1.81107034e+00  -2.09913184e-01   4.50193216e+00]
     [  1.83107034e+00   4.32546349e-01   4.41835785e+00]
     [  1.85107034e+00   7.55783464e-04   4.54527322e+00]
     [  1.87107034e+00   3.93371931e-01   4.22798409e+00]
     [  1.89107034e+00   1.83819411e-02   4.74623216e+00]]




```python

```

2261802d8369d6053e896a9add60a34a7613eaaf
ipynb
Jupyter Notebook
lecture_notebooks/03-least-squares.ipynb
iCarrrot/NN
6c402f086342b52209ae254614b065a4bc6feee2
MIT
lecture_notebooks/03-least-squares.ipynb
iCarrrot/NN
6c402f086342b52209ae254614b065a4bc6feee2
MIT
lecture_notebooks/03-least-squares.ipynb
iCarrrot/NN
6c402f086342b52209ae254614b065a4bc6feee2
MIT
Qwen/Qwen-72B
1. YES
2. YES