-
Notifications
You must be signed in to change notification settings - Fork 0
/
h1code.R
84 lines (73 loc) · 1.77 KB
/
h1code.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
screeplot <- function(p) {
e <- p$sdev ^ 2
e <- e / sum(e)
plot(
1:length(e),
e,
xlab = "Component number",
pch = 20,
ylab = "Variance proportion",
main = "Scree plot",
axes = F,
ylim = c(0, max(e)*1.04)
)
lines(1:length(e), e)
axis(1, at = 1:length(e))
axis(2)
}
sumPartition <- function(dat,z){
k <-length(unique(z))
nobs <- tapply(1:nrow(dat),z,length)
withinSS <- tapply(1:nrow(dat),z,function(i){
if (length(i)==1) 0
else {x<- dat[i,]
xm <- scale(x, scale=F)
sum(xm^2)
}
})
aveD <- tapply(1:nrow(dat),z,function(i){
if (length(i)==1) 0
else {x<- dat[i,]
xm <- scale(x, scale=F)
xm <- apply(xm, 1, function(y) sqrt(sum(y^2)))
mean(xm)
}
})
maxD <- tapply(1:nrow(dat),z,function(i){
if (length(i)==1) 0
else {x<- dat[i,]
xm <- scale(x, scale=F)
xm <- apply(xm, 1, function(y) sqrt(sum(y^2)))
max(xm)
}
})
part<- data.frame("N.obs"=nobs, "Within clus SS" = withinSS, "Ave dist Centroid" = aveD,
"Max dist centroid" =maxD)
rownames(part)<- paste("Cluster", 1:k)
cCen <- cbind(simplify2array(tapply(1:nrow(dat),z,function(i){
x<- dat[i,]
if (length(i)==1) x
else {
colMeans(x)
}
})), colMeans(dat))
colnames(cCen)<- c(paste("Cluster", 1:k), "Grand centrd")
cCenD <- as.matrix(dist(t(cCen[, -ncol(cCen)])))
cat("Final Partition\n")
cat("\n")
cat(paste("Number of clusters ", k))
cat("\n")
cat("\n")
print(part, quote=F)
cat("\n")
cat("\n")
cat("Cluster centroids\n")
cat("\n")
print(cCen,quote=F)
cat("\n")
cat("\n")
cat("Distances between Cluster centroids\n")
cat("\n")
print(cCenD,quote=F)
list(withinSS=withinSS, aveD=aveD, centroids=cCen, centroidDist=cCenD)
}