# Input data d <- read.table("inp_iris.txt",header=TRUE) Ldat1 <- d[1:30,1:4] Ldat2 <- d[51:80,1:4] Ldat3 <- d[101:130,1:4] Edat1 <- d[31:50,1:4] Edat2 <- d[81:100,1:4] Edat3 <- d[131:150,1:4] # Use of R-function 'mahalanobis' m1 <- apply(Ldat1,2,mean) m2 <- apply(Ldat2,2,mean) m3 <- apply(Ldat3,2,mean) s1 <- var(Ldat1) s2 <- var(Ldat2) s3 <- var(Ldat3) D11 <- mahalanobis(Edat1,m1,s1) D21 <- mahalanobis(Edat1,m2,s2) D31 <- mahalanobis(Edat1,m3,s3) D12 <- mahalanobis(Edat2,m1,s1) D22 <- mahalanobis(Edat2,m2,s2) D32 <- mahalanobis(Edat2,m3,s3) D13 <- mahalanobis(Edat3,m1,s1) D23 <- mahalanobis(Edat3,m2,s2) D33 <- mahalanobis(Edat3,m3,s3) # Original function maha_org <- function(ldat,edat){ mm <- ncol(ldat) nn <- nrow(edat) inv <- solve(var(ldat)) datm <- sweep(edat,2,apply(ldat,2,mean),FUN="-") dis <- numeric(0) for(i in 1:nn){ w1 <- datm[i,1:mm] w2 <- as.matrix(inv) %*% t(as.matrix(w1)) w3 <- as.matrix(w1) %*% as.matrix(w2) dis[i]<-(w3) } return (dis) } dis11 <- maha_org(Ldat1,Edat1) dis21 <- maha_org(Ldat2,Edat1) dis31 <- maha_org(Ldat3,Edat1) dis12 <- maha_org(Ldat1,Edat2) dis22 <- maha_org(Ldat2,Edat2) dis32 <- maha_org(Ldat3,Edat2) dis13 <- maha_org(Ldat1,Edat3) dis23 <- maha_org(Ldat2,Edat3) dis33 <- maha_org(Ldat3,Edat3) # Probability pp11 <- pchisq(dis11, ncol(Edat1), lower.tail=FALSE) pp21 <- pchisq(dis21, ncol(Edat1), lower.tail=FALSE) pp31 <- pchisq(dis31, ncol(Edat1), lower.tail=FALSE) pp12 <- pchisq(dis12, ncol(Edat2), lower.tail=FALSE) pp22 <- pchisq(dis22, ncol(Edat2), lower.tail=FALSE) pp32 <- pchisq(dis32, ncol(Edat2), lower.tail=FALSE) pp13 <- pchisq(dis13, ncol(Edat3), lower.tail=FALSE) pp23 <- pchisq(dis23, ncol(Edat3), lower.tail=FALSE) pp33 <- pchisq(dis33, ncol(Edat3), lower.tail=FALSE) # print out data cbind(D11,D21,D31) cbind(dis11,dis21,dis31) cbind(D12,D22,D32) cbind(dis12,dis22,dis32) cbind(D13,D23,D33) cbind(dis13,dis23,dis33) cbind(pp11,pp21,pp31) cbind(pp12,pp22,pp32) cbind(pp13,pp23,pp33) # par(mfrow=c(1,2)) par(pty="s") par(ps=10) smax<-1000 plot( 0.1,0.1, type="n",log="xy", xlim=c(0.1,smax),ylim=c(0.1,smax), pch=21,bg="red", xlab="D squared from center of Group 1", ylab="D squared from center of Group 2" ) polygon(c(0.1,0.1,smax,0.1),c(0.1,smax,smax,0.1),col="springgreen") polygon(c(0.1,smax,smax,0.1),c(0.1,smax,0.1,0.1),col="aliceblue") points(dis11,dis12,pch=0,cex=1.0,col="maroon") points(dis21,dis22,pch=1,cex=1.0,col="blue") text(0.5,5,"Group 1") text(5,0.5,"Group 2") # par(pty="s") plot( 0.1,0.1, type="n",log="xy", xlim=c(0.1,smax),ylim=c(0.1,smax), pch=21,bg="red", xlab="D squared from center of Group 2", ylab="D squared from center of Group 3" ) polygon(c(0.1,0.1,smax,0.1),c(0.1,smax,smax,0.1),col="aliceblue") polygon(c(0.1,smax,smax,0.1),c(0.1,smax,0.1,0.1),col="lightyellow") points(dis22,dis23,pch=1,cex=1.0,col="blue") points(dis32,dis33,pch=2,cex=0.8,col="red") text(0.5,5,"Group 2") text(5,0.5,"Group 3") #