统计建模与R语言笔记
# 判别分析:用以判别个体属于那个总体的一种统计方法
# 1.方法一:距离判断方法(马氏距离)
# 马氏距离较之欧式距离,考虑了随机变量方差的信息
# 分别讨论两总体协方差相同和不同的情况
# 首先讨论两总体的判别分析
discriminiant.distance<-
? function(trnX1,trnX2,testX=NULL,var.equal=FALSE) {
# ? ? 如果测试样本为空,将训练样本1,2合并形成测试样本
? ? ? if (is.null(testX) == TRUE)? testX<-rbind(trnX1,trnX2);
# ? ? 如果测试样本为向量(vector),将其转换为数组(matrix)并转置(t)
? ? ? if (is.vector(testX) == TRUE) testX<-as.matrix(testX);
# ? ? 转换后,如果不是数组,将其转换为数组(此时已经在上一步已经转置)
? ? ? if(is.matrix(testX) != TRUE) testX<-as.matrix(testX);
# ? ? 判断两个训练样本是否为数组,如果不是,将其转换为数组
? ? ? if(is.matrix(trnX1) == TRUE) trnX1<-as.matrix(trnX1);
? ? ? if(is.matrix(trnX2) == TRUE) trnX2<-as.matrix(trnX2);
# ? ? 获取测试样本的行数
? ? ? nrowx<-nrow(testX);
# ? ? 需要增加belong解释
? ? ? belong<-matrix(rep(0,nrowx),nrow=1,byrow=T,dimnames=list("belong",1:nrowx));
# ? ? 获取两个训练样本的列均值mu1,mu2
? ? ? mu1<-colMeans(trnX1);mu2<-colMeans(trnX2);
# ? ? 以下是通过统计学公式计算得来
# ? ? 1.如果两个样本的协方差相同(TRUE or T),则对两个训练样本合并rbind并计算方差(距离的平方)
# ? ? 计算测试样本与样本1和样本2的马氏距离的平方之差,判定测试样本属于哪一类样本
? ? ? if(var.equal == TRUE ||var.equal == T) {
? ? ? ? s<-var(rbind(trnX1,trnX2));
? ? ? ? w<-mahalanobis(testX,mu1,s) - mahalanobis(testX,mu2,s);
? ? ? } else {
# ? ? 2.否则(两样本协方差不相同),则分别计算两样本的方差
? ? ? ? s1<-var(trnX1);s2<-var(trnX2);
# ? ? 再计算测试样本与两训练样本之间的距离差 ? ? ? ?
? ? ? ? w<-mahalanobis(testX,mu1,s1)-mahalanobis(testX,mu2,s2);
? ? ? } ?
# ? ? 开始循环(循环数组的每一行,如果距离差大于0,则属于第一类赋值1,否则小于等于0,属于第二类样本赋值2)
? ? ? for (i in 1:nrowx) {
? ? ? ? if(w[i] > 0) {
? ? ? ? ? belong[i] <- 1;
? ? ? ? } else {
? ? ? ? ? belong[i] <- 2;
? ? ? ? }
? ? ? }
# ? ? 返回赋值结果belong(vector)
? ? ? belong;
? }
?
# 测试数据
data.create<-function () {
? classX1<-data.frame(
? ? x1=c(39.00,39.00, 47.00, 47.00, 32.00, 6.0, 113.00, 52.00,52.00,113.00,172.00,172.00),
? ? x2=c(1.00, 1.00, 1.00, 1.00, 2.00, 1.0, 3.50,1.00,3.50, 0.00, 1.00, 1.50),
? ? x3=c(6.00, 6.00, 6.00, 6.00, 7.50, 7.0, 6.00,6.00,7.50, 7.50, 3.50, 3.00)
? )
?
? classX2<-data.frame(
? ? x1=c(32.0 ,32.00, 32.00, 11.0, 8.00, 8.00, 8.00,161.00, 161.0, 6.0, 6.0, 6.0, 6.00,113.00,113.00, 52.00, 52.00, 97.00, 97.00,89.00,56.00,172.00,283.00),
? ? x2=c(1.00, 2.00, 2.50, 4.5, 4.50, 6.00, 1.50, 1.50, 0.50, 3.5, 1.0, 1.0, 2.50, 3.50, 3.50, 1.00,1.00, 0.00, 2.50, 0.00, 1.50, 1.00, 1.00),
? ? x3=c(5.00, 9.00, 4.00, 7.5, 4.50, 7.50, 6.00, 4.0,2.50, 4.0, 3.0, 6.0, 3.00, 4.50, 4.50, 6.0,7.50, 6.00, 6.00, 6.00, 6.00, 3.50, 4.50)
? )
? print("hello");
? classX1;classX2;
}
?
# 多分类的距离判断
?
distinguish.distance<-
? function (TrnX, TrnG, TstX = NULL, var.equal = FALSE){
# ? TrnG:因子变量,输入样本的分类情况
# ? 1.如果因子变量不是因子,记录训练样本和因子变量的行数,用于生成因子变量
# ? 2.将训练样本和因子变量合并赋值给测试样本变量
# ? 3.对因子变量赋值
? ? if ( is.factor(TrnG) == FALSE){
? ? ? mx<-nrow(TrnX); mg<-nrow(TrnG)
? ? ? TrnX<-rbind(TrnX, TrnG)
? ? ? TrnG<-factor(rep(1:2, c(mx, mg)))
? ? }
# ? 如果测试样本为空,将训练样本赋值给测试样本
? ? if (is.null(TstX) == TRUE) TstX<-TrnX
# ? 如果测试样本为向量,转换为数组并转置;如果转换后依然不是数组,再次转置赋值给测试样本
? ? if (is.vector(TstX) == TRUE)? TstX<-t(as.matrix(TstX))
? ? else if (is.matrix(TstX) != TRUE)
? ? ? TstX<-as.matrix(TstX)
# ? 如果训练样本不是数组,转换为数组
? ? if (is.matrix(TrnX) != TRUE) TrnX<-as.matrix(TrnX)
# ? 记录测试数组的行数 ?
? ? nx<-nrow(TstX)
# ? 给判别结果赋初值,形成数组一行多列 ?
? ? blong<-matrix(rep(0, nx), nrow=1, dimnames=list("blong", 1:nx))
# ? 获取训练样本因子的长度
? ? g<-length(levels(TrnG))
#? ?
? ? mu<-matrix(0, nrow=g, ncol=ncol(TrnX))
? ? for (i in 1:g)
? ? ? mu[i,]<-colMeans(TrnX[TrnG==i,])?
? ? D<-matrix(0, nrow=g, ncol=nx)
? ? if (var.equal == TRUE? || var.equal == T){
? ? ? for (i in 1:g)
? ? ? ? D[i,]<- mahalanobis(TstX, mu[i,], var(TrnX))
? ? }
? ? else{
? ? ? for (i in 1:g)
? ? ? ? D[i,]<- mahalanobis(TstX, mu[i,], var(TrnX[TrnG==i,]))
? ? }
? ? for (j in 1:nx){
? ? ? dmin<-Inf
? ? ? for (i in 1:g)
? ? ? ? if (D[i,j]<dmin){
? ? ? ? ? dmin<-D[i,j]; blong[j]<-i
? ? ? ? }
? ? }
? ? blong
}
data(iris);
X<-iris[,1:4];
G<-gl(3,50);
distinguish.distance(X,G);