R软件做判别分析.doc

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) != 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(rbind(TrnX1,TrnX2)w - mahalanobis(TstX, mu2, S)- mahalanobis(TstX, mu1, S)elseS1 -var(TrnX1); S2 - var(TrnX2)w 0) blongi - 1elseblongi - 2blong 例1: 数据classX1-data.frame(x1=c(6.60, 6.60, 6.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.00, 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.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.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),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, 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, 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)TrnX1-data.frame(X1=c(13.85, 22.31, 28.82,15.29, 28.79),X2=c(2.79, 4.67, 4.63, 3.54, 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)TrnX-data.frame(X1=c(8.85, 28.60),X2=c(3.38, 2.40),X3=c(5.17, 1.20),X4=c(26.10, 127.0)discriminiant.distance(TrnX1, TrnX2)(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 - 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 - 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)Di, - mahalanobis(TstX, mui, var(TrnX)elsefor (i in 1:g)Di, - mahalanobis(TstX, mui, var(TrnXTrnG=i,)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, 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, 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.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.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.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(blong, 1:nx)mu1 - colMeans(TrnX1); mu2 - 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 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, -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 - 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 - 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(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, 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, 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.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(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(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 - 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, 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, byrow=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.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