8.判别分析

#####8 判别分析#####
#判别分析就是根据已掌握的每个类别若干样本的数据信息，总结出客观事物分类的规律性，
#建立判别公式和判别标准；在遇到新的样本点时，再根据已总结的判别公式和判别准则，
#来判断出样本点所属的类别。

#####8.1 概述#####

#费希尔(Fisher)判别:线性判别分析方法(LDA)和二次判别方法(QDA);

#贝叶斯(Bayes)判别:朴素贝叶斯分类(Naive Bayesian Classification)算法;

#距离判别:K最近邻(kNN)及有权重的K最近邻算法

#####8.1.1 费希尔判别#####
#费希尔判别的基本思想就是“投影”，即将高维空间的点向地位空间投影，从而简化问题进
#行处理。

#费希尔判别最重要的就是选择合适的投影轴，对该投影轴方向上的要求是：保证投影后，
#使每一类之内的投影值所形成的类内离差尽可能的小，而不同类之间的投影值所形成的类
#间离差尽可能大，即在该空间内有最佳的可分离性，以此获得较高的判别结果。

#二次判别与线性判别的区别就在于投影面的形状不同，二次判别使用若干二次曲面，而非
#直线或平面来将样本划分至相应的类别中。

#####8.1.2 贝叶斯判别#####
#根据已知的先验概率，利用贝叶斯公式，求出后验概率。即该样本属于某一类的概率，然
#后选择最大的后验概率的类作为该样本所属的类。

#####8.1.3 距离判别#####
#根据待判定样本与已知类别样本之间的距离远近作出的判别。

#K最近邻：如果一个样本在特征空间中的K个最相似/最近邻的样本中的大多数属于一个类别，
#则该样本也属于这个类别。

#有权重的K最近邻算法，即对各已知类别样本点根据其距离未知样本点的远近，赋予了不同
#的权重，即距离越近的权重越大。

######8.2 R中的实现#####
#线性判别分析(LDA)      MASS    lda()
#二次判别分析(QDA)      MASS    qda()
#朴素贝叶斯分类         klaR    NaiveBayes()
#K最近邻(KNN)           class   knn()
#有权重的K最近邻        kknn    kknn()

#####8.2.2 核心函数#####
#1. lda()函数
#有三种使用格式，在默认情况下，即对象为数据框data.frame时，其基本格式为
#lda(x, grouping, prior=proportions, tol=1.0e-4, method, CV=FALSE, nu, ...)
#另有两种适合公式formula形式和矩阵matrix形式的两种格式：
#lda(formula, data, ..., subset, na.action)
#lda(x, grouping, ..., subset, na.action)

#其中，x为该函数要处理的数据框data.frame或数据矩阵matrix；formula则放置用于生产判
#别规则的公式，以y~x1+x2+x3格式呈现；data和subset都用于以formula为对象的函数格式
#中，分别用于指明该formula中遍历所来自的数据集名称及所纳入规则建立过程的样本；
#grouping则指明每个观测样本所属类别；prior可设置各类别的先验概率，在无设置的情况
#下，R默认取训练集中各类别样本的比例；tol用于保证判别效果，可通过设置筛选变量，默
#认取0.0001；na.action用于选择对缺失值的处理，默认情况下，若有缺失值，则该函数无
#法运行，当更改设置为na.omit时，则自动删除在用于判别的特征变量中含有缺失值的观测
#样本。

#2. qda()函数
#该函数与lda()一样，也有三种用于数据框、公式和矩阵对象的函数格式，默认(数据框为
#对象)格式为
#qda(x, grouping, prior=proportions, method, CV=FALSE, nu, ...)
#适用于公式和及矩阵
#qda(formula, data, ..., subset, na.action)
#qda(x, grouping, ..., subset, na.action)

#3. NaiveBayes()函数
#两种格式，一种为默认
#NaiveBayes(x, grouping, prior, usekernel=FALSE, fL=0, ...)
#当对象为公式时，则取
#NaiveBayes(formula, data, .., subset, na.action=na.pass)

#其中x\grouping等参数和之前的是一样的。需要注意的是，虽然该函数中也有na.action参
#数，但与lda()和qda()中的不同，此处的默认为np.pass，表示不将缺失值纳入计算，并不
#会导致函数无法运行，当取值为na.omit时则与lda()函数相同，表示删除相应的含有缺失
#值的观测样本。
#另外，usekernel参数用于选择函数计算过程中，密度估计所采用的算法，默认时取FALSE,
#表示使用标准密度估计，也可以通过取值TRUE，选择使用核密度估计法。
#fL用于设置进行拉普拉斯修正的参数值，默认取0，即不进行修正，该修正在数据较小的情
#况下十分必要。

#4. knn()函数
#knn(train, test, cl, k=1, 1=0, prob=FALSE, ues.all=TRUE)

#knn()函数默认选择欧氏距离来寻找所需的K的最近样本。在可变参数中，train和test参数
#分别代表训练集和测试集；cl用于放置训练集中各已知类别样本的类别取值；k为控制最近
#邻域大小的参数；prob取TRUE时，可输出该待判样本所对应的prob值。use.all用于选择再
#出现“结点”时的处理方式，所为结点即指距离待判样本第K近的样本不止一个，默认取TRUE
#时，会将结点都纳入判别过程，这时，n个K近邻就有K+1个样本，若取FALSE，则会随机选
#一个以保证K近邻中刚好有K个样本。

#5. kknn()函数
#kknn(formula=formula(train), train, test, na.action=na.omit(), k=7, distance=2,
#kernel="optimal", ykernal=NULL, scale=NULL, contrasts=c('unordered'="contr.dummy",
#ordered="contr.ordinal"))

#其中主要参数在之前都已经说明，此处仅解说distance参数。该参数用于设定选择计算样本
#间距离的具体方法，通过设定明氏距离(Minkowski Distance)中的参数来实现，取1或2是最
#常用的，参数取2即为欧氏距离，取1为曼哈顿距离，当取无穷时的极限情况下，可以得到切
#比雪夫距离。

#####8.2.3 数据集#####
#本章选用了kknn软件包中的miete数据集进行算法演示，该数据集记录了1994年慕尼黑的住房
#租金标准中的一些有趣变量，比如房子的面积、是否有浴室、是否有中央供暖、是否供应热
#水等，这些都影响并决定着租金。

#1. 数据概况
library(kknn)
data(miete)
head(miete)
dim(miete)  #显示数据集的维度
summary(miete) #获取数据集中各变量的基本信息

#nm,     定量，净租金(单位：德国马克)，      否(本章不使用)，原因为使用相应的定性变量nmkat
#wfl,    定量，占地面积(单位：平方米),       是
#bj,     定量，建造年份，                    否，            使用相应的定性变量djlat
#bad0,   定性，是否有浴室(1无，0有)，        是，
#zh,     定性，是否有中央空调(1有,0无)，     是，
#ww0,    定性，是否提供热水(1无，0有)，      是
#badkach,定性，是否有铺瓷砖的浴室(1有，0无)，是
#fenster,定性，窗户类型(1普通，0优质)，      是
#kueche, 定性，厨房类型(1设施齐全，0普通)，  是
#mvdauer,定量，可租凭期(单位：年) 是
#bjkat， 定性，按区间划分的建造年份(1:1919年前；2:191-1948年；3:1949-1965年；4:1966-1977年；5:1978-1983年；6:1983年后) ，是
#wfkat,  定性，按区间划分的占地面积wfl(1:少于50平方米；2:51-80平方米；3:至少80平方米),否，使用相应的定量变量wfl
#nmqm,   定量，净租金/每平方米，             是
#rooms,  定性，房间数，                      是
#nmkat,  定性，按区间划分的净租金nm(1:<500马克;2:500-675马克；3:675-850马克；4:850-1150马克；5:至少1150马克)，是
#adr,    定性，地理位置(1:差; 2:中; 3优)，   是
#wohn,   定性，住宅环境(1:差; 2:中; 3优)，   是

#我们在数据集中剔除含义重复的第1、3、12这三个变量，取余下的14个变量进行处理。且其
#中我们选择第15各变量——安区间划分的净租金(nmkat)作为待判别变量。

#2. 数据预处理
#划分训练集和数据集，为提高判别效果，考虑用分层抽样。
library(sampling)         #加载用于获取分层抽样函数strata()函数的软件包sampling
n=round(2/3*nrow(miete)/5)#按照训练集占数据总量的2/3的比例，计算每一等级中应抽取的样本量
n
sub_train = strata(miete,stratanames="nmkat",size=rep(n,5),method="srswor")
#以nmkat变量的5各等级划分层次，进行分层抽样
head(sub_train) #显示训练集抽取情况，包括nmkat变量取值，该样本在数据集中的序号，
#被抽取到的概率，以及所在的层次。

data_train=miete[sub_train$ID_unit,c(-1,-3,-12)]
#获取如上ID_unit所对应样本构成训练集，并删除变量1、3、12
data_test=miete[-sub_train$ID_unit,c(-1,-3,-12)]
#获取如上ID_unit所对应样本之外的数据构成测试集，并删除变量1、3、12
dim(data_train);dim(data_test) #选择训练集、测试集维度，检查抽样结果
head(data_test) 

#####8.3.1 线性判别分析#####
library(MASS)

#1.公式formula格式
fit_ldal=lda(nmkat~.,data_train) #以公式格式执行线性判别
names(fit_ldal)                  #查看lda()可给输出项名称，
#有10个输出结果，分别为执行过程中所使用的先验概率prior、数据集中各类别的
#样本量counts、各变量在每一类别中的均值means等。

fit_ldal$prior            #查看本次执行过程中所使用的先验概率
fit_ldal$counts           #查看数据集data_train中各类别的样本量
#因为之前抽样采用的是nmkat各等级的等概率分层抽样方式，因此如上各类别的先
#验概率和样本量都是相等的。
fit_ldal$means            #查看各变量在每一类别中的值
fit_ldal                  #输出判别分析的各项结果

#2. 数据框data.frame及矩阵matrix格式
fit_lda2=lda(data_train[,-12],data_train[,12])
#分别设置属性变量(除第12个变量nmkat外)与待判别变量(第121各变量nmkat)的取
#值，并记为fit_lda2
fit_lda2

#3. 判别规则可视化
plot(fit_ldal)
#我们可以看到，在所有4个线性判别式(LD)下1至5这5个类别的分布情况，不同类
#别样本已用相应数字标出。
#另外，我们可以通过dimen参数的设定来控制输出图形中所使用的判别式个数。
#当dimen参数取值大于总共的判别式个数时，则默认取所有判别式。
plot(fit_ldal,dimen=1) #对判别规则fit_ldal输出1个判别式的图形,对于图边距过大，拉大图形框即可
plot(fit_ldal,dimen=2) #对判别规则fit_ldal输出2个判别式的图形

#4. 对测试集待判变量取值进行预测

pre_ldal=predict(fit_ldal,data_test) #使用判别规则fit_ldal预测data_test中nmkat变量的类型
pre_ldal$class     #输出data_test中各样本的预测结果
pre_ldal$posterior #输出各样本属于每一类别的后验概率
#后验概率最高者为该样本被判定的类别。
table(data_test$nmkat,pre_ldal$class) #生成nmkat变量的预测值与实际值的混淆矩阵

#行表示实际的类别，列表示预测判定的类别。在362个测试样本中，实际属于第一类的有
#75个，而由判定结果，75个样本中，有55个判定正确，17个被错判2类，1个错判3类，1个
#错判4类，1个错判5类。非对角线上的元素都是被判错的。错误率很容易就计算出来了。
error_ldal=sum(as.numeric(as.numeric(pre_ldal$class)!=as.numeric(
  data_test$nmkat)))/nrow(data_test) #计算错误率
error_ldal

#lda2作不了图，也无法带入预测函数中，懵逼

#####8.3.2 朴素贝叶斯分类#####
library(klaR)

#1. 公式formula格式
fit_Bayes1=NaiveBayes(nmkat~.,data_train)#以nmkat为待判别变量，以data_train来生成贝叶斯判别规则
names(fit_Bayes1)                        #显示fit_Bayes1所包含的输出项名称
#apriori为先验概率
fit_Bayes1$apriori
#tables项存储了用于建立判别规则的所有变量在各类别下的条件概率
fit_Bayes1$tables
#剩余的判别变量等级项levels、判别命令项call、是否使用标准密度估计usekernel，以及
#参与判别规则制定的特征变量名varnames。
fit_Bayes1$levels
fit_Bayes1$call
fit_Bayes1$usekernel
fit_Bayes1$varnames

#2. 各类别下变量密度可视化
plot(fit_Bayes1,vars="wfl",n=50,col=c(1,"darkgrey",1,"darkgrey",1))
                                               #对占地面积whl绘制个类别下的密度图,图例大小不变，拉开窗口即可
plot(fit_Bayes1,vars="mvdauer",n=50,col=c(1,"darkgrey",1,"darkgrey",1)) 
                                               #绘制个类别下的密度图
plot(fit_Bayes1,vars="nmqm",n=50,col=c(1,"darkgrey",1,"darkgrey",1)) 
                                               #绘制个类别下的密度图

#3. 默认格式
fit_Bayes2=NaiveBayes(data_train[,-12],data_train[,12])

#4. 对测试集待判别变量取值进行预测
pre_Bayes1=predict(fit_Bayes1,data_test) #根据fit_Bayes1判别规则对测试集进行预测
pre_Bayes1
table(data_test$nmkat,pre_Bayes1$class) #混淆矩阵
error_Bayes1=sum(as.numeric((as.numeric(pre_Bayes1$class)!=
                               as.numeric(data_test$nmkat)))/nrow(data_test))
error_Bayes1
#判别效果不佳，基本等同于纯猜测得到的正确程度。这很可能是由于不符合朴素贝叶斯判
#别执行的前提条件——各变量相互独立，这些变量都是显著相关的，这在很大程度上影响了
#预测结果的好坏。

######8.3.3 K最近邻#####
library(class)
#按次序向knn()函数依次放入各属性变量(除去待判变量)、测试集(除去待判变量)、训练集，
#并首先取K的默认值1进行判别。
fit_pre_knn=knn(data_train[,-12],data_test[,-12],data_train[,12])
                         #建立K最近邻判别规则，并对测试样本集进行预测
fit_pre_knn
table(data_test$nmkat,fit_pre_knn) #生成nmkat真实值与预测值的混淆矩阵
error_knn=sum(as.numeric((as.numeric(fit_pre_knn)!=as.numeric(data_test$nmkat)))/
                nrow(data_test))   #计算错误率
error_knn

#我们看到K取1时，错误路仅为0.34，判别效果较前几种算法都要好。这与数据集的特点密
#不可分，在其他数据中也可能是另一种算法表现更优，在实际中需要针对不同的数据集选
#取不同的挖掘算法。

#下面我们通过调整K的取值，选择出最适合该数据集的K值，将寻找范围控制在1至20.
error_knn=rep(0,20)#对将用于存储K取1至20时预测错误率的error_knn变量赋初始值为0
for(i in 1:20){
  fit_pre_knn=knn(data_train[,-12],data_test[,-12],cl=data_train[,12],k=i)
                     #构造K最近邻判别规则并预测，预测结果存于fit_pre_knn
  error_knn[i]=sum(as.numeric((as.numeric(fit_pre_knn)!=as.numeric
                               (data_test$nmkat)))/nrow(data_test))
                     #计算取每一个K值时所得到的错误率error_knn
}
error_knn #显示错误率向量的20个取值
plot(error_knn,type="l",xlab="K") #对20个错误率值作折线图
#从图中可以明显的看出当K取几时预测效果最佳
#但K最近邻算法也有其缺陷，当样本不平衡，某些类的养不了很大，而其他样本量很小时，
#容易导致当输出一个新样本时，该样本的K个最近邻中大样本类别的样本占多数。

#####8.3.4 有权重的K最近邻算法#####
library(kknn)
#将规则建立公式nmkat~.,以及训练集和测试集分别放入函数，设置K为5.
fit_pre_kknn=kknn(nmkat~.,data_train,data_test[,-12],k=5)
                              #建立有权重的K最近邻判别规则，并对测试集样本进行预测
summary(fit_pre_kknn)                             #输出在K近邻判别规则下的判别结果
fit=fitted(fit_pre_kknn)                          #输出有权重的K最近邻预测结果
fit
table(data_test$nmkat,fit)                        #混淆矩阵
error_knn=sum(as.numeric((as.numeric(fit)!=as.numeric(data_test$nmkat)))/
                nrow(data_test))                  #计算错误率
error_knn

#####8.4 推荐系统综合实例#####
#基于用户的KNN,找出K个与A对其他电影给予相似评分，且对电影M已经进行评分的用户，然后
#再用这K个用户对M的评分来预测A对M的评分。(User-based KNN)

#基于项目的KNN,找出K个与电影M相似的，并且A评价过的电影，然后再用这K部电影的评分来
#预测A对M的评分。(Item-based KNN)

#####8.4.2 MovieLens数据集说明#####
#http://ww.grouplens.org/node/73

#u.data:含有943位用户对1682部电影总计10万条评分，且每位用户至少记录了其对20部电影
#的评分。用户ID、电影ID、评分、以及时间戳，排列时无序的，主要用前三个变量。

#u.item:记录每部电影的信息，包括电影ID、电影名称、上映时间、视频发布时间、网络电影
#资料库的网址，以及是否为某类型的电影的一系列二分变量。这是探究各电影相似性的重要
#数据资料。

#u.user:记录每位用户的基本信息，包括用户ID、年龄、性别、职业以及右边。这是探究个用
#户间相似性的重要信息来源。

#####8.4.3 综合应用#####
#本次应用User-based KNN算法探究MovieLens数据集
#目的：预测某位用户对某部电影的评分。

setwd("E://books/数据挖掘 R语言实战/ml-100k") #设置常用工作路径
data=read.table("u.data")                     #读取数据文件
data=data[,-4]                                #删除其中不需要使用的第4列
names(data)=c("userid","itemid","rating")     #命名data的各变量名
head(data);dim(data)                          #显示数据集data的前若干项，及其维度

#MovieLens_KNN函数
#为了程序可反复使用，我们来手动编写一个专用于实现“预测某位用户U对某部电影M评分值”
#的函数——MovieLens_KNN()
Userid=1;Itemid=61;n=50;K=10                  #设定参数取值
MovieLens_KNN=function(Userid,Itemid,n,K){
  sub=which(data$userid==Userid)#获取带预测用户在数据集中各条信息所在的行标签，存于sub
  if(length(sub)>=n) sub_n=sample(sub,n)
  if(length(sub)<n) sub_n=sample(sub,length(sub))
                                #获取随机抽取的n个预测用户已评分电影的行标签
  known_itemid=data$itemid[sub_n]#获取预测用户已评分电影的ID
  unknown_itemid=Itemid          #获取带预测电影的ID号
  unknown_sub=which(data$itemid==unknown_itemid)
  user=data$userid[unknown_sub[-1]] #获取已评价电影的用户ID
  data_all=matrix(0,1+length(user),2+length(known_itemid))#设置data.all的行数，列数，所有值暂时取0
  data_all=data.frame(data_all)
  names(data_all)=c("userid",paste("unknown_itemid_",Itemid),
                    paste("itemid_",known_itemid,sep=""))#对data.all的各变量命名
  item=c(unknown_itemid,known_itemid)
  data_all$userid=c(Userid,user) #对data.all中的useris变量赋值
  for (i in 1:nrow(data_all)){   #对data_all按行进行外层循环
    data_temp=data[which(data$userid==data_all$userid[i]),]
    for (j in 1:length(item)){   #对data_all按行进行内层进行循环
      if(sum(as.numeric(data_temp$itemid==item[j]))!=0) #判断该位置是否有取值
      {data_all[i,j+1]=data_temp$rating[which(data_temp$itemid==item[j])]}
    }
  }
  data_test_x=data_all[1,c(-1,-2)] #获取测试集的已知部分
  data_test_y=data_all[1,2]        #获取测试集的待预测值
  data_train_x=data_all[-1,c(-1,-2)]#获取训练集的已知部分
  data_train_y=data_all[-1,2]         #获取训练集的待预测部分
  library(class)
  fit=knn(data_train_x,data_test_x,cl=data_train_y,k=K) #进行knn判别
  list("data_all:"=data_all, "True Rating:"=data_test_y,"Predict Rating:"=fit,
       "UserID:"=Userid,"ItemID:"=Itemid) #设置各项输出结果
}
#下面，我们简单利用MovieLens_KNN()函数来看一看用户1对电影1至电影20这20部电影的
#评分情况。
user1=NULL
for(Item in 1:20)user1=c(user1,MovieLens_KNN(Userid = 1,Itemid = Item,n=50,K=10)$'True Rating:')
user1
which(user1==5) #显示评分为5的电影ID

user2=NULL
for(Item in 1:20)user2=c(user2,MovieLens_KNN(Userid = 1,Itemid = Item,n=50,K=10)$'Predict Rating:')
user2
which(user2==5) #显示评分为5的电影ID
table(user1,user2)