<jupyter_start><jupyter_text>Stats216v: Statistical LearningStanford UniversitySummer 2017<NAME> (<EMAIL>)  Q1.In this problem we investigate a simple example of hierarchical clustering. Problem 4 will use this method on a real data set.Suppose that we have four observations, for which we compute a dissimilarity matrix, given by$\begin{bmatrix}   & 0.3 & 0.4 & 0.7 \\  0.3  & & 0.5 & 0.8 \\  0.4 & 0.5 & & 0.45 \\  0.7 & 0.8 & 0.45 &\end{bmatrix}$For instance, the dissimilarity between the first and second observations is 0.3, and thedissimilarity between the second and fourth observation is 0.8. Gyu-Ho's Answer: I will use *Algorithm 10.2 Hierarchical Clustering at page 395*. (a) Based on this dissimilarity matrix, sketch the dendrogram that results from hierarchically clustering these four observations using **complete linkage**. Indicate on the plot the height at which each fusion occurs and the observations corresponding to each leaf of the dendrogram.<jupyter_code>matrix(c(0, 0.3, 0.4, 0.7,
         0.3, 0, 0.5, 0.8,
         0.4, 0.5, 0.0, 0.45,
         0.7, 0.8, 0.45, 0.0), nrow=4)

# 0.3 is the minimum dissimilarity
# complete linkage records the largest
# thus, fuse first and second observations
# after fuse, new dissimilarity matrix
matrix(c(0,   0.5, 0.8,
         0.5, 0.0, 0.45,
         0.8, 0.45, 0.0), nrow=3)

# now 0.45 is the minimum dissimilarity
# thus, fuse third and fourth observations
# after fuse, new dissimilarity matrix
matrix(c(0, 0.8,
         0.8, 0.0), nrow=2)
# lastly fuse clusters (1,2) and (3,4)
dend = as.dist(matrix(c(0, 0.3, 0.4, 0.7,
                        0.3, 0, 0.5, 0.8,
                        0.4, 0.5, 0.0, 0.45,
                        0.7, 0.8, 0.45, 0.0), nrow=4))

plot(hclust(dend,method="complete"))<jupyter_output><empty_output><jupyter_text>(b) Repeat part (a) using **single linkage** clustering.<jupyter_code>matrix(c(0, 0.3, 0.4, 0.7,
         0.3, 0, 0.5, 0.8,
         0.4, 0.5, 0.0, 0.45,
         0.7, 0.8, 0.45, 0.0), nrow=4)

# 0.3 is the minimum dissimilarity
# single linkage records the smallest
# thus, fuse first and second observations
# after fuse, new dissimilarity matrix
matrix(c(0,   0.4, 0.7,
         0.4, 0.0, 0.45,
         0.7, 0.45, 0.0), nrow=3)

# now 0.4 is the minimum dissimilarity
# thus, fuse third and fourth observations
# after fuse, new dissimilarity matrix
matrix(c(0, 0.45,
         0.45, 0.0), nrow=2)
# lastly fuse clusters ((1,2),3) and (4)
dend = as.dist(matrix(c(0, 0.3, 0.4, 0.7,
                        0.3, 0, 0.5, 0.8,
                        0.4, 0.5, 0.0, 0.45,
                        0.7, 0.8, 0.45, 0.0), nrow=4))

plot(hclust(dend,method="single"))<jupyter_output><empty_output><jupyter_text>(c) Suppose that we cut the dendrogram from (a) in such a way that results in twoclusters. Which observations are in each cluster? Gyu-Ho's Answer: Two clusters of (1,2) and (3,4). (d) Suppose that we cut the dendrogram from (b) in such a way that results in two clusters. Which observations are in each cluster? Gyu-Ho's Answer: Two clusters of ((1,2),3) and (4).  Q2.Here we explore the maximal margin classifier on a toy data set. (a) We are given $n = 7$ observations in $p = 2$ dimensions. For each observation, there is an associated class label. Sketch the observations.| Obs. | $X_1$ | $X_2$ | $Y$ ||---|---|---|---||1|13.58|17.50|blue||2|9.60|9.60|blue||3|17.60|17.60|blue||4|5.56|17.60|blue||5|9.60|5.60|red||6|17.60|13.60|red||7|17.58|5.60|red|<jupyter_code>x1=c(13.58,9.60,17.60,5.56,9.60,17.60,17.58)
x2=c(17.50,9.60,17.60,17.60,5.60,13.60,5.60)
colors=c("blue","blue","blue","blue","red","red","red")
plot(x1,x2,col=colors)<jupyter_output><empty_output><jupyter_text>(b) Sketch the optimal separating hyperplane, and provide the equation for this hyperplane (of the form $\beta_0 + \beta_1 X_1 + \beta_2 X_2 = 0$).<jupyter_code>plot(x1,x2,col=colors)

# between (9.60,9.60),(17.60,17.60)
#     and (9.60,5.60),(17.60,13.60)
#
# line passes (9.6,7.6) and (17.6,15.6)
abline(-2,1) # y=x-2<jupyter_output><empty_output><jupyter_text>(c) Describe the classification rule for the maximal margin classifier. It should be something along the lines of “Classify to Red if $\beta_0 + \beta_1 X_1 + \beta_2 X_2 > 0$, and classify to Green otherwise.” Provide the values for $\beta_0$, $\beta_1$, and $\beta_2$. Gyu-Ho's Answer: Classify to Blue if $X_1 - X_2 - 2 < 0$, and classify to Red otherwise. (d) On your sketch, indicate the margin for the maximal margin hyperplane. How wide is the margin?<jupyter_code>plot(x1,x2,col=colors)
abline(-2,1) # y=x-2
abline(0,1,lty=2) # y=x
abline(-4,1,lty=2) # y=x-4<jupyter_output><empty_output><jupyter_text>Gyu-Ho's Answer: Margin is the smallest distance from the observations to the hyperplane. It's $\sqrt{2}$. (e) Indicate the support vectors for the maximal margin classifier. Gyu-Ho's Answer: Support vectors are (9.60,9.60),(17.60,17.60) and (9.60,5.60),(17.60,13.60). (f) Argue that a slight movement of the seventh observation would not affect the maximal margin hyperplane. Gyu-Ho's Answer: Seventh observation is located at right-bottom (17.58,5.6) and not a support vector. It is far away from the hyperplane, so moving the seventh observation would not change the maximal margin hyperplane. (g) Sketch a hyperplane that is not the optimal separating hyperplane, and provide the equation for this hyperplane. Gyu-Ho's Answer: $X_1 - X_2 - 1 = 0$ is not the optimal separating hyperplane as shown below.<jupyter_code>plot(x1,x2,col=colors)
abline(-1,1)<jupyter_output><empty_output><jupyter_text>(h) Draw an additional observation on the plot so that the two classes are no longer separable by a hyperplane.<jupyter_code>plot(x1,x2,col=colors)
points(c(12),c(6),col=c("blue")) # additional observation
LoadLibraries = function() {
    install.packages("ISLR")
    library(ISLR)
    install.packages("e1071")
    library(e1071)
    install.packages("dendextend")
    library(dendextend)
    install.packages("pls")
    library(pls)
    print("Libraries have been loaded!")
}

LoadLibraries()<jupyter_output>Updating HTML index of packages in '.Library'
Making 'packages.html' ... done
Updating HTML index of packages in '.Library'
Making 'packages.html' ... done
Updating HTML index of packages in '.Library'
Making 'packages.html' ... done
Updating HTML index of packages in '.Library'
Making 'packages.html' ... done

Attaching package: 'pls'

The following object is masked from 'package:stats':

    loadings<jupyter_text>Q3.This problem involves the OJ data set which is part of the ISLR package. (a) Create a training set containing a random sample of 535 observations, and a test set containing the remaining observations. Use the commands `set.seed(2017); train = sample(1:nrow(OJ), 535); test = setdiff(1:nrow(OJ), train)`.<jupyter_code>attach(OJ)
head(OJ)
set.seed(2017)

train=sample(1:nrow(OJ), 535)
test=setdiff(1:nrow(OJ), train)<jupyter_output><empty_output><jupyter_text>(b) Fit a (linear) support vector classifier to the training data using `cost=1`, with Purchase as the response and the other variables as predictors. Use the `summary()` function to produce summary statistics about the SVM, and describe the results obtained.<jupyter_code>OJ.train=OJ[train,]
OJ.test=OJ[-train,]
svmfit=svm(Purchase~., data=OJ.train, kernel="linear", cost=1)
summary(svmfit)<jupyter_output><empty_output><jupyter_text>Gyu-Ho's Answer: This tells us that a linear kernel was used with cost parameter 1. And there are 227 support vectors: 115 in one class, 112 in the other class. Two classes are `CH` and `MM`. (c) What are the training and test error rates?<jupyter_code>OJ.train.pred=predict(svmfit)
OJ.test.pred=predict(svmfit, newdata=OJ.test)
OJ.train.table=table(Predict=OJ.train.pred, Truth=OJ.train$Purchase)
OJ.test.table=table(Predict=OJ.test.pred, Truth=OJ.test$Purchase)
OJ.train.table
OJ.test.table
# training error rate
(OJ.train.table[1,2]+OJ.train.table[2,1])/sum(OJ.train.table)
# test error rate
(OJ.test.table[1,2]+OJ.test.table[2,1])/sum(OJ.test.table)<jupyter_output><empty_output><jupyter_text>(d) Use the `tune()` function to select an optimal cost. Consider values in the range 0.01 to 10.<jupyter_code># cost range
costs=seq(.01, 10, 100)
svmfit.tune=tune(svm,
                 Purchase~.,
                 data=OJ,
                 kernel="linear",
                 ranges=list(cost=costs),
                )

svmfit.tune$best.model<jupyter_output><empty_output><jupyter_text>Gyu-Ho's Answer: 0.01 is already the optimal cost. (e) Compute the training and test error rates using this new value for cost.<jupyter_code>OJ.train=OJ[train,]
OJ.test=OJ[-train,]
svmfit=svm(Purchase~., data=OJ.train, kernel="linear", cost=0.01)

OJ.train.pred=predict(svmfit)
OJ.test.pred=predict(svmfit, newdata=OJ.test)
OJ.train.table=table(Predict=OJ.train.pred, Truth=OJ.train$Purchase)
OJ.test.table=table(Predict=OJ.test.pred, Truth=OJ.test$Purchase)
# training error rate
(OJ.train.table[1,2]+OJ.train.table[2,1])/sum(OJ.train.table)
# test error rate
(OJ.test.table[1,2]+OJ.test.table[2,1])/sum(OJ.test.table)<jupyter_output><empty_output><jupyter_text>(f) Repeat parts (b) through (e) using a support vector machine with a radial kernel. Use the default value for gamma.<jupyter_code>costs=seq(.01, 10, 100)
svmfit.tune.radial=tune(svm,
                        Purchase~.,
                        data=OJ,
                        kernel="radial",
                        ranges=list(cost=costs),
                       )

svmfit.tune.radial$best.model
OJ.radial.train.pred=predict(svmfit.tune.radial$best.model, OJ.train)
OJ.radial.test.pred=predict(svmfit.tune.radial$best.model, OJ.test)

OJ.radial.train.table=table(
    Predict=OJ.radial.train.pred,
    Truth=OJ.train$Purchase)

OJ.radial.test.table=table(
    Predict=OJ.radial.test.pred,
    Truth=OJ.train$Purchase)
OJ.radial.train.table
OJ.radial.test.table # radial dose not fit well in this case<jupyter_output><empty_output><jupyter_text>(g) Repeat parts (b) through (e) using a support vector machine with a polynomial kernel of degree 2.<jupyter_code>costs=seq(.01, 10, 100)
svmfit.tune.poly=tune(svm,
                      Purchase~.,
                      data=OJ,
                      kernel="polynomial",
                      degree=2,
                      ranges=list(cost=costs),
                     )

OJ.poly.train.pred=predict(svmfit.tune.poly$best.model, OJ.train)
OJ.poly.test.pred=predict(svmfit.tune.poly$best.model, OJ.test)

OJ.poly.train.table=table(
    Predict=OJ.poly.train.pred,
    Truth=OJ.train$Purchase)

OJ.poly.test.table=table(
    Predict=OJ.poly.test.pred,
    Truth=OJ.train$Purchase)

# polynomial approach does not fit well either
OJ.poly.train.table
OJ.poly.test.table<jupyter_output><empty_output><jupyter_text>(h) Repeat parts (b) through (e) using a linear support vector machine, applied to an expanded feature set consisting of linear and all possible quadratic terms for the predictors. How does this compare to the polynomial kernel both conceptually and in terms of the results for this problem?<jupyter_code>names(OJ)
costs=seq(.01, 10, 100)
svmfit.tune.all=tune(svm,
                      Purchase~poly(WeekofPurchase,2)
                     +poly(StoreID,2)
                     +poly(PriceCH,2)
                     +poly(PriceMM,2)
                     +poly(DiscCH,2)
                     +poly(DiscMM,2)
                     +poly(LoyalCH,2)
                     +poly(SalePriceMM,2)
                     +poly(SalePriceCH,2)
                     +poly(PriceDiff,2)
                     +poly(PctDiscMM,2)
                     +poly(PctDiscCH,2)
                     +poly(ListPriceDiff,2)
                     +poly(STORE,2),
                      data=OJ,
                      kernel="linear",
                      degree=2,
                      ranges=list(cost=costs),
                     )

OJ.all.train.pred=predict(svmfit.tune.all$best.model, OJ.train)
OJ.all.test.pred=predict(svmfit.tune.all$best.model, OJ.test)

OJ.all.train.table=table(
    Predict=OJ.all.train.pred,
    Truth=OJ.train$Purchase)

OJ.all.test.table=table(
    Predict=OJ.all.test.pred,
    Truth=OJ.train$Purchase)

# linear approach with all possible quadratic terms does not fit well either
OJ.all.train.table
OJ.all.test.table<jupyter_output><empty_output><jupyter_text>Gyu-Ho's Answer: Linear kernel above just adds quadratic terms with existing linear ones. However, decision boundary is still the linear. Polynomial kernel uses polynomial decision boundary. Since linear model performs well in this data, linear kernel with quadratic terms still performs better than polynomial kernel with polynomial decision boundary. (i) Overall, which approach seems to give the best results on this data? Gyu-Ho's Answer: Linear support vector classifier gives the best results (lowest test error rate).  Q4.In this problem we will explore K-means and hierarchical clustering on a wheat seed data set. The data file `SeedData.csv` is available on the class website. There are `n = 210` seeds with `p = 7` real valued features for each observation. The data set is from <NAME>. (2013) from the UCI Machine Learning Repository 1. You can find a full description of the data at http://archive.ics.uci.edu/ml/datasets/seeds.<jupyter_code>seeds=read.table("hw-4-seeds_dataset.txt")
head(seeds)

seeds.labels=seeds[,8]
seeds.data=seeds[,1:7]<jupyter_output><empty_output><jupyter_text>(a) Set `set.seed(2017)`. Using the `kmeans()` function, perform K-means clustering (using all 7 features) with `K = 3` and 20 random starting positions. Plot the data in terms of the variables in the 2nd and 6th columns, with the points colored according to the clusters you just obtained.<jupyter_code>set.seed(2017)
km.out=kmeans(seeds.data,3,nstart=20)
km.out$cluster
plot(seeds.data[,2:6], col=(km.out$cluster+1), main="K-means Clustering with K=3", pch=20, cex=2)<jupyter_output><empty_output><jupyter_text>(b) Using the `hclust()` function, perform single linkage hierarchical clustering on the data. Plot the corresponding dendrogram.<jupyter_code>hc.single=hclust(dist(seeds.data), method="single")
plot(hc.single, main="Single Linkage")<jupyter_output><empty_output><jupyter_text>(c) Using the `cutree()` function, make a cut on the tree that splits the data into 3 clusters. Plot the data again in terms of the variables in the 2nd and 6th columns, coloring the points by the clusters you just obtained. How does it compare to the result you got using k-means clustering? Gyu-Ho's Answer: Most observations are classified as 1, while K-means clustering classifies them quite equally among 1, 2, and 3.<jupyter_code>cutree(hc.single, 3)
dend=as.dendrogram(hclust(dist(seeds.data[,2:6]), method="single"))
dend.color=color_branches(dend, k=3)
plot(dend.color, main="Hierarchical Clustering with Single Linkage")<jupyter_output><empty_output><jupyter_text>(d) Using the `hclust()` function, perform complete linkage hierarchical clustering on the data. Plot the corresponding dendrogram. What do you see? If we wanted to cluster our data into 3 clusters, do you prefer single linkage clustering or complete linkage clustering in this example? Why? Gyu-Ho's Answer: I prefer **complete linkage** because it builds more balanced dendrogram in this data.<jupyter_code>hc.complete=hclust(dist(seeds.data), method="complete")
cutree(hc.complete, 3)
dend=as.dendrogram(hc.complete, method="complete")
dend.color=color_branches(dend, k=3)
plot(dend.color, main="Hierarchical Clustering with Complete Linkage")<jupyter_output><empty_output><jupyter_text>(e) It turns out that there are three different types of wheat seeds: Kama, Rosa, and Canadian. If you had access to the class labels, briefly describe one way to build a model for classifying wheat seeds. Limit your answer to 1 sentence. Gyu-Ho's Answer: Now we have response variable, *wheat seeds*, upon which we can train **logistic regression** model to predict seed class label, thus it's now a supervised learning.  Q5.This problem uses the ALS data set. Recall that in previous problem sets we tried to predict the rate of progression of ALS in patients using linear regression, the Lasso, decision trees, boosting and random forests, and saw varying values of the test RMSE. Now we will see if we can improve on these results by first preprocessing the data in an unsupervised way.<jupyter_code># HW3 Q5
loaded = load("hw-1-als.RData")
head(loaded) # 'train.X' 'train.y' 'test.X' 'test.y'

rData = get(loaded)
head(rData)

ls() # 'Elite' 'college' 'isfar' 'loaded' 'op' 'rData' 'test.X' 'test.y' 'train.X' 'train.y'
head(test.X)<jupyter_output><empty_output><jupyter_text>(a) Use the `prcomp()` function to run a principal components analysis on the training set.<jupyter_code>set.seed(2017)
train.X=train.X
test.X=test.X
pr.out=prcomp(train.X, scale = TRUE)<jupyter_output><empty_output><jupyter_text>(b) Use the standard deviations to calculate the proportion of variance explained (PVE) by each principal component. Plot the PVE and the cumulative PVE.<jupyter_code>pr.var=pr.out$sdev^2
pve=pr.var/sum(pr.var)

par(mfrow=c(1,2))
plot(pve)
plot(cumsum(pve))<jupyter_output><empty_output><jupyter_text>(c) Use the `pcr()` function in the `pls` library to run PCR on the training set us- ing cross-validation; pass along the parameter `ncomp = 100` and call your object `pcr.fit`. Plot the cross-validation RMSE via the function `validationplot()`.<jupyter_code>pcr.fit=pcr(train.y~., data=train.X, scale=TRUE, validation="CV", ncomp=100)
validationplot(pcr.fit)<jupyter_output><empty_output><jupyter_text>(d) What is the number of components that minimize the RMSE? You can extract the errors for each component via `RMSEP(pcr.fit)$val[2,1,1:101]`.<jupyter_code>which.min(RMSEP(pcr.fit)$val[2,1,1:101])<jupyter_output><empty_output><jupyter_text>(e) Compute the test set RMSE using the number of components you found in the previous item. How does this compare to our previous results? Recall that the test RMSE for linear regression was 0.7527, for the Lasso it was 0.5209, for Boosted Trees it was 0.5115 and for Random Forests it was 0.5123.<jupyter_code>pcr.pred=predict(pcr.fit, test.X, ncomp=89)
pcr.RMSE=sqrt(mean((pcr.pred-test.y)^2))
pcr.RMSE<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1230340760157274, 'fraction_numerical': 0.03915465268676278, 'mean_word_length': 3.992639040348964, 'pattern_counts': {'":': 0, '<': 108, '<?xml version=': 0, '>': 108, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '26', 'dir': 'jupyter-structured-clean-dedup', 'id': '515161', 'license': 'MIT', 'n_tokens_mistral': 6755, 'n_tokens_neox': 6180, 'n_words': 1848, 'path': 'notebooks/stats216v/hw-4-gyuho.ipynb', 'repo_name': 'gyuho/dplearn', 'size': '18329'}