STAT 39000: Project 8 — Spring 2022

While data frames are a great way to work with data, they are not the only way. Many high performance parts of code are written using a library like numpy or pytorch. These libraries are optimized to be extremely efficient with computation. Multiplying, transposing, inversing matrices can take time, and these libraries can make you code run blazingly fast.

This project is the first project in a series of projects focused on the pytorch library. It is difficult to understand why a library like pytorch is useful without introducing some math. This series of projects will involve some math, however, only at a very high level. Some intuition will be presented as notes, but what is really needed is the ability to read some formulas, and perform the appropriate computations. Throughout this series of projects, we will do our best to ensure that math or statistics is not at all a barrier to completing these projects and getting familiar with pytorch. If it does end up an issue, please post in Piazza and we will do our best to address any issues as soon as possible.

This first project will start slowly, and only focus on the numpy -like functionality of pytorch. We’ve provided you with a set of 100 observations. 75 of the observations are in the train.csv file, 25 are in the test.csv file. We will build a regression model using the data in the train.csv file. In addition, we will calculate some other statistics. Finally, we will (optionally) test our model out on new data in the test.csv dataset.

Start by reading the train.csv file into a pytorch tensor.

Use matplotlib or plotly to plot the data on a scatterplot — x_train on the x-axis, and y_train on the y-axis. After talking to your colleague, you agreed that the data is clearly following a 2nd order polynomial. Something like:

$y = \beta_0 + \beta_1x + \beta_2x^2$

Our goal will be to estimate the values of $\beta_0$, $\beta_1$, and $\beta_2$ using the data in x_train and y_train. Then, we will have a model that could look something like:

$y = 1.2 + .4x + 2.2x^2$

Then, for any given value of x, we can use our model to predict the value of y.

In order to build our model, we need to estimate our parameters: $\beta_0$, $\beta_1$, and $\beta_2$. Luckily, linear regression has a closed form solution, so we can calculate these values directly with the following equation.

$\hat{\beta} = (X^{T} X)^{-1} X^{T} y$

What do these symbols mean? X is a matrix (or tensor), where each column is a term in the polynomial, and each row is an observation. So, for our polynomial, if our X data was simply: 1, 2, 3, 4, the X matrix (or design matrix) would be the following:

Here, the first column is the constant term, the second column is the term of x, the third column is the term of $x^2$, and so on.

When we raise the matrix to the "T" this means to transpose the matrix. The transpose of X, for example, would look like:

When we raise the matrix to the "-1" this means to invert the matrix.

Finally, placing these matrices next to each other means we need to perform matrix multiplication.

pytorch has built in functions to do all of these operations: torch.mm, mat.T, and torch.inverse.

Lastly, y is the tensor containing the observations in y_train.

Start by creating a new tensor called x_mat that is 75 rows and 3 columns. The first column should be filled with 1’s (using torch.ones(x_train.shape[0]).reshape(75,1)), the second column should be the values in x_train, the third column should be the values in x_train squared. Use torch.cat to combine the 75x1 tensors into a single 75x3 tensor (x_mat).

Calculate our estimates for $\beta_0$, $\beta_1$, and $\beta_2$, and save the values in a tensor called betas. The following should be the successful result.

Now that you know the values for $\beta_0$, $\beta_1$, and $\beta_2$, what is our model (as an equation)? It should be:

$y = 4.3677-1.7885x+.4840x^2$

That is pretty cool, and very fast. Now, for any given value of x, we can predict a value of y. Of course, we could write a predict function that accepts a value x, and returns our prediction y, and apply that function to each of the x values in our x_train tensor, however, this can be accomplished even faster and more flexibly using matrix multiplication — simply use the following formula:

$\hat{y} = X\hat{\beta}$

Where X is the x_mat tensor from earlier, and $\hat{\beta}$ is the betas tensor from question (2). Use torch.mm to multiply the two matrices together. Save the resulting tensor to a variable called y_predictions. Finally, create side by side scatterplots. In the first scatterplot, put the values in x_train on the x-axis and the values of y_train on the y-axis. In the second scatterplot put the values of x_train on the x-axis, and your predictions (y_predictions) on the y-axis.

Very cool! Your model should be killing it (after all, we generated this data to follow a known distribution).

To better understand our model, let us create one of the most common forms of uncertainty quantification, confidence intervals. Confidence intervals (95% confidence intervals) show you the range of values (for each x) where we are 95% confident that the average value y for a given x is within the range.

The formula is the following:

$\hat{y_h} \pm t_{(\alpha/2, n-p)} * \sqrt{MSE * diag(x_h(X^{T} X)^{-1} x_h^{T})}$

$MSE = \frac{1}{n-p}\sum_{i=1}^{n}(Y_i - \hat{Y_i})^2$

Since we are calculating the 95% confidence interval for the values of x in our x_train tensor, we can simplify this to:

$\hat{Y} \pm 1.993464 * \sqrt{MSE * diag(X(X^{T} X)^{-1} X^{T})}$

$\frac{1}{72}\sum_{i=1}^{n}(Y_i - \hat{Y_i})^2$

Where:

Create a scatterplot of x_train on the x-axis, and y_predictions on the y-axis. Add the confidence intervals to the plot.

Great! It is unsurprising that our model is a great fit.