R软件做判别分析

1、R 软件做判别分析:1. 距离判别（1）两总体discriminiant.distance <- function(TrnX1, TrnX2, TstX = NULL, var.equal = FALSE)if (is.null(TstX) = TRUE) TstX <- rbind(TrnX1,TrnX2)if (is.vector(TstX) = TRUE) TstX <- t(as.matrix(TstX)else if (is.matrix(TstX) != TRUE)TstX <- as.matrix(TstX)if (is.matrix(TrnX1) !=

2、TRUE) TrnX1 <- as.matrix(TrnX1)if (is.matrix(TrnX2) != TRUE) TrnX2 <- as.matrix(TrnX2)nx <- nrow(TstX)blong <- matrix(rep(0, nx), nrow=1, byrow=TRUE, dimnames=list("blong", 1:nx)mu1 <- colMeans(TrnX1); mu2 <- colMeans(TrnX2)if (var.equal = TRUE | var.equal = T)S <- var

3、(rbind(TrnX1,TrnX2)w <- mahalanobis(TstX, mu2, S)- mahalanobis(TstX, mu1, S)elseS1 <-var(TrnX1); S2 <- var(TrnX2)w <- mahalanobis(TstX, mu2, S2)- mahalanobis(TstX, mu1, S1)for (i in 1:nx)if (wi > 0) blongi <- 1elseblongi <- 2blong 例1: 数据classX1<-data.frame(x1=c(6.60, 6.60, 6.

4、10, 6.10, 8.40, 7.2, 8.40, 7.50,7.50, 8.30, 7.80, 7.80),x2=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),x3=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),x4=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),x5=c(6.0

5、0, 12.00, 6.00, 12.00, 19.00, 28.0, 18.00, 12.00,6.00, 35.00, 14.00, 15.00),x6=c(0.12, 0.12, 0.08, 0.08, 0.35, 0.3, 0.15, 0.16,0.16, 0.12, 0.21, 0.21),x7=c(20.00,20.00, 12.00, 12.00, 75.00, 30.0, 75.00, 40.00,40.00,180.00, 45.00, 45.00)classX2<-data.frame(x1=c(8.40, 8.40, 8.40, 6.3, 7.00, 7.00, 7

6、.00, 8.30,8.30, 7.2, 7.2, 7.2, 5.50, 8.40, 8.40, 7.50,7.50, 8.30, 8.30, 8.30, 8.30, 7.80, 7.80),x2=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),x3=c(1.00, 2.00, 2.50, 4.5, 4.50, 6.00, 1.50, 1.50,0.5

7、0, 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),x4=c(5.00, 9.00, 4.00, 7.5, 4.50, 7.50, 6.00, 4.00,2.50, 4.0, 3.0, 6.0, 3.00, 4.50, 4.50, 6.00,7.50, 6.00, 6.00, 6.00, 6.00, 3.50, 4.50),x5=c(4.00, 10.00, 10.00, 3.0, 9.00, 4.00, 1.00, 4.00,1.00, 12.0, 3.0, 5.0, 7.00,

8、 6.00, 8.00, 6.00,8.00, 5.00, 5.00,10.00,13.00, 6.00, 6.00),x6=c(0.35, 0.35, 0.35, 0.2, 0.25, 0.25, 0.25, 0.08,0.08, 0.30, 0.3, 0.3, 0.18, 0.15, 0.15, 0.16,0.16, 0.15, 0.15, 0.16, 0.25, 0.21, 0.18),x7=c(75.00,75.00, 75.00, 15.0,30.00, 30.00, 30.00, 70.00,70.00, 30.0, 30.0, 30.0,18.00, 75.00, 75.00,

9、40.00,40.00,180.00,180.00,180.00,180.00,45.00,45.00)source("discriminiant.distance.R")discriminiant.distance(classX1, classX2, var.equal=TRUE)discriminiant.distance(classX1, classX2)课本应用多元统计分析P182 例5.1.1TrnX1<-data.frame(x1=c(13.85, 22.31, 28.82,15.29, 28.79),x2=c(2.79, 4.67, 4.63, 3.54

10、, 4.90),x3= c(7.80, 12.31, 16.18, 7.50, 16.12),x4=c(49.60, 47.80, 62.15, 43.20, 58.10)TrnX2<-data.frame(x1=c(2.18, 3.85, 11.40, 3.66, 12.10),x2=c(1.06, 0.80, 0.00, 2.42, 0.00),x3=c(1.22, 4.06, 3.50, 2.14,5.68),x4=c(20.60, 47.10, 0.00, 15.10, 0.00)TstX<-data.frame(x1=c(8.85,28.60,20.70,7.90,3.1

11、9,12.40,16.80,15.00),x2=c(3.38, 2.40,6.70,2.40,3.20,5.10,3.40,2.70),x3=c(5.17, 1.20,7.60,4.30,1.43,4.43,2.31,5.02),x4=c(26.10,127.00,30.20,33.20,9.90,24.60,31.30,64.00)discriminiant.distance(TrnX1,TrnX2,TstX,var.equal=TRUE)discriminiant.distance(TrnX1,TrnX2,TstX,var.equal=FALSE)#默认var.equal=FALSE即是：

12、discriminiant.distance(TrnX1,TrnX2,TstX)(2)多总体distinguish.distance <- function(TrnX, TrnG, TstX = NULL, var.equal = FALSE)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 <- TrnXi

13、f (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 <- ma

14、trix(0, nrow=g, ncol=ncol(TrnX)for (i in 1:g)mui, <- colMeans(TrnXTrnG=i,)D <-matrix(0, nrow=g, ncol=nx)if (var.equal = TRUE | var.equal = T)for (i in 1:g)Di, <- mahalanobis(TstX, mui, var(TrnX)elsefor (i in 1:g)Di, <- mahalanobis(TstX, mui, var(TrnXTrnG=i,)for (j in 1:nx)dmin <- Inff

15、or (i in 1:g)if (Di,j < dmin)dmin <- Di,j; blongj <- iblong多总体距离判别:(要求数据矩阵或数据框)X<- matrix(c(34.16, 7.44, 1.12, 7.87, 95.19, 69.30, 33.06, 6.34, 1.08, 6.77, 94.08, 69.70, 36.26, 9.24, 1.04, 8.97, 97.30, 68.80, 40.17, 13.45, 1.43, 13.88,101.2, 66.2, 50.06, 23.03, 2.83, 23.74, 112.52, 63.3,

16、 33.24, 6.24, 1.18, 22.9, 160.01, 65.4, 32.22, 4.22, 1.06, 20.7, 124.7, 68.7, 41.15, 10.08, 2.32, 32.84, 172.06, 65.85, 53.04, 25.74, 4.06, 34.87, 152.03, 63.5, 38.03, 11.2, 6.07, 27.84, 146.32, 66.8, 34.03, 5.41, 0.07, 5.2, 90.10, 69.5, 32.11, 3.02, 0.09, 3.14,85.15, 70.8, 44.12, 15.12, 1.08, 15.15

17、, 103.12, 64.8, 54.17, 25.03, 2.11, 25.15, 110.14, 63.7, 28.07, 2.01, 0.07, 3.02, 81.22, 68.3),nrow=15,ncol=6, byrow=TRUE)G<-gl(3,5)Y<- matrix(c(50.22, 6.66, 1.08, 22.54,170.6, 65.2, 34.64, 7.33, 1.11, 7.78, 95.16, 69.3, 33.42, 6.22, 1.12, 22.95, 160.31, 68.3, 44.02,15.36,1.07,16.45,105.3,64.2

18、), nrow=4, ncol=6, byrow=TRUE)distinguish.distance(X, G, Y,var.equal=TRUE)Bayes判别:(1)两总体discriminiant.bayes <- function(TrnX1, TrnX2, rate = 1, TstX = NULL, var.equal = FALSE)if (is.null(TstX) = TRUE) TstX<-rbind(TrnX1,TrnX2)if (is.vector(TstX) = TRUE) TstX <- t(as.matrix(TstX)else if (is.m

19、atrix(TstX) != TRUE)TstX <- as.matrix(TstX)if (is.matrix(TrnX1)!= TRUE) TrnX1 <- as.matrix(TrnX1)if (is.matrix(TrnX2)!= TRUE) TrnX2 <- as.matrix(TrnX2)nx <- nrow(TstX)blong <- matrix(rep(0, nx), nrow=1, byrow=TRUE,dimnames=list("blong", 1:nx)mu1 <- colMeans(TrnX1); mu2 &l

20、t;- colMeans(TrnX2)if (var.equal = TRUE | var.equal = T)S <- var(rbind(TrnX1,TrnX2); beta <- 2*log(rate)w <- mahalanobis(TstX, mu2, S)- mahalanobis(TstX, mu1, S)elseS1 <- var(TrnX1); S2 <- var(TrnX2)beta <- 2*log(rate) + log(det(S1)/det(S2)w <- mahalanobis(TstX, mu2, S2)- mahala

21、nobis(TstX, mu1, S2)for (i in 1:nx)if (wi > beta)blongi <- 1elseblongi <- 2blong例2: 两总体TrnX1<-matrix(c(24.8, 24.1, 26.6, 23.5, 25.5, 27.4,-2.0, -2.4, -3.0, -1.9, -2.1, -3.1),ncol=2)TrnX2<-matrix(c(22.1, 21.6, 22.0, 22.8, 22.7, 21.5, 22.1, 21.4,-0.7, -1.4, -0.8, -1.6, -1.5, -1.0, -1.2,

22、 -1.3),ncol=2)source("discriminiant.bayes.R")discriminiant.bayes(TrnX1, TrnX2, rate=8/6, var.equal=TRUE)(2)多总体distinguish.bayes <- function(TrnX, TrnG, p = rep(1, length(levels(TrnG),TstX =NULL, var.equal = FALSE)if ( is.factor(TrnG) = FALSE)mx <- nrow(TrnX); mg <- nrow(TrnG)TrnX

23、<- rbind(TrnX, TrnG)TrnG <- factor(rep(1:2, c(mx, mg)if (is.null(TstX) = TRUE) TstX <- TrnXif (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 <-

24、 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)mui, <- colMeans(TrnXTrnG=i,)D <- matrix(0, nrow=g, ncol=nx)if (var.equal = TRUE | var.equal = T)for (i in 1:g)d2 <- mahalanobis(TstX, mui, var

25、(TrnX)Di, <- d2 - 2*log(pi)Elsefor (i in 1:g)S <- var(TrnXTrnG=i,)d2 <- mahalanobis(TstX, mui, S)Di, <- d2 - 2*log(pi)-log(det(S)for (j in 1:nx)dmin <- Inffor (i in 1:g)if (Di,j < dmin)dmin <- Di,j; blongj <- iblong X<- matrix(c(34.16,7.44,1.12,7.87,95.19,69.30,33.06,6.34,

26、1.08,6.77,94.08,69.70,36.26,9.24,1.04,8.97,97.30,68.80,40.17,13.45,1.43,13.88,101.2,66.2,50.06,23.03,2.83,23.74,112.52,63.3,33.24,6.24,1.18,22.9,160.01,65.4,32.22,4.22,1.06,20.7,124.7,68.7,41.15,10.08,2.32,32.84,172.06,65.85,53.04,25.74,4.06,34.87,152.03,63.5,38.03,11.2,6.07,27.84,146.32,66.8,34.03,

27、5.41,0.07,5.2,90.10,69.5,32.11,3.02,0.09,3.14,85.15,70.8,44.12,15.12,1.08,15.15,103.12,64.8,54.17,25.03,2.11,25.15,110.14,63.7,28.07,2.01,0.07,3.02,81.22,68.3),nrow=15,ncol=6,byrow=TRUE)G<-gl(3,5)Y<- matrix(c(50.22, 6.66, 1.08, 22.54,170.6, 65.2, 34.64, 7.33, 1.11, 7.78, 95.16, 69.3, 33.42, 6.

28、22, 1.12, 22.95, 160.31, 68.3, 44.02,15.36,1.07,16.45,105.3,64.2),nrow=4,ncol=6,byrow=TRUE)distinguish.bayes(X,G,p= rep(1, length(levels(G),Y,var.equal =TRUE)(3)Fisher判别discriminiant.fisher <- function(TrnX1, TrnX2, TstX = NULL)if (is.null(TstX) = TRUE) TstX <- rbind(TrnX1,TrnX2)if (is.vector(

29、TstX) = TRUE) TstX <- t(as.matrix(TstX)else if (is.matrix(TstX) != TRUE)TstX <- as.matrix(TstX)if (is.matrix(TrnX1) != TRUE) TrnX1 <- as.matrix(TrnX1)if (is.matrix(TrnX2) != TRUE) TrnX2 <- as.matrix(TrnX2)nx <- nrow(TstX)blong <- matrix(rep(0, nx), nrow=1, byrow=TRUE, dimnames=list

30、("blong", 1:nx)n1 <- nrow(TrnX1); n2 <- nrow(TrnX2)mu1 <- colMeans(TrnX1); mu2 <- colMeans(TrnX2)S <- (n1-1)*var(TrnX1) + (n2-1)*var(TrnX2)mu <- n1/(n1+n2)*mu1 + n2/(n1+n2)*mu2w <- (TstX-rep(1,nx) %o% mu) %*% solve(S, mu2-mu1);for (i in 1:nx)if (wi <= 0)blongi <-

31、 1elseblongi <- 2blong多总体Fisher判别:(要求数据框)library(MASS)X<- matrix(c(34.16,7.44,1.12,7.87,95.19,69.30,33.06,6.34,1.08,6.77,94.08,69.70,36.26,9.24,1.04,8.97,97.30,68.80,40.17,13.45,1.43,13.88,101.2,66.2,50.06,23.03,2.83,23.74,112.52,63.3,33.24,6.24,1.18,22.9,160.01,65.4,32.22,4.22,1.06,20.7,124.7

32、,68.7,41.15,10.08,2.32,32.84,172.06,65.85,53.04,25.74,4.06,34.87,152.03,63.5,38.03,11.2,6.07,27.84,146.32,66.8,34.03,5.41,0.07,5.2,90.10,69.5,32.11,3.02,0.09,3.14,85.15,70.8,44.12,15.12,1.08,15.15,103.12,64.8,54.17,25.03,2.11,25.15,110.14,63.7,28.07,2.01,0.07,3.02, 81.22,68.3), nrow=15, ncol=6, byro

33、w=TRUE)Y<- data.frame(X, Sp = rep(c("1","2","3"), rep(5,3)Train<-data.frame(matrix(c(50.22,6.66,1.08,22.54,170.6,65.2,34.64,7.33,1.11,7.78,95.16,69.3,33.42,6.22,1.12,22.95,160.31,68.3,44.02,15.36,1.07,16.45,105.3,64.2),nrow=4,ncol=6,byrow=TRUE)z <- lda(Sp ., Y)predict(z, Train )多总体Bayes判别:(要求数据矩阵或数据框)X<- matrix(c(34.16,7.44,1.12,7.87,95.19,69.30,33.06,6.34,1.08,6.77,94.08,69.70,36.26,9.24,1.04,8.97,97.30,68.80,40.17,13.45,1.43,13.88,101.2,66.2,50.06,23.03,2.83,23.74,112.52,63.3,33.24,6.24,1.