<jupyter_start><jupyter_text>MAC0460 / MAC5832 (2021) EP2: Linear regression - analytic solution Objectives:- to implement and test the analytic solution for the linear regression task (see, for instance, Slides of Lecture 03 and Lecture 03 of *Learning from Data*)- to understand the core idea (*optimization of a loss or cost function*) for parameter adjustment in machine learning  Linear regressionGiven a dataset $\{(\mathbf{x}^{(1)}, y^{(1)}), \dots ,(\mathbf{x}^{(N)}, y^{(N)})\}$ with $\mathbf{x}^{(i)} \in \mathbb{R}^{d}$ and $y^{(i)} \in \mathbb{R}$, we would like to approximate the unknown function $f:\mathbb{R}^{d} \rightarrow \mathbb{R}$ (recall that $y^{(i)} =f(\mathbf{x}^{(i)})$) by means of a linear model $h$:$$h(\mathbf{x}^{(i)}; \mathbf{w}, b) = \mathbf{w}^\top  \mathbf{x}^{(i)} + b$$Note that $h(\mathbf{x}^{(i)}; \mathbf{w}, b)$ is, in fact, an  [affine transformation](https://en.wikipedia.org/wiki/Affine_transformation) of  $\mathbf{x}^{(i)}$. As commonly done, we will use the term "linear" to refer to an affine transformation.The output of $h$ is a linear transformation of $\mathbf{x}^{(i)}$. We use the notation $h(\mathbf{x}^{(i)}; \mathbf{w}, b)$ to make clear that $h$ is a parametric model, i.e., the transformation $h$ is defined  by the parameters $\mathbf{w}$ and $b$. We can view vector $\mathbf{w}$ as a *weight* vector that controls the effect of each *feature* in the prediction.By adding one component with value equal to 1 to the observations $\mathbf{x}$ (an artificial coordinate), we have:$$\tilde{\mathbf{x}} = (1, x_1, \ldots, x_d) \in \mathbb{R}^{1+d}$$and then we can simplify the notation:$$h(\mathbf{x}^{(i)}; \mathbf{w}) = \hat{y}^{(i)} = \mathbf{w}^\top  \tilde{\mathbf{x}}^{(i)}$$We would like to determine the optimal parameters $\mathbf{w}$ such that prediction $\hat{y}^{(i)}$ is as closest as possible to $y^{(i)}$ according to some error metric. Adopting the *mean square error* as such metric we have the following cost function:\begin{equation}J(\mathbf{w}) = \frac{1}{N}\sum_{i=1}^{N}\big(\hat{y}^{(i)} - y^{(i)}\big)^{2}\end{equation}Thus, the task of determining a function $h$ that is closest to $f$ is reduced to the task of finding the values $\mathbf{w}$ that minimize $J(\mathbf{w})$.**Now we will explore this model, starting with a simple dataset.**  Auxiliary functions<jupyter_code># some imports
import numpy as np
import time
import matplotlib.pyplot as plt

%matplotlib inline
# An auxiliary function
def get_housing_prices_data(N, verbose=True):
    """
    Generates artificial linear data,
    where x = square meter, y = house price

    :param N: data set size
    :type N: int
    
    :param verbose: param to control print
    :type verbose: bool
    :return: design matrix, regression targets
    :rtype: np.array, np.array
    """
    cond = False
    while not cond:
        x = np.linspace(90, 1200, N)
        gamma = np.random.normal(30, 10, x.size)
        y = 50 * x + gamma * 400
        x = x.astype("float32")
        x = x.reshape((x.shape[0], 1))
        y = y.astype("float32")
        y = y.reshape((y.shape[0], 1))
        cond = min(y) > 0
        
    xmean, xsdt, xmax, xmin = np.mean(x), np.std(x), np.max(x), np.min(x)
    ymean, ysdt, ymax, ymin = np.mean(y), np.std(y), np.max(y), np.min(y)
    if verbose:
        print("\nX shape = {}".format(x.shape))
        print("y shape = {}\n".format(y.shape))
        print("X: mean {}, sdt {:.2f}, max {:.2f}, min {:.2f}".format(xmean,
                                                               xsdt,
                                                               xmax,
                                                               xmin))
        print("y: mean {:.2f}, sdt {:.2f}, max {:.2f}, min {:.2f}".format(ymean,
                                                                 ysdt,
                                                                 ymax,
                                                                 ymin))
    return x, y
# Another auxiliary function
def plot_points_regression(x,
                           y,
                           title,
                           xlabel,
                           ylabel,
                           prediction=None,
                           legend=False,
                           r_squared=None,
                           position=(90, 100)):
    """
    Plots the data points and the prediction,
    if there is one.

    :param x: design matrix
    :type x: np.array
    :param y: regression targets
    :type y: np.array
    :param title: plot's title
    :type title: str
    :param xlabel: x axis label
    :type xlabel: str
    :param ylabel: y axis label
    :type ylabel: str
    :param prediction: model's prediction
    :type prediction: np.array
    :param legend: param to control print legends
    :type legend: bool
    :param r_squared: r^2 value
    :type r_squared: float
    :param position: text position
    :type position: tuple
    """
    fig, ax = plt.subplots(1, 1, figsize=(8, 8))
    line1, = ax.plot(x, y, 'bo', label='Real data')
    if prediction is not None:
        line2, = ax.plot(x, prediction, 'r', label='Predicted data')
        if legend:
            plt.legend(handles=[line1, line2], loc=2)
        ax.set_title(title,
                 fontsize=20,
                 fontweight='bold')
    if r_squared is not None:
        bbox_props = dict(boxstyle="square,pad=0.3",
                          fc="white", ec="black", lw=0.2)
        t = ax.text(position[0], position[1], "$R^2 ={:.4f}$".format(r_squared),
                    size=15, bbox=bbox_props)

    ax.set_xlabel(xlabel, fontsize=20)
    ax.set_ylabel(ylabel, fontsize=20)
    plt.show()<jupyter_output><empty_output><jupyter_text>The dataset The first dataset we will use is a toy dataset. We will generate $N=100$ observations with only one *feature* and a real value associated to each of them. We can view these observations as being pairs *(area of a real state in square meters, price of the real state)*. Our task is to construct a model that is able to predict the price of a real state, given its area.<jupyter_code>X, y = get_housing_prices_data(N=100)<jupyter_output>X shape = (100, 1)
y shape = (100, 1)

X: mean 645.0, sdt 323.65, max 1200.00, min 90.00
y: mean 44221.50, sdt 16674.25, max 79182.43, min 12196.20<jupyter_text>Ploting the data<jupyter_code>plot_points_regression(X,
                       y,
                       title='Real estate prices prediction',
                       xlabel="m\u00b2",
                       ylabel='$')<jupyter_output><empty_output><jupyter_text>The solutionGiven $f:\mathbb{R}^{N\times M} \rightarrow \mathbb{R}$ and $\mathbf{A} \in \mathbb{R}^{N\times M}$, we define the gradient of $f$ with respect to $\mathbf{A}$ as:\begin{equation*}\nabla_{\mathbf{A}}f = \frac{\partial f}{\partial \mathbf{A}} =  \begin{bmatrix}\frac{\partial f}{\partial \mathbf{A}_{1,1}} & \dots & \frac{\partial f}{\partial \mathbf{A}_{1,m}} \\\vdots &  \ddots & \vdots \\\frac{\partial f}{\partial \mathbf{A}_{n,1}} &  \dots & \frac{\partial f}{\partial \mathbf{A}_{n,m}}\end{bmatrix}\end{equation*}Let $\mathbf{X} \in \mathbb{R}^{N\times d}$ be a matrix (sometimes also called the *design matrix*) whose rows are the observations of the dataset and let $\mathbf{y} \in \mathbb{R}^{N}$ be the vector consisting of all values $y^{(i)}$ (i.e., $\mathbf{X}^{(i,:)} = \mathbf{x}^{(i)}$ and $\mathbf{y}^{(i)} = y^{(i)}$). It can be verified that: \begin{equation}J(\mathbf{w}) = \frac{1}{N}(\mathbf{X}\mathbf{w} - \mathbf{y})^{T}(\mathbf{X}\mathbf{w} - \mathbf{y})\end{equation}Using basic matrix derivative concepts we can compute the gradient of $J(\mathbf{w})$ with respect to $\mathbf{w}$:\begin{equation}\nabla_{\mathbf{w}}J(\mathbf{w}) = \frac{2}{N} (\mathbf{X}^{T}\mathbf{X}\mathbf{w} -\mathbf{X}^{T}\mathbf{y})   \end{equation}Thus, when $\nabla_{\mathbf{w}}J(\mathbf{w}) = 0$ we have \begin{equation}\mathbf{X}^{T}\mathbf{X}\mathbf{w} = \mathbf{X}^{T}\mathbf{y}\end{equation}Hence,\begin{equation}\mathbf{w} = (\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{y}\end{equation}Note that this solution has a high computational cost. As the number of variables (*features*) increases, the cost for matrix inversion becomes prohibitive. See  [this text](https://sgfin.github.io/files/notes/CS229_Lecture_Notes.pdf) for more details.  Exercise 1Using only **NumPy** (a quick introduction to this library can be found  [here](http://cs231n.github.io/python-numpy-tutorial/)), complete the two functions below. Recall that $\mathbf{X} \in \mathbb{R}^{N\times d}$; thus you will need to add a component of value 1 to each of  the observations in $\mathbf{X}$ before performing the computation described above.NOTE: Although the dataset above has data of dimension $d=1$, your code must be generic (it should work for $d\geq1$)  1.1. Weight computation function<jupyter_code>def normal_equation_weights(X, y):
    """
    Calculates the weights of a linear function using the normal equation method.
    You should add into X a new column with 1s.

    :param X: design matrix
    :type X: np.ndarray(shape=(N, d))
    :param y: regression targets
    :type y: np.ndarray(shape=(N, 1))
    :return: weight vector
    :rtype: np.ndarray(shape=(d+1, 1))
    """
    
    # START OF YOUR CODE:
    N = X.shape[0]
    X_tilde = np.column_stack((np.ones((N, 1)), X))
    X_cross = np.dot(np.linalg.inv(np.dot(X_tilde.T, X_tilde)), X_tilde.T)
    w = np.dot(X_cross, y)
    
    return w
    # END OF YOUR CODE
# test of function normal_equation_weights()

w = 0  # this is not necessary
w = normal_equation_weights(X, y)
print("Estimated w =\n", w)<jupyter_output>Estimated w =
 [[12043.47818971]
 [   49.88841379]]<jupyter_text>1.2. Prediction function<jupyter_code>def normal_equation_prediction(X, w):
    """
    Calculates the prediction over a set of observations X using the linear function
    characterized by the weight vector w.
    You should add into X a new column with 1s.

    :param X: design matrix
    :type X: np.ndarray(shape=(N, d))
    :param w: weight vector
    :type w: np.ndarray(shape=(d+1, 1))
    :param y: regression prediction
    :type y: np.ndarray(shape=(N, 1))
    """
    
    # START OF YOUR CODE:
    N = X.shape[0]
    X_tilde = np.column_stack((np.ones((N, 1)), X))
    y = np.dot(X_tilde, w)
    
    return y
    # END OF YOUR CODE<jupyter_output><empty_output><jupyter_text>1.3. Coefficient of determinationWe can use the [$R^2$](https://pt.wikipedia.org/wiki/R%C2%B2) metric (Coefficient of determination) to evaluate how well the linear model fits the data.**Which $𝑅^2$ value would you expect to observe ?**<jupyter_code>from sklearn.metrics import r2_score

# test of function normal_equation_prediction()
prediction = normal_equation_prediction(X, w)

# compute the R2 score using the r2_score function from sklearn
# Replace 0 with an appropriate call of the function

# START OF YOUR CODE:
r_2 = r2_score(y, prediction)
# END OF YOUR CODE

plot_points_regression(X,
                       y,
                       title='Real estate prices prediction',
                       xlabel="m\u00b2",
                       ylabel='$',
                       prediction=prediction,
                       legend=True,
                       r_squared=r_2)<jupyter_output><empty_output><jupyter_text>Additional testsLet us compute a prediction for $x=650$<jupyter_code># Let us use the prediction function
x = np.asarray([650]).reshape(1,1)
prediction = normal_equation_prediction(x, w)
print("Area = %.2f  Predicted price = %.4f" %(x[0], prediction))<jupyter_output>Area = 650.00  Predicted price = 44470.9472<jupyter_text>1.4. Processing timeExperiment with different nummber of samples $N$ and observe how processing time varies.Be careful not to use a too large value; it may make jupyter freeze ...<jupyter_code># Add other values for N
# START OF YOUR CODE:
N = [100, 200, 500, 1000, 5000, 10000, 100000, 10000000]
# END OF YOUR CODE

for i in N:
    X, y = get_housing_prices_data(N=i)
    init = time.time()
    w = normal_equation_weights(X, y)
    prediction = normal_equation_prediction(X,w)
    init = time.time() - init
    
    print("\nExecution time = {:.8f}(s)\n".format(init))<jupyter_output>X shape = (100, 1)
y shape = (100, 1)

X: mean 645.0, sdt 323.65, max 1200.00, min 90.00
y: mean 44131.31, sdt 16462.71, max 77546.66, min 10220.15

Execution time = 0.00058222(s)


X shape = (200, 1)
y shape = (200, 1)

X: mean 645.0, sdt 322.04, max 1200.00, min 90.00
y: mean 44643.07, sdt 16741.32, max 77416.36, min 11328.75

Execution time = 0.00007701(s)


X shape = (500, 1)
y shape = (500, 1)

X: mean 645.0, sdt 321.07, max 1200.00, min 90.00
y: mean 44582.73, sdt 16791.18, max 83074.05, min 9172.27

Execution time = 0.00033569(s)


X shape = (1000, 1)
y shape = (1000, 1)

X: mean 645.0, sdt 320.75, max 1200.00, min 90.00
y: mean 44441.20, sdt 16539.07, max 79329.80, min 10420.06

Execution time = 0.00034523(s)


X shape = (5000, 1)
y shape = (5000, 1)

X: mean 645.0, sdt 320.49, max 1200.00, min 90.00
y: mean 44266.45, sdt 16549.26, max 80368.79, min 5516.63

Execution time = 0.00048804(s)


X shape = (10000, 1)
y shape = (10000, 1)

X: mean 645.0000610351562, sdt 320.46, max 1[...]<jupyter_text>Exercise 2Let us test the code with $𝑑>1$. We will use the data we have collected in our first class. The [file](https://edisciplinas.usp.br/pluginfile.php/5982803/course/section/6115454/QT1data.csv) can be found on e-disciplinas. Let us try to predict the weight based on one or more features.<jupyter_code>import pandas as pd

# load the dataset
df = pd.read_csv('QT1data.csv')
df.head()
df.describe()
# Our target variable is the weight
y = df.pop('Weight').values
y<jupyter_output><empty_output><jupyter_text>2.1. One feature ($d=1$)We will use 'Height' as the input feature and predict the weight<jupyter_code>feature_cols = ['Height']
X = df.loc[:, feature_cols]
X.shape<jupyter_output><empty_output><jupyter_text>Write the code for computing the following- compute the regression weights using $\mathbf{X}$ and $\mathbf{y}$- compute the prediction- compute the $R^2$ value- plot the regression graph (use appropriate values for the parameters of function plot_points_regression())<jupyter_code># START OF YOUR CODE:
w = normal_equation_weights(X, y)
prediction = normal_equation_prediction(X, w)
r_2 = r2_score(y, prediction)
plot_points_regression(X,
                       y,
                       title='Weight based on height prediction',
                       xlabel="Height",
                       ylabel='Weight',
                       prediction=prediction,
                       legend=True,
                       r_squared=r_2)
# END OF YOUR CODE<jupyter_output><empty_output><jupyter_text>2.2 - Two input features ($d=2$)Now repeat the exercise with using as input the features 'Height' and 'Shoe number'- compute the regression weights using $\mathbf{X}$ and $\mathbf{y}$- compute the prediction- compute and print the $R^2$ valueNote that our plotting function can not be used. There is no need to do plotting here.<jupyter_code># START OF YOUR CODE:
feature_cols = ['Height', 'Shoe number']
X = df.loc[:, feature_cols]

w = normal_equation_weights(X, y)
prediction = normal_equation_prediction(X, w)
r_2 = r2_score(y, prediction)

print(r_2)
# END OF YOUR CODE<jupyter_output>0.45381183096658584<jupyter_text>2.3 - Three input features ($d=3$)Now try with three features. There is no need to do plotting here.- compute the regression weights using $\mathbf{X}$ and $\mathbf{y}$- compute the prediction- compute and print the $R^2$ value<jupyter_code># START OF YOUR CODE:
feature_cols = ['Height', 'Shoe number', 'Age']
X = df.loc[:, feature_cols]

w = normal_equation_weights(X, y)
prediction = normal_equation_prediction(X, w)
r_2 = r2_score(y, prediction)

print(r_2)
# END OF YOUR CODE<jupyter_output>0.4776499498669615
{'fraction_non_alphanumeric': 0.13863863863863865, 'fraction_numerical': 0.04504504504504504, 'mean_word_length': 3.061229674796748, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 49, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '22673', 'license': 'MIT', 'n_tokens_mistral': 5718, 'n_tokens_neox': 5126, 'n_words': 1674, 'path': 'ep02_linreg_analytic.ipynb', 'repo_name': 'ricardokojo/MAC0460-2021', 'size': '15984'}