---
title: "第11章: 支持向量机和神经网络"
output:
  html_document:
    df_print: default
---

## 统一设置ggplot2的绘图风格
```{r}
library(ggplot2)
theme_set(theme_bw(base_family = "STKaiti"))

```


    
## 11.1：支持向量机分类

### Breast Cancer Wisconsin (Diagnostic) Data Set
### Predict whether the cancer is benign or malignant

https://www.kaggle.com/uciml/breast-cancer-wisconsin-data

```{r}
library(e1071)
library(Metrics)
library(readr)
library(gmodels)
## 读取数据
Cancer <- read_csv("data/chap11/Breast Cancer Wisconsin.csv")
Cancer$diagnosis <- as.factor(Cancer$diagnosis)
head(Cancer)
summary(Cancer)

table(Cancer$diagnosis)

```
### 主成分分析对数据进行降维

```{r}
library(psych)
library(ggplot2)
library(pheatmap)
## 可视化碎石图，选择合适的主成分数
parpca <- fa.parallel(Cancer[,3:32],fa = "pc")

## 主成分分析,提取前2个主成分
can_cor <- cor(Cancer[,3:32],Cancer[,3:32])
dim(can_cor)
canpca <- principal(can_cor,nfactors = 5,rotate = "cluster")
## 可视化主成分和每个变量的相关性
pheatmap(canpca$loadings,cluster_rows = F,cluster_cols = F)
## 使用pca模型获取数据集的5个主成分
can_score<- predict.psych(canpca,Cancer[,3:32])
dim(can_score)
head(can_score)
## 可视化前两个维度的数据分布
plotdata <- data.frame(PC1 = can_score[,1],PC2 = can_score[,2],
                       diagnosis = Cancer$diagnosis)

ggplot(plotdata,aes(x=PC1,y= PC2,colour = diagnosis,shape = diagnosis))+
  geom_point()+theme(legend.position = "top")


```

### 支持向量机分类

```{r}
## 数据70%用于训练，30%用于测试模型效果
## 为了方便可视化，只使用前两个主成分
set.seed(123)
index <- sample(nrow(can_score),round(0.7*nrow(can_score)))
cancerdata <- data.frame(can_score[,1:2])
cancerdata$diagnosis <- as.factor(Cancer$diagnosis)
colnames(cancerdata) <- c("PC1","PC2","y")
train_data <- cancerdata[index,]
test_data <- cancerdata[-index,]


## 线性核的支持向量机分类模型
can_svm <- svm(y~.,data = train_data,kernel ="linear")
summary(can_svm)
## 在训练数据集上，可视化支持向量和分类面
plot(can_svm,train_data)
## 查看模型在测试集上的预测效果
test_pre <- as.character(predict(can_svm,test_data))
sprintf("base SVM model acc: %f",accuracy(test_data$y,test_pre))
table(test_data$y,test_pre)
CrossTable(test_data$y,test_pre,prop.r=F, prop.c=F,prop.t=F, prop.chisq=F)

```


### 参数搜索

```{r}
## 模型的精度能够提高吗?
## 调整参数提升模型的精度
set.seed(1)
svm_opt <- tune.svm(y~.,data = train_data,kernel ="linear",cost = seq(1,20,1))

## 可视化参数搜索的结果
plot(svm_opt,main = "SVM best cost parameter")
svm_opt$best.parameters
svm_opt$best.performance
svmbest <- svm_opt$best.model
summary(svmbest)

## 查看模型在测试集上的预测效果
test_pre <- as.character(predict(svmbest,test_data))
sprintf("best SVM model acc: %f",accuracy(test_data$y,test_pre))
table(test_data$y,test_pre)
CrossTable(test_data$y,test_pre,prop.r=F, prop.c=F,prop.t=F, prop.chisq=F)

## 针对线性核，虽然不同的cost参数下支持向量的数量不一样，但是在测试集上的预测精度是一样的
```


### 更改支持向量机模型的分类kernal，查看模型的效果

```{r}

## 进一步的提高模型精度，搜索更多的参数
set.seed(1)
Optvm <- tune.svm(y~.,data = train_data,kernel ="radial",
                    cost = seq(1,20,1),gamma = seq(0.1,1,0.1))
plot(Optvm)

Optvm$best.parameters
Optvm$best.performance
svmbest2 <- Optvm$best.model
summary(svmbest2)
## 可视化新的支持向量机分类面
plot(svmbest2,data = train_data)
## 查看模型在测试集上的预测效果
test_pre <- as.character(predict(svmbest2,test_data))
sprintf("best SVM model acc: %f",accuracy(test_data$y,test_pre))
table(test_data$y,test_pre)
CrossTable(test_data$y,test_pre,prop.r=F, prop.c=F,prop.t=F, prop.chisq=F)

## 使用新的核函数得到了新的数据切分曲面，模型的分类效果比使用线性核效果好

```



## 11.2： SVM识别异常值

  识别二维数据中的异常值

  异常值识别垃圾邮件

```{r}
## 二维数据中的异常值
library(ggplot2)
library(DMwR)
library(e1071)
library(Rlof)


# ## 生成随机数据
# set.seed(1234)
# samp1 <- mvrnorm(n = 200,mu=c(5,5),
#                  Sigma = matrix(c(1,0,0,1),nrow = 2,ncol = 2))
# samp1 <- samp1[samp1[,1] < 7,]
# samp2 <- mvrnorm(n = 3,mu=c(0,0),
#                  Sigma = matrix(c(5,0,0,5),nrow = 2,ncol = 2))
# samp3 <- mvrnorm(n = 2,mu=c(6,12),
#                  Sigma = matrix(c(1,0,0,2),nrow = 2,ncol = 2))
# samp4 <- mvrnorm(n = 2,mu=c(12,3),
#                  Sigma = matrix(c(3,0,0,4),nrow = 2,ncol = 2))
# samp5 <- mvrnorm(n = 1,mu=c(12,6),
#                  Sigma = matrix(c(1,0,0,1.5),nrow = 2,ncol = 2))
# samp6 <- mvrnorm(n = 2,mu=c(11,10),
#                  Sigma = matrix(c(2,0,0,2),nrow = 2,ncol = 2))
# samp7 <- mvrnorm(n = 3,mu=c(-1,10),
#                  Sigma = matrix(c(6,0,0,6),nrow = 2,ncol = 2))
# samp <- rbind(samp1,samp2,samp3,samp4,samp5,samp6,samp7)
# par(pty = c("s"))
# plot(samp[,1],samp[,2])
# ## 随机分布样本
# scatter2dim <- data.frame(x = samp[,1],y = samp[,2])
# index <- sample(nrow(scatter2dim),nrow(scatter2dim))
# scatter2dim <- scatter2dim[index,]
# plot(scatter2dim$x,scatter2dim$y)
# write.csv(scatter2dim,"data/chap10/outlierdata.csv",row.names = F)

## 读取异常值数据
outlierdata <- read.csv("data/chap11/outlierdata.csv")
## 可视化数据的分布
qplot(outlierdata$x,outlierdata$y,xlab = "X",ylab = "Y",main = "散点图")

## Lof方法识别异常值
lof_score <- lofactor(outlierdata,k=6)
outlierdata$lof_score <- lof_score
outlierdata$scorecut <- cut_width(lof_score,2)
## 可视化LOF得分
ggplot(outlierdata,aes(x=x,y=y,colour=lof_score))+
  geom_point(aes(shape = scorecut))+
  ggtitle("LOF得分")

## 可以指定距离，多核并行计算lof得分
lof_score2 <- Rlof::lof(outlierdata[,c("x","y")],k=6,cores = 2,
                        method = "maximum")
outlierdata$lof_score2 <- lof_score2
outlierdata$scorecut2 <- cut_width(lof_score2,2)
## 可视化LOF得分
ggplot(outlierdata,aes(x=x,y=y,colour=lof_score2))+
  geom_point(aes(shape = scorecut2))+
  ggtitle("LOF得分")

## 使用支持向量机来识别异常值SVM
out_svm <- svm(x = outlierdata[,c("x","y")],y = NULL,
               type = "one-classification",nu = 0.08)
out_svm
## 预测是否为异常值
svmpre <- predict(out_svm,outlierdata[,c("x","y")])
outlierdata$isoutlier <- !svmpre
## 可视化SVM的识别结果
ggplot(outlierdata,aes(x=x,y=y))+
  geom_point(aes(shape = isoutlier))+
  ggtitle("SVM识别异常值")




```

## 使用SVM异常值处理方法识别垃圾邮件

```{r}
library(tm)
library(gmodels)
load("data/chap8/spam_cp.RData")
load("data/chap8/spam.RData")
## 构建TF－IDF矩阵
controls = list(weighting = function(x) weightTfIdf(x,normalize = FALSE))
spam_dtm <- DocumentTermMatrix(spam_cp,control = controls)
## 通常情况下，文档－词项的tf-idf矩阵非常的稀疏，可以删除一些不重要的词来缓解矩阵的稀疏性，同时提高计算效率
spam_dtm <- removeSparseTerms(spam_dtm,0.999)
spam_dtm
## 数据随机切分为75%训练集和25%测试集
set.seed(123)
index <- sample(nrow(spam),nrow(spam)*0.75)
spam_dtm2mat <- as.matrix(spam_dtm)
train_x <- spam_dtm2mat[index,]
train_y <- spam$label[index]
test_x <- spam_dtm2mat[-index,]
test_y <- spam$label[-index]

# train_x <- apply(train_x, 2, function(x) ifelse(x>0,1,0))
# test_x <- apply(test_x, 2, function(x) ifelse(x>0,1,0))
## 使用支持向量机来识别异常值SVM
out_svm <- svm(x = train_x,y = NULL,
               type = "one-classification",nu = 0.1)
## 预测是否为异常值
svmpre <- predict(out_svm,test_x)
table(svmpre)
table(test_y)
svmpre2 <- ifelse(svmpre == TRUE, "ham","spam")

CrossTable(test_y,svmpre2,prop.r = F,prop.t = F,prop.chisq = F)
table(test_y,svmpre2)
sprintf("异常值检测识别精度为%.4f",accuracy(test_y,svmpre2))



```


## 11.3：支持向量回归

### House Sales in King County, USA

https://www.kaggle.com/harlfoxem/housesalesprediction


### 读取数据，数据特征处理，数据切分
```{r}
library(e1071)
library(caret)
library(Metrics)
library(readr)
library(gmodels)
library(outliers)
library(GGally)
library(corrplot)
library(lubridate)
kc_house <- read_csv("data/chap11/kc_house_data.csv")
kc_house$date <- year(kc_house$date)
kc_house <- kc_house[,2:21]
## 剔除价格异常的数据
outlierval <- boxplot.stats(kc_house$price)$out
kc_house <- kc_house[!(kc_house$price %in% outlierval),]
dim(kc_house)
## 剔除 bedrooms 异常的数据(==0,>6)
index <- c(which(kc_house$bedrooms == 0),which(kc_house$bedrooms > 6))
kc_house <- kc_house[-index,]
dim(kc_house)
## 剔除 bathrooms 取值为异常值的数据
outlierval <- boxplot.stats(kc_house$bathrooms)$out
kc_house <- kc_house[!(kc_house$bathrooms %in% outlierval),]
dim(kc_house)
## 剔除 sqft_lot 取值为异常值的数据
outlierval <- boxplot.stats(kc_house$sqft_lot)$out
kc_house <- kc_house[!(kc_house$sqft_lot %in% outlierval),]
dim(kc_house)
## 剔除特征 waterfront and view
kc_house$waterfront <- NULL
kc_house$view <- NULL
## 剔除 condition 取值为异常值的数据
outlierval <- boxplot.stats(kc_house$condition)$out
kc_house <- kc_house[!(kc_house$condition %in% outlierval),]
dim(kc_house)
## 剔除 sqft_basement 取值为异常值的数据
outlierval <- boxplot.stats(kc_house$sqft_basement)$out
kc_house <- kc_house[!(kc_house$sqft_basement %in% outlierval),]
dim(kc_house)
## 房子建立的时间长短,并剔除取值小于0的数据
kc_house$yr_built <- kc_house$date - kc_house$yr_built
index <- which(kc_house$yr_built < 0)
kc_house <- kc_house[-index,]
dim(kc_house)
## 房子是否装修过,装修过为1，没装修过为0
index <- which(kc_house$yr_renovated != 0)
kc_house$yr_renovated[index] <- 1
dim(kc_house)
kc_house$zipcode <- as.factor(kc_house$zipcode)
## 剔除不需要的特征
kc_house <- kc_house[,2:14]
## 可视化房子的售价分布
ggplot(kc_house)+geom_histogram(aes(price),fill = "lightblue",bins = 60)

## 可视化特征之间的相关系数
kc_cor <- cor(kc_house[1:11],kc_house[1:11])
corrplot.mixed(kc_cor,tl.col="black",tl.pos = "lt",
               tl.cex = 0.8,number.cex = 0.7)

head(kc_house)
summary(kc_house)
dim(kc_house)
anyNA(kc_house)

write.csv(kc_house,"data/chap11/kc_house_clean.csv",row.names = FALSE)
## 切分数据
set.seed(123)
index <- sample(nrow(kc_house),round(0.7*nrow(kc_house)))
train_data <- kc_house[index,]
test_data <- kc_house[-index,]
dim(train_data)
dim(test_data)
```

### 支持向量机回归模型

```{r}

## 支持向量机回归模型
system.time(
  svmreg <- svm(price~.,data = train_data,kernel = "radial")
)

summary(svmreg)
## 在训练集上的误差
train_pre <- predict(svmreg,train_data)
train_mae <- mae(train_data$price,train_pre)
sprintf("训练集上的绝对值误差: %f",train_mae)
test_pre <- predict(svmreg,test_data)
test_mae <- mae(test_data$price,test_pre)
sprintf("测试集上的绝对值误差: %f",test_mae)


```


### 使用参数搜索，寻找更好的模型

```{r}

## 该部分程序运行时间较长，此处先将其注释掉 
#############################################
# set.seed(234)
# tunecontrols <- tune.control(sampling = "bootstrap",boot.size = 0.8)
# system.time(

#   svmregopt <- tune(svm,price~.,data = train_data,kernel = "radial",
#                     ranges = list(gamma = seq(0.01,0.1,0.02),
#                                   cost= seq(11,20,2)),
#                     tunecontrol = tunecontrols)
# )
## 找到的模型参数
# plot(svmregopt,main = "Performance of svm regression")
# svmregopt$best.parameters
# svmregopt$best.performance
#############################################

## 得到的最优参数为cost=19,gamma=0.01

 
 
## 使用找到的最好参数进行SVM回归
system.time(
  svmregbest <- svm(price~.,data = train_data,kernel = "radial",
                    cost=19,gamma=0.01)
)

## 在训练集上的误差
train_pre <- predict(svmregbest,train_data)
train_mae <- mae(train_data$price,train_pre)
sprintf("训练集上的绝对值误差: %f",train_mae)
test_pre <- predict(svmregbest,test_data)
test_mae <- mae(test_data$price,test_pre)
sprintf("测试集上的绝对值误差: %f",test_mae)

## 模型和原始的支持向量回归相比，效果提升了5000多
```


## 11.4：MLP神经网络分类

### 单隐藏层神经网络，分类肺癌数据

```{r}
library(neuralnet)
library(NeuralNetTools)
library(dplyr)
library(readr)
library(psych)
library(Metrics)
library(ggplot2)
## 读取数据
Cancer <- read_csv("data/chap11/Breast Cancer Wisconsin.csv")
Cancer$diagnosis <- as.integer(as.factor(Cancer$diagnosis))-1
## 主成分分析,提取前5个主成分
can_cor <- cor(Cancer[,3:32],Cancer[,3:32])
canpca <- principal(can_cor,nfactors = 5,rotate = "cluster")
can_score<- predict.psych(canpca,Cancer[,3:32])
## 数据70%用于训练，30%用于测试模型效果
set.seed(123)
index <- sample(nrow(can_score),round(0.7*nrow(can_score)))
train_data <- as.data.frame(can_score[index,])
train_data$y <- Cancer$diagnosis[index]
test_data <- as.data.frame(can_score[-index,])
test_data$y <- Cancer$diagnosis[-index]
head(train_data)
colnames(train_data)
## 
## 单隐藏层 MLP model
sigmlp <- neuralnet(y~RC1+RC5+RC2+RC3+RC4,data = train_data,
                    hidden = c(10),linear.output = FALSE,
                    act.fct = "logistic",algorithm = "rprop+")

par(cex = 0.6)
plotnet(sigmlp,pos_col = "red", neg_col = "grey")

summary(sigmlp)
## 测试集上的精度
test_pre <- neuralnet::compute(sigmlp,test_data[,1:5])
test_pre <- as.vector(ifelse(test_pre$net.result >0.5,1,0))
accuracy(test_data$y,test_pre)


## 训练单隐藏层神经网络,分析隐藏层的不同在测试集上的精度
set.seed(123)
size <- seq(2,20,1)
netacc <- size
for (ii in 1:length(size)) {
  signet <- neuralnet(y~RC1+RC5+RC2+RC3+RC4,data = train_data,
                    hidden = c(size[ii]),linear.output = FALSE,
                    act.fct = "logistic",algorithm = "rprop+")
  test_pre <- neuralnet::compute(signet,test_data[,1:5])
  test_pre <- as.vector(ifelse(test_pre$net.result >0.5,1,0)) 
  netacc[ii] <- accuracy(test_data$y,test_pre)
}

data.frame(size = size,netacc = netacc) %>%
  ggplot(aes(x=size,y = netacc))+
  geom_line()+geom_point(colour="red")

```


### 多隐藏层神经网络，fashion-mnist数据分类

### 使用h2o.deeplearning：feed-forward multilayer artificial neural network

https://www.kaggle.com/zalando-research/fashionmnist

```{r}
library(h2o)
library(caret)
library(ggcorrplot)
library(RColorBrewer)
## 启动初始化一个h2o实例; 定义为2核同时计算；
h2o.init(nthreads=4,max_mem_size='8G')
## 因为数据量很大，为了节省计算时间,这里只使用部分数据集，完成训练和测试的工作
fashion <- read_csv("data/chap11/fashion-mnist_train.csv")
fashion_test <- read_csv("data/chap11/fashion-mnist_test.csv")
fashion$label <- as.factor(fashion$label)
fashion[2:785] <- fashion[2:785] / 255.0
fashion_test[2:785] <- fashion_test[2:785] / 255.0
dim(fashion)
dim(fashion_test)
table(fashion$label)
table(fashion_test$label)
## 可视化部分数据
## 可视化部分样本的图像
set.seed(123)
index <- sample(nrow(fashion),50)
par(mfrow = c(5,10),mai=c(0.05,0.05,0.05,0.05))
for(ii in seq_along(index)){
  im <- as.numeric(fashion[index[ii],2:785])
  im <- matrix(im,nrow=28,ncol = 28,byrow = F)
  im <- apply(apply(im, 1, rev),1,rev)
  image(im,col = gray(seq(0, 1, length = 256)),xaxt= "n", yaxt= "n")
}

```



```{r}

## 数据切分为训练集和验证集,80%训练集，20%验证集
fashion_h2o <- as.h2o(fashion)
splits <- h2o.splitFrame(data = fashion_h2o, ratios = 0.8,seed = 1234)
train_data <- splits[[1]]
val_data <- splits[[2]]

dim(train_data)
dim(val_data)

## MLP分类模型
mlpclf <- h2o.deeplearning(x = 2:785, y = 1, train_data, activation = "Tanh",
                           hidden=c(100,100,100,100),epochs = 5,
                           loss = "CrossEntropy",score_each_iteration = TRUE,
                           mini_batch_size = 500,validation_frame = val_data,
                           stopping_rounds=0,seed = 123)


summary(mlpclf)
h2o.mse(mlpclf)
## 可视化模型误差的收敛过程
plot(mlpclf)


## 查看5张图像经过MLP模型后的特征的变化情况
index <- seq(200,nrow(fashion),6000)[1:5]
par(mfrow = c(5,5),mai=c(0.1,0.1,0.1,0.1))
for(ii in seq_along(index)){
  imone <- fashion_h2o[index[ii],2:785]
  imlay1 <- h2o.deepfeatures(mlpclf,imone,layer = 1)
  imlay2 <- h2o.deepfeatures(mlpclf,imone,layer = 2)
  imlay3 <- h2o.deepfeatures(mlpclf,imone,layer = 3)
  imlay4 <- h2o.deepfeatures(mlpclf,imone,layer = 4)
  imone <- matrix(imone,nrow = 28,ncol = 28)
  imone <- apply(apply(imone, 1, rev),1,rev) 
  image(imone,axes = FALSE,main = "image")
  image(matrix(imlay1,nrow = 10,ncol = 10),axes = FALSE,main = "layer1")
  image(matrix(imlay2,nrow = 10,ncol = 10),axes = FALSE,main = "layer2")
  image(matrix(imlay3,nrow = 10,ncol = 10),axes = FALSE,main = "layer3")
  image(matrix(imlay4,nrow = 10,ncol = 10),axes = FALSE,main = "layer4")
}


## 计算模型在测试集上的预测值和性能
test_data <- as.h2o(fashion_test)
mlp_pre <- as.data.frame(h2o.predict(mlpclf, newdata = test_data))
head(mlp_pre)

acc <- accuracy(fashion_test$label,mlp_pre$predict)
sprintf("MLP model acc: %f",acc)

## 在测试集上的混淆矩阵热力图
confum <- caret::confusionMatrix(as.factor(fashion_test$label),mlp_pre$predict)
confumat <- as.data.frame(confum$table)
confumat[,1:2] <- apply(confumat[,1:2],2,as.integer)

ggplot(confumat,aes(x=Reference,y = Prediction))+
  geom_tile(aes(fill = Freq))+
  geom_text(aes(label = Freq))+
  scale_x_continuous(breaks = c(0:9))+
  scale_y_continuous(breaks = unique(confumat$Prediction),
                     trans = "reverse")+
  scale_fill_gradient2(low="darkblue", high="lightgreen", 
                       guide="colorbar")+
  ggtitle("MLP分类器在测试集结果")


```



## 11.5：MLP神经网络回归


```{r}
library(RSNNS)
kc_house <- read_csv("data/chap11/kc_house_clean.csv")
## 将zipcode进行label编码
zipcode <- decodeClassLabels(kc_house$zipcode)
kc_house <- cbind(kc_house[1:12],zipcode)
## 数据max-min归一化到0-1之间
kc_house[,2:11] <- normalizeData(kc_house[,2:11],"0_1")
## 房价归一化到0-1之间，并获取归一化参数
price <- normalizeData(kc_house$price,type = "0_1")
NormParameters <- getNormParameters(price)
head(kc_house)
dim(kc_house)
hist(price)
## 数据切分
set.seed(123)
datasplist <- splitForTrainingAndTest(kc_house[,2:82],price,
                                      ratio = 0.3)

## MLP回归模型
system.time(
mlpreg <- mlp(datasplist$inputsTrain,datasplist$targetsTrain,maxit = 200,
              size = c(100,50,100,50),learnFunc = "Rprop",
              inputsTest=datasplist$inputsTest,
              targetsTest = datasplist$targetsTest,
              metric = "RSME")
)
## 可视化模型的效果
plotIterativeError(mlpreg,main = "MLP Iterative Erro")
plotRegressionError(datasplist$targetsTrain,mlpreg$fitted.values,
                    main="MLP train fit")
plotRegressionError(datasplist$targetsTest,mlpreg$fittedTestValues,
                    main="MLP test fit")

## 在训练集上的误差
mlppre_train <- denormalizeData(mlpreg$fitted.values,NormParameters)
mlp_train <- denormalizeData(datasplist$targetsTrain,NormParameters)
train_mae <- mae(mlp_train,mlppre_train)
sprintf("训练集上的绝对值误差: %f",train_mae)
mlppre_test <- denormalizeData(mlpreg$fittedTestValues,NormParameters)
mlp_test <- denormalizeData(datasplist$targetsTest,NormParameters)
test_mae <- mae(mlp_test,mlppre_test)
sprintf("测试集上的绝对值误差: %f",test_mae)



```

