统计建模与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);