---
title: "Data Science and Predictive Analytics (UMich HS650)"
subtitle: "<h2><u>Black Box Machine-Learning Methods: Neural Networks and Support Vector Machines</u></h2>"
author: "<h3>SOCR/MIDAS (Ivo Dinov)</h3>"
date: "`r format(Sys.time(), '%B %Y')`"
tags: [DSPA, SOCR, MIDAS, Big Data, Predictive Analytics] 
output:
  html_document:
    theme: spacelab
    highlight: tango
    includes:
      before_body: SOCR_header.html
      after_body: SOCR_footer_tracker.html
    toc: true
    number_sections: true
    toc_depth: 2
    toc_float:
      collapsed: false
      smooth_scroll: true
---

In this chapter, we are going to cover two very powerful machine-learning algorithms. These techniques have complex mathematical formulations, however, efficient algorithms and reliable software packages have been developed to utilize them for various practical applications. We will (1) describe Neural Networks as analogues of biological neurons, (2) develop hands-on a neural net that can be trained to compute the square-root function, (3) describe support vector machine (SVM) classification, and (4) complete several case-studies, including optical character recognition (OCR), the Iris flowers, Google Trends and the Stock Market, and Quality of Life in chronic disease.

Later, in [Chapter 22](http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/notes/22_DeepLearning.html), we will provide more details and additional examples of deep neural network learning. For now, let's start by exploring the *magic* inside the machine learning black box.

# Understanding Neural Networks

## From biological to artificial neurons

An `Artificial Neural Network (ANN)` model mimics the biological brain response to multisource (sensory-motor) stimuli (inputs). ANN simulate the brain using a network of interconnected neuron cells to create a massive parallel processor. Of course, it uses a network of artificial nodes, not brain cells, to train data.

When we have three signals (or inputs) $x_1$, $x_2$ and $x_3$, the first step is weighting the features ($w$'s) according to their importance. Then, the weighted signals are summed by the "neuron cell" and this sum is passed on according to an **activation function** denoted by **f**. The last step is generating an output **y** at the end of the process. A typical output will have the following mathematical relationship to the inputs.
$$y(x)=f\left (\sum_{i=1}^n w_ix_i\right ).$$

There are three important components for building a neural network:

* **Activation function**: transforms weighted and summed inputs to output.
* **Network topology**: describes the number of "neuron cells", the number of layers and manner in which the cells are connected.
* **Training algorithm**: how to determine weights $w_i$

Let's look at each of these components one by one.

## Activation functions

One of the functions is known as threshold activation function that results in an output signal once a specified input threshold has been attained.

```{r echo=FALSE}
x<-runif(1000, -10, 10)
y<-ifelse(x>=0, 1, 0)
plot(x, y, xlab = "Sum of input signals", ylab = "Output signal", main = "Threshold")
abline(v=0, lty=2)
```

$$
f(x)= \left\{
\begin{array}{ll}
      0 & x<0 \\
      1 & x\geq 0 \\
\end{array} 
\right. .
$$

This is the simplest form for activation functions. It seldom used in real world situations. Most commonly used alternative is the sigmoid activation function where $f(x)=\frac{1}{1+e^{-x}}$. The *Euler number* `e` is the base of natural logarithms. The output signal is no longer binary but can be any real number ranged from 0 to 1. 

```{r echo=F}
x<-runif(100000, -10, 10)
y<-1/(1+exp(-x))
plot(x, y, xlab = "Sum of input signals", ylab = "Output signal", main = "Sigmoid")
```

Other activation functions might also be useful:

```{r echo=FALSE, fig.height=6}
x<-runif(100000, -10, 10)
y1<-x
y2<-ifelse(x<=-5, -5, ifelse(x>=5, 5, x))
y3<-(exp(x)-exp(-x))/(exp(x)+exp(-x))
y4<-exp(-x^2/2)
par(mfrow=c(2, 2))
plot(x, y1, main="Linear", xlab="", ylab="")
plot(x, y2, main="Saturated Linear", xlab="", ylab="")
plot(x, y3, main="Hyperbolic tangent", xlab="", ylab="")
plot(x, y4, main = "Gaussian", xlab="", ylab="")
```

Basically, we can chose a proper activation function based on the corresponding codomain of the function. For example, with hyperbolic tangent activation function, we can only have outputs ranging from -1 to 1 regardless of what input do we have. With linear function we can go from $-\infty$ to $+\infty$. Our Gaussian activation function will give us a model called *Radial Basis Function* network.

## Network topology

The number of layers: The $x$'s or features in the dataset is called **input nodes** while the predicted values are called the **output nodes**. Multilayer networks include multiple hidden layers. The following graph represents a two layer neural network: 

![](http://wiki.socr.umich.edu/images/3/36/DSPA_Figs_TwoLayer_Chapter6.png)
$$\text{Two Layer network.}$$

When we have multiple layers, the information flow could be complicated.

## The direction of information travel

The arrows in the last graph (with multiple layers) suggest a feed forward network. In such network, we can also have multiple outcomes modeled simultaneously.

![](http://wiki.socr.umich.edu/images/0/07/DSPA_Figs_multioutput_Chapter6.png)
$$\text{Multi-output}$$

Alternatively, in a recurrent network (feedback network), information can also travel backwards in loops (or delay). This is illustrated in the following graph.  

![](http://wiki.socr.umich.edu/images/f/f9/DSPA_Figs_delay_chapter6.png)
$$\text{Delay(feedback network)}$$

This short-term memory increases the power of recurrent networks dramatically. However, in practice, recurrent networks are seldom used.

## The number of nodes in each layer

Number of input nodes and output nodes are predetermined by the dataset and predictive variables. The number we can edit is the hidden nodes in the model. Our goal is to add fewer hidden nodes as possible to simplify the model when the model performance is pleasant. 

## Training neural networks with backpropagation

This algorithm could determine the weights in the model using a strategy of back-propagating errors. First, we assign random numbers for weights (but all weights must be non-trivial, i.e., $\not= 0$). For example, we can use normal distribution, or any other random process, to assign initial weights. Then we adjust the weights iteratively by repeating the process until until certain convergence or stopping criterion is met. Each iteration contains two phases. 

* Forward phase: from input layer to output layer using current weights. Outputs are produced at the end of this phase, and
* Backward phase: compare the outputs and true target values. If the difference is significant, we change the weights and go through the forward phase, again.

In the end, we pick a set of weights, corresponding to the least total error, to be the weights in our network.

# Case Study 1: Google Trends and the Stock Market - **Regression**

## Step 1 - collecting data

In this case study, we are going to use the [Google trends and stock market dataset](http://wiki.socr.umich.edu/index.php/SOCR_Data_GoogleTrends_2005_2011). A [doc file with the meta-data](https://umich.instructure.com/files/416273/download?download_frd=1) and the [CSV data](https://umich.instructure.com/files/416274/download?download_frd=1) are available on the [Case-Studies Canvas Site](https://umich.instructure.com/courses/38100/files/folder/Case_Studies). These daily data (between 2008 and 2009) can be used to examine the associations between Google search trends and the daily marker index - Dow Jones Industrial average.

### Variables

* **Index**: Time Index of the Observation
* **Date**: Date of the observation (Format: YYYY-MM-DD)
* **Unemployment**: The Google Unemployment Index tracks queries related to "unemployment, social, social security, unemployment benefits" and so on.
* **Rental**: The Google Rental Index tracks queries related to "rent, apartments, for rent, rentals", etc.
RealEstate: The Google Real Estate Index tracks queries related to "real estate, mortgage, rent, apartments" and so on.
* **Mortgage**: The Google Mortgage Index tracks queries related to "mortgage, calculator, mortgage calculator, mortgage rates".
* **Jobs**: The Google Jobs Index tracks queries related to "jobs, city, job, resume, career, monster" and so forth.
* **Investing**: The Google Investing Index tracks queries related to "stock, finance, capital, yahoo finance, stocks", etc.
* **DJI_Index**: The Dow Jones Industrial (DJI) index. These data are interpolated from 5 records per week (Dow Jones stocks are traded on week-days only) to 7 days per week to match the constant 7-day records of the Google-Trends data.
* **StdDJI**: The standardized-DJI Index computed by: StdDJI = 3+(DJI-11091)/1501, where m=11091 and s=1501 are the approximate mean and standard-deviation of the DJI for the period (2005-2011).
* **30-Day Moving Average Data Columns**: The 8 variables below are the 30-day moving averages of the 8 corresponding (raw) variables above: *Unemployment30MA, Rental30MA, RealEstate30MA, Mortgage30MA, Jobs30MA, Investing30MA, DJI_Index30MA*, and *StdDJI_30MA*.
* **180-Day Moving Average Data Columns**: The 8 variables below are the 180-day moving averages of the 8 corresponding (raw) variables: *Unemployment180MA, Rental180MA, RealEstate180MA, Mortgage180MA, Jobs180MA, Investing180MA, DJI_Index180MA*, and *StdDJI_180MA*.

Here we use the `RealEstate` as our dependent variable. Let's see if Google Real Estate Index could be predicted by other variables in the dataset.

## Step 2 - exploring and preparing the data

First thing first, we need to load the dataset into R.

```{r echo=T}
google<-read.csv("https://umich.instructure.com/files/416274/download?download_frd=1", stringsAsFactors = F)
```

Let's delete the first two columns, since the only goal is to predict Google Real Estate Index with other indexes and DJI.

```{r echo=T}
google<-google[, -c(1, 2)]
str(google)
```

As we can see from the structure of the data, these indexes and DJI have different ranges. We should `rescale` the data. In [Chapter 6](http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/notes/06_LazyLearning_kNN.html), we learned that normalizing these features using our own `normalize()` function could fix the problem. We can use `lapply()` to apply the `normalize()` function to each column.

```{r echo=T}
normalize <- function(x) {
return((x - min(x)) / (max(x) - min(x)))
}
google_norm<-as.data.frame(lapply(google, normalize))
summary(google_norm$RealEstate)
```

The last step clearly normalizes all feature vectors into the 0 to 1 range.

The next step would be separating our `google` dataset into training and test subsets. This time we will use the `sample()` and `floor()` function to separate training and test dataset. `sample()` is a function to create a set of indicators for row numbers. We can subset the original dataset with random rows using these indicators. `floor()` takes a number *x* and returns the closest integer to *x*

`sample(row, size)`

* row: rows in the dataset that you want to select from. If you want to select all the rows, you can use $nrow(data)$ or $1:nrow(data)$ (for a single number or a vector).
* size: how many rows you want for your subset.

```{r echo=T}
sub<-sample(nrow(google_norm), floor(nrow(google_norm)*0.75))
google_train<-google_norm[sub, ]
google_test<-google_norm[-sub, ]
```

We are good to go! Let's move forward to training phase.

## Step 3 - training a model on the data

Here, we use the function `neuralnet()` in package `neuralnet`. `neuralnet` returns a NN object containing:

* call; the matched call.
* response; extracted from the data argument.
* covariate; the variables extracted from the data argument.
* model.list; a list containing the covariates and the response variables extracted from the formula argument.
* err.fct and act.fct; the error and activation functions.
* net.result; a list containing the overall result of the neural network for every repetition.
* weights; a list containing the fitted weights of the neural network for every repetition.
* result.matrix; a matrix containing the reached threshold, needed steps, error, AIC, BIC, and weights for every repetition. Each column represents one repetition.

`m<-neuralnet(target~predictors, data=mydata, hidden=1)`, where:

* target: variable we want to predict.
* predictors: predictors we want to use. Note that we cannot use "." to denote all the variables in this function. We have to add all predictors one by one to the model.
* data: training dataset.
* hidden: number of hidden nodes that we want to use in the model. By default, it is set to one.

```{r fig.height=15, fig.width=15}
# install.packages("neuralnet")
library(neuralnet)
google_model<-neuralnet(RealEstate~Unemployment+Rental+Mortgage+Jobs+Investing+DJI_Index+StdDJI, data=google_train)
plot(google_model)
```
![](http://wiki.socr.umich.edu/images/e/e5/DSPA_Figs_plot_google_model_Ch6.png)

The above graph shows that we have only one hidden node. `Error` is stand for the sum of squared errors and `Steps` is how many iterations the model has go through. Do note that these outputs could be different when you run exact same codes twice because the weights are randomly generated.

## Step 4 - evaluating model performance

Similar to the `predict()` function that we have mentioned in previous chapters, `compute()` is an alternative method that could help us to generate the model predictions.

`p<-compute(m, test)`

* m: a trained neural networks model.
* test: the test dataset. This dataset should only contain the same type of predictors in the neural network model.

In our model we picked `Unemployment, Rental, Mortgage, Jobs, Investing, DJI_Index, StdDJI` as our predictors. So we need to find these corresponding column numbers in the test dataset (1, 2, 4, 5, 6, 7, 8 respectively).

```{r echo=T}
google_pred<-compute(google_model, google_test[, c(1:2, 4:8)])
pred_results<-google_pred$net.result
cor(pred_results, google_test$RealEstate)
```

As mentioned in [Chapter 9](http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/notes/09_RegressionForecasting.html), we can still use the correlation between predicted results and observed Real Estate Index to evaluate the algorithm. A correlation over 0.9 is very good for real world datasets. Could this be improved further?

## Step 5 - improving model performance

This time we managed to include 4 hidden nodes in the model. Let's see what results we can get from this more complicated model.

```{r echo=T}
google_model2<-neuralnet(RealEstate~Unemployment+Rental+Mortgage+Jobs+Investing+DJI_Index+StdDJI, data=google_train, hidden = 4)
plot(google_model2)
```
![](http://wiki.socr.umich.edu/images/d/d1/DSPA_Figs_plot_google_model2_ch6.png)

Although the graph looks very complicated, we have smaller `Error` or sum of squared errors. Actually, it can be used both for *classification* and *regression*, which we will see in the next part. Let's first try regression.

```{r echo=T}
google_pred2<-compute(google_model2, google_test[, c(1:2, 4:8)])
pred_results2<-google_pred2$net.result
cor(pred_results2, google_test$RealEstate)
```

We get an even higher correlation. This is almost an ideal result! The predicted and observed indexes have a strong linear relationship. Nevertheless, too many hidden nodes might even decrease the correlation between predicted and true values, which will be examined in the practice problems later in this chapter.


## Step 6 - adding additional layers

We observe an even lower `Error` by use three hidden layer with nodes 4,3,3, respectively.

```{r}
google_model2<-neuralnet(RealEstate~Unemployment+Rental+Mortgage+Jobs+Investing+DJI_Index+StdDJI, data=google_train, hidden = c(4,3,3))
google_pred2<-compute(google_model2, google_test[, c(1:2, 4:8)])
pred_results2<-google_pred2$net.result
cor(pred_results2, google_test$RealEstate)
```

# Simple NN demo - learning to compute $\sqrt {\ \ \ }$

This simple example demonstrates the foundation of the neural network prediction of a basic mathematical function: $\sqrt {\ \ \ }: \Re^+ \longrightarrow \Re^+$.

```{r eval=T, warning=F, message=FALSE}
# generate random training data: 1,000 |X_i|, where X_i ~ Uniform (0,10) or perhaps ~ N(0,1)
rand_data <- abs(runif(1000, 0, 100))
 
# create a 2 column data-frame (and_data, sqrt_data)
sqrt_df <- data.frame(rand_data, sqrt_data=sqrt(rand_data)) 
plot(rand_data, sqrt_df$sqrt_data)
```

```{r eval=T, warning=F, message=FALSE}
# Train the neural net
net.sqrt <- neuralnet(sqrt_data ~ rand_data,  sqrt_df, hidden=10, threshold=0.01)
```

```{r eval=T, warning=F, message=FALSE}
 # report the NN
 # print(net.sqrt)
 
 # generate testing data seq(from=0.1, to=N, step=0.1)
 N <- 200.0   # out of range [100: 200] is also included in the testing!
 test_data <- seq(0, N, 0.1); test_data_sqrt <- sqrt(test_data)
 
 # try to predict the square-root values using 10 hidden nodes
 # Compute or predict for test data, test_data
 pred_sqrt <- compute(net.sqrt, test_data)$net.result
 
 # compute uses the trained neural net (net.sqrt), 
 # to estimate the square-roots of the testing data 
 
 # compare real (test_data_sqrt) and NN-predicted (pred_sqrt) square roots of test_data
 plot(pred_sqrt, test_data_sqrt, xlim=c(0, 12), ylim=c(0, 12)); abline(0,1, col="red", lty=2)
 legend("bottomright",  c("Pred vs. Actual SQRT", "Pred=Actual Line"), cex=0.8, lty=c(1,2), lwd=c(2,2),col=c("black","red"))
 compare_df <-data.frame(pred_sqrt, test_data_sqrt); # compare_df
 
 plot(test_data, test_data_sqrt)
 lines(test_data, pred_sqrt, pch=22, col="red", lty=2)
 legend("bottomright",  c("Actual SQRT","Predicted SQRT"), lty=c(1,2), lwd=c(2,2),col=c("black","red"))
``` 

We observe that the NN, `net.sqrt` actually learns and predicts pretty close the complex square root function. Of course, everyone's results may vary as we randomly generate the training data (rand_data) and the NN construction (net.sqrt) is also stochastic.

#Case Study 2: Google Trends and the Stock Market - **Classification**

In practice, NN may be more useful as a classifier. Let's demonstrate this by using again the Stock Market data. We mark the samples according to their `RealEstate`. For those higher than 75% percentile, we give them label 0; For those lower than 0.25 percentile, we label them as 2; Otherwise, label 1. Even in the classification set, response still must be numeric.

```{r}
google_class = google_norm
id1 = which(google_class$RealEstate>quantile(google_class$RealEstate,0.75))
id2 = which(google_class$RealEstate<quantile(google_class$RealEstate,0.25))
id3 = setdiff(1:nrow(google_class),union(id1,id2))
google_class$RealEstate[id1]=0
google_class$RealEstate[id2]=1
google_class$RealEstate[id3]=2
summary(as.factor(google_class$RealEstate))
```

Here, we divide the data to training and testing sets. We need 3 more column indicators, which correspond to the 3 outcomes labels.

```{r}
set.seed(2017)
train = sample(1:nrow(google_class),0.7*nrow(google_class))
google_tr = google_class[train,]
google_ts = google_class[-train,]
train_x = google_tr[,c(1:2,4:8)]
train_y = google_tr[,3]
colnames(train_x)
test_x = google_ts[,c(1:2,4:8)]
test_y = google_ts[3]
train_y_ind = model.matrix(~factor(train_y)-1)
colnames(train_y_ind) = c("High","Median","Low")
train = cbind(train_x, train_y_ind)
```

We use non-linear output and display every 2,000 iterations.

```{r}
nn_single = neuralnet(High+Median+Low~Unemployment+Rental+Mortgage+Jobs+Investing+DJI_Index+StdDJI,
    data = train,
    hidden=4,
    linear.output=FALSE,
    lifesign='full', lifesign.step=2000)
```

Below is the prediction function translating this model to the forecasting results.

```{r}
pred = function(nn, dat) {
    # compute uses the trained neural net (nn=nn_single), and 
    # new testing data (dat=google_ts) to generate predictions (y_hat)
    # compute returns a list containing: 
    #     (1) neurons: a list of the neurons' output for each layer of the neural network, and
    #     (2) net.result: a matrix containing the overall result of the neural network.
    yhat = compute(nn, dat)$net.result
    
    # find the maximum in each row (1) in the net.result matrix
    # to determine the first occurrence of a specific element in each row (1)
    # we can use the apply function with which.max
    yhat = apply(yhat, 1, which.max)-1
    return(yhat)
}

mean(pred(nn_single, google_ts[,c(1:2,4:8)]) != as.factor(google_ts[,3]))
```

Now let's inspect the structure of the Neural Network.

```{r}
plot(nn_single)
```

Similarly, we can change `hidden` to utilize multiple hidden layers, however, a more complicated model won't necessarily guarantee an improved performance.

```{r}
nn_single = neuralnet(High+Median+Low~Unemployment+Rental+Mortgage+Jobs+Investing+DJI_Index+StdDJI,
    data = train,
    hidden=c(4,5),
    linear.output=FALSE,
    lifesign='full', lifesign.step=2000)
mean(pred(nn_single, google_ts[,c(1:2,4:8)]) != as.factor(google_ts[,3]))
```

#  Support Vector Machines (SVM)

## Motivation

Recall that Lazy learning method in [Chapter 6](http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/notes/06_LazyLearning_kNN.html) assign classes using geometrical distances of different features. In multidimensional space (multiple features), it is actually using spheres that centered by training dataset to assign the test dataset with classes. Can we make a fancier shape in the $nD$ space to do the classification?

## Classification with hyperplanes

The easiest shape would be a plane. Support Vector Machine (SVM) can use a hyperplane to separate data into several groups or classes. This is used for datasets that are linearly separable.

Assume that we have only two features, will you use A or B hyperplane to separate the data? Or even another plane C?

```{r echo=F}
A<-c(1, 4, 3, 2, 4, 8, 6, 10, 9)
B<-c(1, 5, 3, 2, 3, 8, 8, 7, 10)
plot(A, B, xlab="", ylab="", pch=16, cex=2)
abline(v=5, col="red", lty=2)
text(5.4, 9, labels="A")
abline(12, -1, col="red", lty=2)
text(6, 5.4, labels="B")
```

### Finding the maximum margin

To answer the above question, we need to search for the **Maximum Margin Hyperplane (MMH)**. That is the hyperplane that creates greatest separation between the two closest observations.

We define support vectors as the points from each class that are closest to MMH. Each class must have at least one observation as support vector.

Using support vectors along is not sufficient for finding the MMH. Although tricky mathematical calculations are involved, the principal of the process is fairly simple. Let's look at linearly separable data and non-linearly separable data individually.

### Linearly separable data

If the dataset is linearly separable, we can find the outer boundaries of our two groups of data points. These boundaries are called convex hull (red lines in the following graph). The MMH (black solid line) is just the line that perpendicular to the shortest line between the two convex hulls.

```{r echo=F}
plot(A, B, xlab="", ylab="", pch=16, cex=2)
segments(1, 1, 4, 5, lwd=1, col = "red")
segments(1, 1, 4, 3, lwd = 1, col = "red")
segments(4, 3, 4, 5, lwd = 1, col = "red")
segments(6, 8, 10, 7, lwd = 1, col = "red")
segments(6, 8, 9, 10, lwd = 1, col = "red")
segments(10, 7, 9, 10, lwd = 1, col = "red")
segments(6, 8, 4, 5, lwd = 1, lty=2)
abline(9.833, -2/3, lwd=2)
```

An alternative way would be picking two parallel planes that can separate the data into two groups while the distance between two planes is as far as possible.

We can use vector notation to mathematically define planes. In *n*-dimensional space, a plane could be expressed by the following equation:
$$\vec{w}\cdot\vec{x}+b=0,$$
where the vectors $\vec{w}$ (weights) and $\vec{x}$ (unknowns) both have *n* coordinates and *b* is a (scalar) constant. 

To clarify this notation let's look at the situation in a 3D space where we can express (embed) 2D Euclidean planes using a point ($(x_o,y_o,z_o)$) and normal-vector ($(a,b,c)$) form. This is just a linear equation, where $d=-(ax_o + by_o + cz_0)$:
$$ax + by + cz + d = 0,$$
or equivalently
$$w_1x_1+w_2x_2+w_3x_3+b=0$$
We can see that it is equivalent to the vector notation

Using the vector notation, we can specify two hyperplanes as follows:
$$\vec{w}\cdot\vec{x}+b\geq+1$$
and
$$\vec{w}\cdot\vec{x}+b\leq-1$$
We require that all of the observations in the first class fall above the first plane and all observations in the other class fall below the second plane.

The distance between two planes is calculated as:
$$\frac{2}{\lVert \vec{w}\rVert}$$
where $\lVert \vec{w}\rVert$ is the Euclidean norm. To maximize the distance we need to minimize the Euclidean norm.

To sum up we are going to find $min\frac{\lVert \vec{w}\rVert}{2}$ with the following constrain:
$$y_i(\vec{w}\cdot\vec{x}-b)\geq1, \, \forall\vec{x}_i$$
where $\forall$ means "for all".

For each nonlinear programming problem, the **primal problem**, there is related nonlinear programming problem, the **Lagrangian dual problem**. Under certain assumptions for convexity and suitable constraints, the [primal and dual problems have equal optimal objective values](http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/notes/21_FunctionOptimization.html#31_dual_problem_optimization). Primal optimization problems are typically described as:
$$\begin{array}{rcl} \min_x{f(x)} \\ 
\text{subject to}  \\
g_i(x) \leq 0\\
h_j(x) = 0 \\
\end{array}.$$

Then the Lagrangian dual problem is defined as a parallel nonlinear programming problem:
$$\begin{array}{rcl} & \min_{u,v}{\theta(u,v)} & \\ 
& \text{subject to} & \\
& u \geq 0\\
\end{array},$$
where 
$$ \theta(u,v)= \inf_{x}{ \left ( f(x)+\displaystyle\sum_i {u_i g_i(x)} +\displaystyle\sum_j {v_j h_j(x)} \right )}.$$

[Chapter 21](http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/notes/21_FunctionOptimization.html#31_dual_problem_optimization) provides additional technical details about optimization duality. 

Suppose the Lagrange primal is $$L_p = \frac{1}{2}||w||^2-\sum_{i=1}^{n}\alpha_i[y_i(w_0+x_i^{T}w)-1],\ \text{where}\ \alpha_i\geq 0.$$

To optimize that objective funciton, we can set the partial derivatives equal to zero:
  $$\frac{\partial}{\partial w}\ :\ w = \sum_{i=1}^{n}\alpha_iy_ix_i$$

  $$\frac{\partial}{\partial b}\ :\ 0 = \sum_{i=1}^{n}\alpha_iy_i.$$

Substituting into the Lagrange primal, we obtain the Lagrange dual:

$$L_D = \sum_{i=1}^{n}\alpha_i - \frac{1}{2} \sum_{i=1}^{n}\alpha_i\alpha_i^{'}y_iy_i^{'}x_i^Tx_i^{'}.$$

We maximize $L_D$ subject to $\alpha_i \geq 0$ and $\sum_{i=1}^{n}\alpha_iy_i =0$.

Since it follows the [Karush-Kuhn-Tucker optimization conditions](https://en.wikipedia.org/wiki/Karush%E2%80%93Kuhn%E2%80%93Tucker_conditions), we have $\hat\alpha[y_i(\hat{b}+x_i^T\hat{w})-1]=0.$

Which implies that if $y_i\hat{f}(x_i)>1$, then $\hat{\alpha}_i=0$.

The **support** of a function ($f$) is the smallest subset of the domain containing only arguments ($x$) which are not mapped to zero ($f(x)\not=0$). In our case, the solution $\hat{w}$ is defined in terms of a linear combination of the **support points**:

$$\hat{f}(x)=w^Tx = w = \sum_{i=1}^{n}\alpha_iy_ix_i. $$

That's where the name of Support Vector Machines (SVM) comes from.

### Non-linearly separable data

For non-linearly separable data, we need to use a small trick. Still, we use a plane but allow some of the points to be misclassified into the wrong class. To penalize for that, we add a cost term after the Euclidean norm function that we need to minimize.

Therefore, the solution will optimize the following objective (cost) function:
$$\min \left (\frac{\lVert \vec{w}\rVert}{2} \right )+C\sum_{i=1}^{n} \xi_i$$
$$\text{subject to}$$
$$y_i(\vec{w}\cdot\vec{x}-b)\geq1, \, \forall\vec{x}_i, \, \xi_i\geq0,$$
where $C$ is the cost and $\xi_i$ is the distance between the misclassified observation *i* and the plane.

We have Lagrange primal problem: 
$$L_p = \frac{1}{2}||w||^2 + C\sum_{i=1}^{n}\xi_i-\sum_{i=1}^{n}\alpha_i[y_i(b+x_i^{T}w)-(1-\xi_i)] - \sum_{i=1}^{n}\gamma_i\xi_i,$$
where $$\alpha_i,\gamma_i \geq 0.$$

Similar to what we did above for the separable case, we can use the derivatives of the primal problem to solve the dual problem. 

Notice the inner product in the final expression. We can replace this inner product with a kernel function that maps the feature space into a higher dimensional space (e.g., using a polynomial kernel) or an infinite dimensional space (e.g., using a Gaussian kernel).

### Using kernels for non-linear spaces

An alternative way to solve for the non-linear separable is called the kernel trick. That is to add some dimensions (or features) to make these non-linear separable data to be separable in a higher dimensional space.

How can we do that? We transform our data using kernel functions. A general form for kernel functions would be:
$$K(\vec{x_i}, \vec{x_j})=\phi(\vec{x_i})\cdot\phi(\vec{x_j}),$$
where $\phi$ is a mapping of the data into another space.

The linear kernel would be the simplest one that is just the dot product of the features.
$$K(\vec{x_i}, \vec{x_j})=\vec{x_i}\cdot\vec{x_j}.$$
The polynomial kernel of degree *d* transform the data by adding a simple non-linear transformation of the data.
$$K(\vec{x_i}, \vec{x_j})=(\vec{x_i}\cdot\vec{x_j}+1)^d.$$
The sigmoid kernel is very similar to neural network. It uses a sigmoid activation function.
$$K(\vec{x_i}, \vec{x_j})=tanh(k\vec{x_i}\cdot\vec{x_j}-\delta).$$
The Gaussian RBF kernel is similar to RBF neural network and is a good place to start investigating a dataset.
$$K(\vec{x_i}, \vec{x_j})=exp \left (\frac{-\lVert \vec{x_i}-\vec{x_j}\rVert^2}{2\sigma^2}\right ) .$$

# Case Study 2: Optical Character Recognition (OCR)

This example illustrates management and transferring of handwritten notes (text) and converting it to typeset or printed text representing the characters in the original notes (unstructured image data).

* Protocol:
 * Divides the image (typically optical image of handwritten notes on paper) into a fine grid where each cell contains 1 glyph (symbol, letter, number).
 * Match the glyph in each cell to 1 of the possible characters in a dictionary.
 * Combine individual characters together into words to reconstitute the digital representation of the optical image of the handwritten notes.
 
In this example, we use an optical document image (data) that has already been pre-partitioned into rectangular grid cells containing 1 character of the [26 English letters, A through Z](http://wiki.socr.umich.edu/index.php/SOCR_LetterFrequencyData).

The resulting gridded dataset is distributed by the [UCI Machine Learning Data Repository](http://archive.ics.uci.edu/ml). The [dataset contains 20, 000 examples of 26 English capital letters printed using 20 different randomly reshaped and morphed fonts](https://umich.instructure.com/courses/38100/files/folder/Case_Studies/16_HandwrittenLetters).

![Example of the preprocessed gridded handwritten letters](http://wiki.socr.umich.edu/images/2/2d/SOCR_Data_Dinov_EnglishLetterFrequency_Fig2.png)

## Step 1: Prepare and explore the data
```{r eval=T}
# read in data and examine structure
hand_letters <- read.csv("https://umich.instructure.com/files/2837863/download?download_frd=1", header = T)
str(hand_letters)

# divide into training (3/4) and testing (1/4) data
hand_letters_train <- hand_letters[1:15000, ]
hand_letters_test  <- hand_letters[15001:20000, ]
```

## Step 3: Training an SVM model

We can specify `vanilladot` as a linear kernel, or alternatively:


* `rbfdot` Radial Basis kernel i.e, "Gaussian" 
* `polydot` Polynomial kernel
* `tanhdot` Hyperbolic tangent kernel
* `laplacedot` Laplacian kernel
* `besseldot` Bessel kernel
* `anovadot` ANOVA RBF kernel
*  `splinedot` Spline kernel 
* `stringdot` String kernel


```{r eval=T}
# begin by training a simple linear SVM
library(kernlab)
set.seed(123)
hand_letter_classifier <- ksvm(letter ~ ., data = hand_letters_train, kernel = "vanilladot")

# look at basic information about the model
hand_letter_classifier
```

## Step 4: Evaluating model performance

```{r eval=T}
# predictions on testing dataset
hand_letter_predictions <- predict(hand_letter_classifier, hand_letters_test)

head(hand_letter_predictions)

table(hand_letter_predictions, hand_letters_test$letter)

# look only at agreement vs. non-agreement
# construct a vector of TRUE/FALSE indicating correct/incorrect predictions
agreement <- hand_letter_predictions == hand_letters_test$letter # check if characters agree
table(agreement)
prop.table(table(agreement))
```

## Step 5: Improving model performance

Replacing the `vanilladot` linear kernel with `rbfdot` Radial Basis Function kernel, i.e., "Gaussian" kernel may improve the OCR prediction.
    
```{r eval=T}
hand_letter_classifier_rbf <- ksvm(letter ~ ., data = hand_letters_train, kernel = "rbfdot")
hand_letter_predictions_rbf <- predict(hand_letter_classifier_rbf, hand_letters_test)

agreement_rbf <- hand_letter_predictions_rbf == hand_letters_test$letter
table(agreement_rbf)
prop.table(table(agreement_rbf))
```

Note the improvement of automated (SVM) classification accuracy ($0.928$) for `rbfdot` compared to the previous (`vanilladot`) result ($0.844$).

# Case Study 3: Iris Flowers

Let's have another look at the *iris data* that we saw in [Chapter 2](http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/notes/02_ManagingData.html).

## Step 1 - collecting data

SVM require all features to be numeric and each feature has to be scaled into a relative small interval. We are using the Edgar Anderson's Iris Data in R for this case study. This dataset measures the length and width of sepals and petals from three Iris flower species.

## Step 2 - exploring and preparing the data

Let's load the data first. In this case study we want to explore the variable `Species`.

```{r echo=T}
data(iris)
str(iris)
table(iris$Species)
```

The data looks good. However, recall that we need a fairly normalized data. We could normalize the data by hand. Luckily, the R package we are going to use will normalized the dataset automatically.

Now we can separate the training and test dataset using 75%-25% rule.

```{r echo=T}
sub<-sample(nrow(iris), floor(nrow(iris)*0.75))
iris_train<-iris[sub, ]
iris_test<-iris[-sub, ]
```

Let's first try a toy (iris data) example.

```{r,fig.height=6,fig.align='center',fig.width=8}
require(e1071)
iris.svm_1 <- svm(Species~Petal.Length+Petal.Width, data=iris_train,
                kernel="linear", cost=1)
iris.svm_2 <- svm(Species~Petal.Length+Petal.Width, data=iris_train,
                kernel="radial", cost=1)
par(mfrow=c(2,1))
plot(iris.svm_1, iris[,c(5,3,4)]); legend("center", "Linear")
plot(iris.svm_2, iris[,c(5,3,4)]); legend("center", "Radial")
```

## Step 3 - training a model on the data

We are going to use `kernlab` for this case study. However other packages like `e1071` and `klaR` are available if you are quite familiar with C++.

Let's break down the function `ksvm()`

`m<-ksvm(target~predictors, data=mydata, kernel="rbfdot", c=1)`

* target: the outcome variable that we want to predict.
* predictors: features that the prediction based on. In this function we can use the "." to represent all the variables in the dataset again.
* data: the training dataset that the *target* and *predictors* can be find.
* kernel: is the kernel mapping we want to use. By default it is the radio basis function (`rbfdot`).
* C is a number that specifies the cost of misclassification.

Let's in stall the package and play with the data now.

```{r echo=T}
# install.packages("kernlab")
library(kernlab)
iris_clas<-ksvm(Species~., data=iris_train, kernel="vanilladot")
iris_clas
```

Here, we used all the variables other than the `Species` in the dataset as predictors. We also used kernel `vanilladot` that is the linear kernel in this model. We get a training error less than 0.02. 

## Step 4 - evaluating model performance

Our old friend `predict()` function is used again to make predictions. Here we have a factor outcome, so we need the command `table()` to show us how well do the predictions and actual data match.

```{r echo=T}
iris.pred<-predict(iris_clas, iris_test)
table(iris.pred, iris_test$Species)
```

We can see that only 1-2 cases of may be misclassifies as Iris versicolor. The species of the majority of the flowers are all correctly identified.

To see the results more clearly, we can use the proportional table to show the agreements of the categories.

```{r echo=T}
agreement<-iris.pred==iris_test$Species
prop.table(table(agreement))
```

Here `==` means "equal to". Over 90% of predictions are correct. Nevertheless, is there any chance that we can improve the outcome? What if we try a Gaussian kernel?

## Step 5 - RBF kernel function

Linear kernel is the simplest one but usually not the best one. Let's try the **RBF (Radial Basis "Gaussian" Function)** kernel instead.

```{r echo=T}
iris_clas1<-ksvm(Species~., data=iris_train, kernel="rbfdot")
iris_clas1
iris.pred1<-predict(iris_clas1, iris_test)
table(iris.pred1, iris_test$Species)
agreement<-iris.pred1==iris_test$Species
prop.table(table(agreement))
```

Unfortunately, the model performance is actually worse than the previous one (you might get slightly different results). This is because this Iris dataset has a linear feature. In practice, we could try alternative kernel functions and see which one fits the dataset the best.

## Parameter Tuning

We can tune the SVM using the `tune.svm` function in the package `e1071`.

```{r}
costs = exp(-5:8)
tune.svm(Species~., kernel = "radial", data = iris_train, cost = costs)
```

Further, we can draw a **cv** plot to gauge the model performance:

```{r}
set.seed(2017)
require(sparsediscrim); require (reshape); require(ggplot2)
folds = cv_partition(iris$Species, num_folds = 5)
train_cv_error_svm = function(costC) {
  #Train
  ir.svm = svm(Species~., data=iris,
                 kernel="radial", cost=costC)
  train_error = sum(ir.svm$fitted != iris$Species) / nrow(iris)
  #Test
  test_error = sum(predict(ir.svm, iris_test) != iris_test$Species) / nrow(iris_test)
  #CV error
  ire.cverr = sapply(folds, function(fold) {
    svmcv = svm(Species~.,data = iris, kernel="radial", cost=costC, subset = fold$training)
    svmpred = predict(svmcv, iris[fold$test,])
    return(sum(svmpred != iris$Species[fold$test]) / length(fold$test))
  })
  cv_error = mean(ire.cverr)
  return(c(train_error, cv_error, test_error))
}

costs = exp(-5:8)
ir_cost_errors = sapply(costs, function(cost) train_cv_error_svm(cost))
df_errs = data.frame(t(ir_cost_errors), costs)
colnames(df_errs) = c('Train', 'CV', 'Test', 'Logcost')
dataL <- melt(df_errs, id="Logcost")
ggplot(dataL, aes_string(x="Logcost", y="value", colour="variable",
                         group="variable", linetype="variable", shape="variable")) +
  geom_line(size=1) + labs(x = "Cost",
                           y = "Classification error",
                           colour="",group="",
                           linetype="",shape="") + scale_x_log10()
```

## Improving the performance of Gaussian kernels

Now, let's attempt to improve the performance of a Gaussian kernel by tuning:

```{r}
set.seed(2020)
gammas = exp(-5:5)
tune_g = tune.svm(Species~., kernel = "radial", data = iris_train, cost = costs, gamma = gammas)
tune_g
```

We observe that the model achieves a better prediction now.

```{r}
iris.svm_g <- svm(Species~., data=iris_train,
                kernel="radial", gamma=0.0183, cost=20)
table(iris_test$Species, predict(iris.svm_g, iris_test))
agreement<-predict(iris.svm_g, iris_test)==iris_test$Species
prop.table(table(agreement))
```

[Chapter 22](http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/notes/22_DeepLearning.html) provides more details about neural networks and deep learning.

# Practice

## Problem 1: Google Trends and the Stock Market

Use the Google trend data. Fit a neural network model with the same training data as case study 1. This time use `Investing` as target and `Unemployment, Rental, RealEstate, Mortgage, Jobs, DJI_Index, StdDJI` as predictors. Use 3 hidden nodes. 

**Note**: remember to change the columns you want to include in the test dataset when predicting.

The following number is the correlation between predicted and observed values.

```{r echo=FALSE, fig.height=15, fig.width=15}
google_model3<-neuralnet(Investing~Unemployment+Rental+RealEstate+Mortgage+Jobs+DJI_Index+StdDJI, data=google_train, hidden = 3)
plot(google_model3)
google_pred3<-compute(google_model3, google_test[, c(1:5, 7:8)])
pred_results3<-google_pred3$net.result
cor(pred_results3, google_test$Investing)
```

You might get slightly different results since the weights are generated randomly.

## Problem 2: Quality of Life and Chronic Disease

Use the same data we showed in [Chapter 8](http://www.socr.umich.edu/people/dinov/2017/Spring/DSPA_HS650/notes/08_DecisionTreeClass.html), [Quality of life and chronic disease](https://umich.instructure.com/files/481332/download?download_frd=1) and the corresponding [meta-data doc](https://umich.instructure.com/files/399150/download?download_frd=1&verifier=aq3okIt8HwhhGPPWmbz7yd7vYGOKILe0Wvu13E6d). 

Let's load the data first. In this case study, we want to use the variable `CHARLSONSCORE` as our target variable.

```{r echo=T}
qol<-read.csv("https://umich.instructure.com/files/481332/download?download_frd=1")
str(qol)
```

Delete the first two columns (we don't need ID variables) and rows that have missing values in `CHARLSONSCORE`(where `CHARLSONSCORE`equals "-9")
`!qol$CHARLSONSCORE==-9` means we want all the rows that have CHARLSONSCORE not equal to -9. The exclamation sign indicates "exclude". Also, we need to convert our categorical variable `CHARLSONSCORE` into a factor.

```{r echo=T}
qol<-qol[!qol$CHARLSONSCORE==-9 , -c(1, 2)]
qol$CHARLSONSCORE<-as.factor(qol$CHARLSONSCORE)
str(qol)
```

Now the dataset is ready. First, separate the dataset into training and test datasets using 75%-25% rule. Then, build a SVM model using all other variables in the dataset to be predictor variables. Try to add different cost of misclassification to the model. Rather than the default `C=1` we use `C=2` and `C=3`. See how the model behaves. Here we utilize the radio basis kernel.

Output for `C=2`
```{r echo=F}
sub<-sample(nrow(qol), floor(nrow(qol)*0.75))
qol_train<-qol[sub, ]
qol_test<-qol[-sub, ]
qol_clas2<-ksvm(CHARLSONSCORE~., data=qol_train, kernel="rbfdot", C=2)
qol_clas2
qol.pred2<-predict(qol_clas2, qol_test)
table(qol.pred2, qol_test$CHARLSONSCORE)
agreement<-qol.pred2==qol_test$CHARLSONSCORE
prop.table(table(agreement))
```
Output for `C=3`
```{r echo=F}
qol_clas3<-ksvm(CHARLSONSCORE~., data=qol_train, kernel="rbfdot", C=3)
qol_clas3
qol.pred3<-predict(qol_clas3, qol_test)
table(qol.pred3, qol_test$CHARLSONSCORE)
agreement<-qol.pred3==qol_test$CHARLSONSCORE
prop.table(table(agreement))
```

Can you reproduce or improve these results?

# Appendix
Below is some additional `R` code demonstrating various results reported in this Chapter.

```{r eval=F}
#Picture 1
x<-runif(1000, -10, 10)
y<-ifelse(x>=0, 1, 0)
plot(x, y, xlab = "Sum of input signals", ylab = "Output signal", main = "Threshold")
abline(v=0, lty=2)
#Picture 2
x<-runif(100000, -10, 10)
y<-1/(1+exp(-x))
plot(x, y, xlab = "Sum of input signals", ylab = "Output signal", main = "Sigmoid")
#Picture 3
x<-runif(100000, -10, 10)
y1<-x
y2<-ifelse(x<=-5, -5, ifelse(x>=5, 5, x))
y3<-(exp(x)-exp(-x))/(exp(x)+exp(-x))
y4<-exp(-x^2/2)
par(mfrow=c(2, 2))
plot(x, y1, main="Linear", xlab="", ylab="")
plot(x, y2, main="Saturated Linear", xlab="", ylab="")
plot(x, y3, main="Hyperbolic tangent", xlab="", ylab="")
plot(x, y4, main = "Gaussian", xlab="", ylab="")
#Picture 4
A<-c(1, 4, 3, 2, 4, 8, 6, 10, 9)
B<-c(1, 5, 3, 2, 3, 8, 8, 7, 10)
plot(A, B, xlab="", ylab="", pch=16, cex=2)
abline(v=5, col="red", lty=2)
text(5.4, 9, labels="A")
abline(12, -1, col="red", lty=2)
text(6, 5.4, labels="B")
#Picture 5
plot(A, B, xlab="", ylab="", pch=16, cex=2)
segments(1, 1, 4, 5, lwd=1, col = "red")
segments(1, 1, 4, 3, lwd = 1, col = "red")
segments(4, 3, 4, 5, lwd = 1, col = "red")
segments(6, 8, 10, 7, lwd = 1, col = "red")
segments(6, 8, 9, 10, lwd = 1, col = "red")
segments(10, 7, 9, 10, lwd = 1, col = "red")
segments(6, 8, 4, 5, lwd = 1, lty=2)
abline(9.833, -2/3, lwd=2)
```

Try to replicate these results with [other data from the list of our Case-Studies](https://umich.instructure.com/courses/38100/files/).
