Python for Data Science Workshop @VeloCity

1.1 Jupyter Notebook

Jupyter notebook is often used by data scientists who work in Python. It is loosely based on Mathematica and combines code, text and visual output in one page.

Some relevant short cuts:

• SHIFT + ENTER executes 1 block of code called a cell
• Tab-completion is omnipresent after the import of a package has been executed
• SHIFT + TAB gives you extra information on what parameters a function takes
• Repeating SHIFT + TAB multiple times gives you even more information

To get used to these short cuts try them out on the cell below.

1.2 Numpy arrays

We'll be working often with numpy arrays so here's a very short introduction.

Basic arithmetical operations on arrays of the same shape are done elementwise:

1.3 Parts to be implemented

In cells like the following example you are expected to implement some code. The remainder of the tutorial won't work if you skip these.

Sometimes assertions are added as a check.

2. Anomaly Detection

2.1 Load Data

First we will load the data using a pickle format. (We use import cPickle as pickle because cPickle is faster.)

The data we use contains the pageviews of one of our own websites and for convenience there is only 1 data point per hour.

2.2 Plot past data

To plot the past data we will use matplotlib.pyplot. For convenience we import it as plt.

% matplotlib inline makes sure you can see the output in the notebook.

(Use % matplotlib notebook if you want to make it ineractive. Don't forget to click the power button to finish the interaction and to be able to plot a new figure.)

2.3 Find the minimum and maximum

Use np.nanmax() and np.nanmin() to find the minmum and maximum while ignoring the NaNs.

And plot these together with the data using the plt.axhline() function.

2.4 Testing the model on unseen data

Now plot all the data instead of just the past data

You can clearly see now that this model does not detect any anomalies. However, the last day of data clearly looks different compared to the other days.

In what follows we will build a better model for anomaly detection that is able to detect these 'shape shifts' as well.

2.5 Building a model with seasonality

To do this we are going to take a step by step approach. Maybe it won't be clear at first why every step is necessary, but that will become clear throughout the process.

First we are going to reshape the past data to a 2 dimensional array with 24 columns. This will give us 1 row for each day and 1 column for each hour. For this we are going to use the np.reshape() function. The newshape parameter is a tuple which in this case should be (-1, 24). If you use a -1 the reshape function will automatically compute that dimension. Pay attention to the order in which the numbers are repositonned (the default ordering should work fine here).

Now we are going to compute the average over all days. For this we are going to use the np.mean() with the axis variable set to the first dimension (axis=0). Next we are going to plot this.

What you can see in the plot above is the average number of pageviews for eacht hour of the day.

Now let's plot this together with the past data on 1 plot. Use a for loop and the np.concatenate() function to concatenate this average 6 times into the variable model.

In the next step we are going to compute the maximum (= positive) and minimum (= negative) deviations from the average to determine what kind of deviations are normal. (Just substract the average/model from the past and take the min and the max of that)

Now let's plot this.

Now let's test this on all data

Now you can clearly see where the anomaly is detected by this more advanced model. The code below can gives you the exact indices where an anomaly is detected. The functions uses are the following np.where() and np.logical_or().

3. Modeling

It is often desired to understand the relationship between different sources of information. As an example we'll consider the historical request rate of a web server and compare it to its CPU usage. We'll try to predict the CPU usage of the server based on the request rates of the different pages. First some imports:

3.1 Data import and inspection

Pandas is a popular library for data wrangling, we'll use it to load and inspect a csv file that contains the historical web request and cpu usage of a web server:

The head command allows one to quickly see the structure of the loaded data:

We can select the CPU column and plot the data:

Next we plot the request rates, leaving out the CPU column as it has another unit:

Now to continue and start to model the data, we'll work with basic numpy arrays. By doing this we also drop the time-information as shown in the plots above.

We extract the column labels as the request_names for later reference:

We extract the request rates as a 2-dimensional numpy array:

and the cpu usage as a one-dimensional numpy array

3.2 Simple linear regression

First, we're going to work with the total request rate on the server, and compare it to the CPU usage. The numpy function np.sum can be used to calculate the total request rate when selecting the right direction (axis=1) for the summation.

Let's plot the total request rate to check:

We can make use of a PyPlot's scatter plot to understand the relation between the total request rate and the CPU usage:

There clearly is a strong correlation between the request rate and the CPU usage. Because of this correlation we can build a model to predict the CPU usage from the total request rate. If we use a linear model we get a formula like the following:

$$ \text{cpu} = c_0 + c_1 \text{total_request_rate} $$

Since we don't know the exact values for $c_0$ and $c_1$ we will have to compute them. For that we'll make use of the scikit-learn machine learning library for Python and use least-squares linear regression

Now we need to feed the data to the model to fit it. The model.fit(X,y) method in general takes an array X and vector y as arguments:

      X = [[x_11, x_12, x_13, ...],                  y = [y_1,
           [x_21, x_22, x_23, ...],                       y_2,  
           [x_31, x_32, x_33, ...],                       y_3,
           ...]                                           ...]

and tries to find coefficients that allow to predict the y_i's from the x_ij's. In our case the matrix X will consist of only 1 column containing the total request rates. Our total_request_rate variable however, is still only a one-dimensional vector, so we need to np.reshape it into a two-dimensional array:

Then we fit our model using the the total request rate and cpu. The coefficients found are automatically stored in the simple_lin_model object.

We can now inspect the coefficient $c_1$ and constant term (intercept) $c_0$ of the model:

So this means that each additional request adds about 0.11% CPU load to the server and all the other processes running on the server consume on average 0.72% CPU.

Once the model is trained we can use it to predict the outcome for a given input (or array of inputs). Note that the predict function requires a 2-dimensional array similar to the fit function.

What is the expected CPU usage when we have 880 requests per second?

Now we plot the linear model together with our data to verify it captures the relationship correctly (the predict method can accept the entire total_request_rate_M array at once).

Our model can calculate a score indicating how well the linear model captures the data. A score of 1 means the data is perfectly linear, a score of 0 (or lower) means the data is not linear at all (and it does not make sense to try to model it that way). The score method takes the same arguments as the fit method:

3.3 Multiple linear regression

Now let us consider the separate request rates instead and build a linear model for that. The model we try to fit takes the form:

$$\text{cpu} = c_0 + c_1 \text{request_rate}_1 + c_2 \text{request_rate}_2 + \ldots + c_n \text{request_rate}_n$$

where the $\text{request_rate}_i$'s correspond the our different requests:

We start again by creating a LinearRegression model.

Next we fit the model on the data, using multi_lin_model.fit(X,y). In contrast to the case above our request_rates variable already has the correct shape to pass as the X matrix: it has one column per request type.

Now, given the coefficients calculated by the model, which capture the contribution of each request type to the total CPU usage, we can start to answer some interesting questions. For example, which request causes most CPU usage, on a per visit basis?

For this we can generate a table of request names with their coefficients in descending order:

From this table we see that 'resources/js/basket.js' consumes the most per CPU per request. It generates about 0.30% CPU load for each additional request. 'products/science.html' on the other hand is much leaner and only consumes about 0.04% CPU per request.

Now let us investigate the constant term again.

As you can see this term is very similar to the result achieved in single linear regression, but it is not entirely the same. This means that these models are not perfect. However, the seem to be able to give a reliable estimate.

3.4 Multiple linear regression 'advanced'

In the previous section we have modeled how much load each individual request generates. But in some cases you might want to transfer one of the request to another server. Now, suppose we want to minimize average CPU usage on this server by deviating traffic of only one webpage to another server, which page should we choose?

For this we simulate diverting the traffic of one page to another server. This means that for the request that is diverted the rate becomes 0, for the other requests we use the average rate.

We implement this by first calculating the average_request_rates using np.mean. These average_request_rates are then fed to the multi_lin_model.predict() method but with setting each individual request rate to 0 once.

(For linear models you can also compute the result based on the coefficients, but this approach also works for non-linear models.)

As you can see in the table above, it is best to divert the traffic of 'api/product/get.php' (Why is the result different than the table based on the coefficient)?

4. Forecasting

For the forecasting we are going to use page views data, very similar to the data used in the anomaly detection section. It is also page view data and contains 1 sample per hour.

In the graph above you can clearly see that there is a rising trend in the data.

4.1 One-step ahead prediction

This forecasting section will describe the one-step ahead prediction. This means in this case that we will only predict the next data point which is in this case the number of pageviews in the next hour.

Now let's first build a model that tries to predict the next data point from the previous one.

As you can see from the score above, the model is not perfect but it seems to get a relatively high score. Now let's make a prediction into the future and plot this.

To predict the datapoint after that we will use the predicted data to make a new prediction. The code below shows how this works for this data set using the linear model you used earlier. Don't forget to fill out the missing code.

As you can see from the image above the model doesn't quite seem to fit the data well. Let's see how we can improve this.

4.2 Multiple features

If your model is not smart enough there is a simple trick in machine learning to make your model more intelligent (but also more complex). This is by adding more features.

To make our model better we will use more than 1 sample from the past. To make your life easier there is a simple function below that will create a data set for you. The width parameter sets the number of hours in the past that will be used.

As you can see from the print above both X and y contains 303 datapoints. For X you see that there are now 5 features which contain the pageviews from the 5 past hours.

So let's have a look what the increase from 1 to 5 features results to.

Now change the width parameter to see if you can get a better score.

4.3 Over-fitting

Now execute the code below to see the prediction of this model.

As you can see in the image above the prediction is not what you would expect from a perfect model. What happened is that the model learned the training data by heart without 'understanding' what the data is really about. This fenomenon is called over-fitting and will always occur if you make your model more complex.

Now play with the width variable below to see if you can find a more sensible width.

As you will have noticed by now is that it is better to have a non-perfect score which will give you a much better outcome. Now try the same thing for the following models:

• sklearn.linear_model.RidgeCV()
• sklearn.linear_model.LassoCV()
• sklearn.gaussian_process.GaussianProcess()

The first 2 models also estimate the noise that is present in the data to avoid overfitting. RidgeCV will keep the weights that are found small, but it won't put them to zero. LassoCV on the other hand will put several weights to 0. Execute model.coef_ to see the actual coefficients that have been found.

GaussianProcess is a non-linear method. This makes this method a lot more complex and therefore it will need significantly less features to be able to learn the data by hart (and thus to over-fit). In many cases however this additional complexity allows to better understand the data. Additionally it has the advantage that it can estimate confidance intervals similar to the red lines used in the anomaly detection.

4.4 Automation

What we have done up to now is manually selecting the best outcome based on the test result. This can be considered cheating because you have just created a self-fulfilling prophecy. Additionally it is not only cheating it is also hard to find the exact width that gives the best result by just visually inspecting it. So we need a more objective approach to solve this.

To automate this process you can use a validation set. In this case we will use the last 48 hours of the training set to validate the score and select the best parameter value. This means that we will have to use a subset of the training set to fit the model.

If everything is correct the LassoCV methods was selected.

Now we are going to train this best model on all the data. In this way we use all the available data to build a model.

Altough the optimal result found here might not be the best visually, it is a far better result than the one you selected manually just because there was no cheating involved ;-).

Some additional info:

• This noise level of RidgeCV and LassoCV is estimated by automatically performing train and validation within the method itself. This will make them much more robust against over-fitting. The actual method used is Cross-validation) which is a better approach of what we do here because it repeats the training and validation multiple times for different training and validation sets. The parameter that is set for these methods is often called the regularization parameter in literature and is well suited to avoid over-fitting.
• Although sklearn supports estimating the noise level in Gaussian Processes it is not implemented within the method itself. Newer versions of sklearn seem to entail a lot of changes in this method so possibly it will be integrated in the (near) future. If you want to implement this noise level estimation yourself you can use their cross-validation tool to set the alpha parameter in version 0.18 of sklearn. (The version used here is 0.17.)

5. Challenge

Details will be given near the end of the tutorial.