Nothing
R/fpc.R
In fpc: Flexible Procedures for Clustering
Defines functions
print.summary.mfpc
print.mfpc
plot.mfpc
fpclusters.mfpc
summary.mfpc
mahalconf
wfu
cmahal
solvecov
simmatrix
con.comp
fpclusters
Documented in
cmahal
con.comp
fpclusters
fpclusters.mfpc
mahalconf
plot.mfpc
print.mfpc
print.summary.mfpc
simmatrix
solvecov
summary.mfpc
wfu
# .packageName <- "fpc"
# .First.lib <- function(lib,pkg){
# require(cluster)
# # if (as.numeric(R.version$major)<=1 & as.numeric(R.version$minor)<9){
# # require(lqs)
# # require(mva)
# # }
# require(MASS)
# }
#
#
# generic functions
#
fpclusters <- function(object, ...)
UseMethod("fpclusters")
#
# general/set distance utilities
#
#
# sseg= position of distance between i and j in smatrix/imatrix
sseg <- function (i,j){
if (i>=j)
out <- (i-1)*i/2+j
else
out <- (j-1)*j/2+i
out
}
# Connectivity components of a graph by depth-first search
# Input: Boolean coincidence matrix, output: Vector of component numbers
con.comp <- function(comat){
nc <- ncol(comat)
ccn <- rep(0, times=nc) # con.comp number for each point
fhist <- rep(FALSE, times=nc) # indicator if point had been under consideration
stn <- 0 # current cc number
pn <- 1 # point no. to which similar objects are looked for
while(pn>0){
stn <- stn+1
repeat{
sm <- 0 # smallest new point no.
ccn[pn] <- stn
fhist[pn] <- TRUE
if(nc>1)
{
for(i in 2:nc)
{
# cat(i, pn, ccn[i], comat[i,pn],"\n")
if((ccn[i]==0) & (comat[i,pn]))
ccn[i] <- stn
if ((sm==0) & (ccn[i]==stn) & (fhist[i]==FALSE))
sm <- i
} # for i
# cat("stn=", stn, "sm=", sm, "\n")
} # if nc>1
if (sm>0)
pn <- sm
else
break
} # repeat
# print("repeat terminated")
pn <- 0
i <- 2
while(i<=nc){
if(ccn[i]==0){
pn <- i
i <- nc
} # if
i <- i+1
} # while i
} # while pn>0 (stn-loop)
ccn
}
# Intersection matrix between significant fpc groups
simmatrix <- function(fpcobj){
sim <- c()
# print(stn)
for(i in 1:fpcobj$stn)
for(j in i:fpcobj$stn){
# cat(i,j,fpcobj$sfpc[fpcobj$sto[i]],fpcobj$sfpc[fpcobj$sto[j]])
sim[sseg(i,j)] <- fpcobj$imatrix[sseg(fpcobj$sfpc[fpcobj$sto[i]],fpcobj$sfpc[fpcobj$sto[j]])]
}
sim
}
#
# fixmahal utilities
#
# Weighted ML-covariance matrix
cov.wml <- function (x, wt = rep(1/nrow(x), nrow(x)),
cor = FALSE, center = TRUE)
{
if (is.data.frame(x))
x <- as.matrix(x)
else if (!is.matrix(x))
stop("x must be a matrix or a data frame")
if (!all(is.finite(x)))
stop("x must contain finite values only")
n <- nrow(x)
if (with.wt <- !missing(wt)) {
if (length(wt) != n)
stop("length of wt must equal the number of rows in x")
if (any(wt < 0) || (s <- sum(wt)) == 0)
stop("weights must be non-negative and not all zero")
wt <- wt/s
}
if (is.logical(center)) {
center <- if (center)
colSums(wt * x)
else 0
}
else {
if (length(center) != ncol(x))
stop("length of center must equal the number of columns in x")
}
x <- sqrt(wt) * sweep(x, 2, center)
cov <- (t(x) %*% x)
y <- list(cov = cov, center = center, n.obs = n)
if (with.wt)
y$wt <- wt
if (cor) {
sdinv <- diag(1/sqrt(diag(cov)), nrow(cov))
y$cor <- sdinv %*% cov %*% sdinv
}
y
}
# inversion of cov-matrices; if singular, Eigenvalues below 1/cmax are set
# to 1/cmax.
solvecov <- function(m, cmax=1e10){
options(show.error.messages = FALSE)
covinv <- try(solve(m))
if(!inherits(covinv,"try-error"))
coll=FALSE
else{
p <- nrow(m)
cove <- eigen(m, symmetric=TRUE)
coll <- TRUE
if (min(cove$values)<1/cmax){
covewi <- diag(p)
for (i in 1:p)
if (cove$values[i]<1/cmax) covewi[i,i] <- cmax
else covewi[i,i] <- 1/cove$values[i]
}
else covewi <- diag(1/cove$values,nrow=length(cove$values))
# cat("covewi ",covewi)
covinv <- cove$vectors %*% covewi %*% t(cove$vectors)
}
options(show.error.messages = TRUE)
out <- list(inv=covinv,coll=coll)
}
# Generation of ca by "recursion for convergence"
cmahal <- function(n,p,nmin,cmin,nc1,c1=cmin,q=1){
cs <- rep(c1,max(n,nc1))
for (i in (nc1-1):nmin)
cs[i] <- cs[i+1]+q*(cs[i+1]-p)/(i+1-q)
cs[1:(nmin-1)] <- rep(cs[nmin],nmin-1)
if (n>nc1)
for (i in (nc1+1):n){
cs[i] <- cs[i-1]-q*(cs[i-1]-p)/i
if (cs[i]<=cmin){
cs[i] <- cmin
i <- n
}
}
cs
}
wfu <- function(md,ca,ca2,a1=1/(ca-ca2),a0=-a1*ca2){
v <- rep(0,length(md))
v[md<=ca] <- 1
mind <- (md>ca) & (md<=ca2)
v[mind] <- a1*md[mind]+a0
v
}
# mahalanofix=vector of mahalanobis distances from all points to
# center of points indexed by gv (default=all)
mahalanofix <- function (x, n=nrow(as.matrix(x)), p=ncol(as.matrix(x)),
gv=rep(1, times=n),
cmax=1e+10, method="ml") {
# print("Mahalanofix")
gv <- as.logical(gv)
ng <- sum(gv)
xg <- x[gv,1:p,drop=FALSE]
if (method=="ml"){
mg <- colMeans(xg)
if (method=="ml")
covg <- (ng-1)*cov(xg)/ng
else
covg <- cov(xg)
}
else{
# if (as.numeric(R.version$major)<=1 & as.numeric(R.version$minor)<9)
# require(lqs)
# require(MASS)
grob <- cov.rob(xg, method=method)
mg <- grob$center
covg <- grob$cov
}
cm <- solvecov(covg,cmax=cmax)
# print("robmahal")
md <- mahalanobis(x,mg,cm$inv,inverted=TRUE)
# print("end robmahal")
out <- list(md=md, mg=mg, covg=covg, covinv=cm$inv, coll=cm$coll)
}
# vector of mahalanobis distances from all points to
# center of points indexed by "fuzzy" gv (default=all)
mahalanofuz <- function (x, n=nrow(as.matrix(x)), p=ncol(as.matrix(x)),
gv=rep(1, times=n),
cmax=1e+10) {
# print("Mahalanofix")
mg <- cov.wml(x,gv)
covg <- mg$cov
mg <- mg$center
cm <- solvecov(covg, cmax=cmax)
# print("Mahalanobis")
md <- mahalanobis(x,mg,cm$inv,inverted=TRUE)
# print("End Mahalanofix")
out <- list(md=md, mg=mg, covg=covg, covinv=cm$inv, coll=cm$coll)
}
# Start configuration for fixmahal, point no. no
# startn points, cov/center of whole data
mahalconf <- function(x, no, startn, covall, plot){
p <- ncol(x)
n <- nrow(x)
gv <- rep(FALSE, n)
# print("mahalconf")
# print(p)
# print(n)
# print(no)
# print(startn)
# print(covall)
om <- order(mahalanobis(x,center=x[no,],cov=solvecov(covall)$inv, inverted=TRUE))
gv[om[1:(p+1)]] <- TRUE
# cat("conf ",gv," ",om,"\n")
if (plot=="start" || plot=="both")
plot(x,col=1+gv+((1:n)==no),pch=1+gv+((1:n)==no))
if (startn>p+1)
for (pn in (p+2):startn){
om <- order(mahalanofix(x,n,p,gv)$md)
j <- 0
while (sum(gv)<pn){
j <- j+1
if (!gv[om[j]])
gv[om[j]] <- TRUE
}
}
# print(sum(gv))
gv
}
#
# fixmahal iteration
#
fpmi <- function (dat, n=nrow(as.matrix(dat)), p=ncol(as.matrix(dat)),
gv, ca, ca2, method="ml", plot,
maxit=5*n, iter=n*1e-6) {
gva <- rep(0, times=n)
change <- TRUE
coll <- FALSE
ir <- 0
while (change) {
ir <- ir+1
if (method=="fuzzy")
change <- sum((gv-gva)^2>iter)
else
change <- !identical(gv,gva)
if (ir>maxit){
change <- FALSE
warning("Maximum number of iterations exceeded.")
}
gva <- gv
# print("mahalanobis")
if (method=="fuzzy")
mfix <- mahalanofuz(x=dat, n=n, p=p, gv=gv)
else
mfix <- mahalanofix(x=dat, n=n, p=p, gv=gv, method=method)
# print(ca)
md <- mfix$md
mg <- mfix$mg
if (mfix$coll) coll=TRUE
covg <- mfix$covg
# print(md)
if (method=="fuzzy"){
gv <- wfu(md, ca, ca2)
# print(gv)
cn <- ca
# plot(x,col=1+gv)
}
else{
ng <- sum(gv)
if (length(ca)==1)
cn <- ca
else
cn <- ca[ng]
# print(ng)
# cat("Target function = ",
# det(covg*(ng-1)/ng)*prod(exp((2-ca)/(1:ng))), "\n")
gv <- md<=cn
}
# print(gv)
if (plot=="iteration" || plot=="both"){
if (method=="fuzzy")
plotgv <- as.integer(gv>0.99)+
3*as.integer(gv>0.01 & gv<=0.99)
else
plotgv <- gv
plot(dat,col=1+plotgv,pch=1+plotgv)
}
} # while change
out <- list(mg=mg, covg=covg, md=md, gv=gv, coll=coll, method=method, ca=cn)
out
}
#
# fixmahal main and output
#
fixmahal <- function (dat, n=nrow(as.matrix(dat)), p=ncol(as.matrix(dat)),
method="fuzzy", cgen="fixed",
ca=NA, ca2=NA,
calpha=ifelse(method=="fuzzy",0.95,0.99),
calpha2=0.995,
pointit=TRUE, subset=n,
nc1=100+20*p,
startn = 18+p, mnc = floor(startn/2),
mer=ifelse(pointit,0.1,0), distcut=0.85,
maxit=5*n, iter=n*1e-5,
init.group=list(),
ind.storage=TRUE, countmode=100,
plot="none"){
# ca2, calpha2 are ignored unless method="fuzzy"
# Initializations, WDC
# print("fixmahal")
if (startn>n) startn <- n
dat <- as.matrix(dat)
if (method=="fuzzy"){
if (is.na(ca)){
ca <- qchisq(calpha,p)
ca2 <- qchisq(calpha2,p)
}
}
else{
if (cgen=="fixed"){
if (is.na(ca))
ca <- qchisq(calpha,p)
else
calpha <- pchisq(ca,p)
}
else{
if (is.na(ca))
ca <- qchisq(calpha,p)
ca <- cmahal(n,p,nmin=startn,cmin=ca,nc1=nc1,q=1)
}
ca2 <- 0
}
dat <- dat[1:n,1:p,drop=FALSE]
tsc <- 0 # number of too small FPCs
ncoll <- 0 # number of iterations leading to collinear data
gv <- rep (1, times=n)
# print("fpmi")
fpc <- fpmi(dat, n, p, gv, ca, ca2, method, plot, maxit, iter) # WDC-iteration
# print("1Ende")
if (cgen=="auto")
cv <- c(fpc$ca)
imatrix <- c(sum(fpc$g)) # matrix of FPC size and intersections
smatrix <- c(sum(fpc$g)) # matrix of FPC size and similarities
nc <- 1 # number of FPCs
clist <- list(fpc$mg) # list of FPC locations
vlist <- list(fpc$covg) # list of FPC covariance matrices
nfound <- c(1) # times found per FPC
expectratio <- c() # nfound/size
if (ind.storage)
glist <- list(fpc$g) # list of FPC indicator vectors
# Standard iterations
if (pointit){
if(subset < n)
itpoints <- sample(1:n,subset)
else{
itpoints <- 1:n
subset <- n
}
if (method=="fuzzy" || method=="classical" || method=="ml")
covall <- cov(dat)
else
covall <- cov.rob(dat, method=method)$cov
for (ni in 1:subset){
i <- itpoints[ni]
if(countmode*round(ni/countmode)==i)
cat("Iteration run no. ", ni, " of ",subset, "\n")
# print(covall)
gv <- mahalconf(dat,i,startn,covall=covall, plot)
# cat("i= ",i," gv= ",(1:366)[gv]," \n")
# print("iteration")
# cat(sum(gv[1:50]),sum(gv[51:100]))
fpc <- fpmi(dat, n, p, gv, ca, ca2, method, plot, maxit, iter)
# print(fpc$mg)
# print("neu?")
neu <- TRUE
j <- 1
while(j<=nc) {
if (sum((clist[[j]]-fpc$mg)^2)<1e-6 && sum((vlist[[j]]-fpc$covg)^2)<1e-5){
neu <- FALSE
cnum <- j
} # if j found
j <- j+1
} # while j
ng <- sum(fpc$g)
if (neu & (ng>=mnc)){
nc <- nc +1
nfound[nc] <- 1
clist[[nc]] <- fpc$mg
vlist[[nc]] <- fpc$covg
if (cgen=="auto")
cv <- c(cv,fpc$ca)
# cat("ind.storage= ",ind.storage, "nc= ", nc, "\n")
if (ind.storage){
glist[[nc]] <- fpc$g
j <- 1
while( j<=(nc-1)) {
imatrix[sseg(nc,j)] <- sum(pmin(fpc$g,glist[[j]]))
# cat(" (storage T) smatrix: ",nc,j,smatrix[sseg(nc,j)],"\n")
j<- j+1
} # for j
} # if ind.storage
else{
for(j in 1:(nc-1)){
if (method=="fuzzy"){
mah <- mahalanobis(dat,center=clist[[j]],
cov=solvecov(vlist[[j]])$inv,inverted=TRUE)
glistj <- wfu(mah, ca, ca2)
}
else{
if (cgen=="fixed"){
glistj <- (mahalanobis(dat,center=clist[[j]],
cov=solvecov(vlist[[j]])$inv,inverted=TRUE) <= ca)
}
else{
glistj <- (mahalanobis(dat,center=clist[[j]],
cov=solvecov(vlist[[j]])$inv,inverted=TRUE) <= cv[j])
}
}
imatrix[sseg(nc,j)] <- sum(pmin(fpc$g,glistj))
# cat(" smatrix: ",nc,j,smatrix[sseg(nc,j)],"\n")
} # for j
} # else (!ind.storage)
imatrix[sseg(nc,nc)] <- sum(fpc$g)
# cat("Neu: Punkt ",i,", Cluster mit ",sum(fpc$g)," Punkten. \n")
} # if neu & ng>=mnc
if (ng<mnc)
tsc <- tsc+1
if (!neu){
nfound[cnum] <- nfound[cnum]+1
}
if (neu & fpc$coll)
ncoll <- ncoll+1
} # for i
} # if pointit
# Iterations with init.group
# print(length(init.group))
grfpc <- FALSE
if(length(init.group)>0){
i <- 1
grfpc <- rep(0, times=length(init.group))
while(i<=length(init.group)){
gv <- init.group[[i]]
# print(sum(gv))
# print(length(gv))
fpc <- fpmi(dat, n, p, gv, ca, ca2, method, plot, maxit, iter)
# print("fpc")
neu <- TRUE
j <- 1
while(j<=nc) {
if (sum((clist[[j]]-fpc$mg)^2)<1e-6 && sum((vlist[[j]]-fpc$covg)^2)<1e-5){
neu <- FALSE
cnum <- j
} # if j found
j <- j+1
} # while j
ng <- sum(fpc$g)
if (neu & (ng>=mnc)){
# print("neu! smatrix")
nc <- nc +1
nfound[nc] <- 1
clist[[nc]] <- fpc$mg
vlist[[nc]] <- fpc$covg
if (cgen=="auto")
cv <- c(cv,fpc$ca)
# cat("ind.storage= ",ind.storage, "nc= ", nc, "\n")
if (ind.storage){
glist[[nc]] <- fpc$g
j <- 1
while( j<=(nc-1)) {
imatrix[sseg(nc,j)] <- sum(pmin(fpc$g,glist[[j]]))
# cat(" (storage T) smatrix: ",nc,j,smatrix[sseg(nc,j)],"\n")
j<- j+1
} # for j
} # if ind.storage
else{
for(j in 1:(nc-1)){
if (method=="fuzzy"){
mah <- mahalanobis(dat,center=clist[[j]],
cov=solvecov(vlist[[j]])$inv,inverted=TRUE)
glistj <- wfu(mah, ca, ca2)
}
else{
if (cgen=="fixed"){
glistj <- (mahalanobis(dat,center=clist[[j]],
cov=solvecov(vlist[[j]])$inv,inverted=TRUE) <= ca)
}
else{
glistj <- (mahalanobis(dat,center=clist[[j]],
cov=solvecov(vlist[[j]])$inv,inverted=TRUE) <= cv[j])
}
}
imatrix[sseg(nc,j)] <- sum(pmin(fpc$g,glistj))
# cat(" smatrix: ",nc,j,smatrix[sseg(nc,j)],"\n")
} # for j
} # else (!ind.storage)
imatrix[sseg(nc,nc)] <- sum(fpc$g)
} # if neu & ng>=mnc
if (ng<mnc)
tsc <- tsc+1
if (!neu)
nfound[cnum] <- nfound[cnum]+1
else
cnum <- nc
if (neu & fpc$coll)
ncoll <- ncoll+1
grfpc[i] <- cnum
i <- i+1
} # while i
} # if init.group
# Find structures
# print(nc)
if(nc>1){
for (i in 1:(nc-1)){
smatrix[sseg(i,i)] <- imatrix[sseg(i,i)]
for (j in (i+1):nc){
smatrix[sseg(i,j)] <- 2*imatrix[sseg(i,j)]/(imatrix[sseg(i,i)]+
imatrix[sseg(j,j)])
# cat("smatrix: ",i,j,smatrix[sseg(i,j)],"\n")
} # for j
}
smatrix[sseg(nc,nc)] <- imatrix[sseg(nc,nc)]
comat <- matrix(nrow=nc,ncol=nc)
for (i in 1:(nc-1))
for(j in (i+1):nc)
comat[i,j] <- comat[j,i] <- (smatrix[sseg(i,j)]>=distcut)
} #if nc>1
else
comat <- matrix(1,1)
# print(comat)
struc <- con.comp(comat)
# print(struc)
stn <- max(struc)
rm(comat)
# print("description of structures")
tf <- rep(0, times=stn) # structure: times found
maxf <- rep(0, times=stn) # structure, rep. fpc: expectratio
sfpc <- rep(0, times=stn) # structure: representative fpc
ser <- rep(0, times=stn) # structure: expectation ratio
for(i in 1:nc){
expectratio[i] <- nfound[i] / imatrix[sseg(i,i)]
tf[struc[i]] <- tf[struc[i]] + nfound[i]
if(expectratio[i]==maxf[struc[i]])
if(imatrix[sseg(i,i)]<imatrix[sseg(sfpc[struc[i]],sfpc[struc[i]])])
sfpc[struc[i]] <- i
if(expectratio[i]>maxf[struc[i]]){
sfpc[struc[i]] <- i
maxf[struc[i]] <- expectratio[i]
}
} # for i
for (i in 1:stn)
ser[i] <- (tf[i]*n)/(imatrix[sseg(sfpc[i],sfpc[i])]*subset)
skc <- stn
stn <- sum(ser>=mer)
skc <- skc-stn
# print(stn)
# Output
cvec <- ca
if (cgen=="auto")
ca <- cv
if (!ind.storage)
glist=FALSE
out <- list(nc=nc, g=glist, means=clist, covs=vlist,
nfound=nfound, er=expectratio,
tsc=tsc, ncoll=ncoll, skc=skc, grto=grfpc,
imatrix=imatrix, smatrix=smatrix,
stn=stn, stfound=tf, ser=ser,
sfpc=sfpc, ssig=sfpc[ser>=mer], sto=order(-ser),
struc=struc,
n=n, p=p, method=method, cgen=cgen,
ca=ca, ca2=ca2, cvec=cvec, calpha=calpha,
pointit=pointit, subset=subset,
mnc=mnc, startn=startn, mer=mer, distcut=distcut)
class(out) <- "mfpc"
out
}
summary.mfpc <- function(object, ...){
# print("Beginn reducemahal")
clist <- list(0)
vlist <- list(0)
expectratio <- c()
tf <- c()
sn <- c()
ca <- object$ca
method <- object$method
strn <- object$stn
if (strn>0)
for(i in 1:strn){
# cat(i, object$coefs[[object$sfpc[object$sto[i]]]], "\n")
clist[[i]] <- object$means[[object$sfpc[object$sto[i]]]]
# cat(i, object$vars[[object$sfpc[object$sto[i]]]], "\n")
vlist[[i]] <- object$covs[[object$sfpc[object$sto[i]]]]
tf[i] <- object$stfound[object$sto[i]]
sn[i] <- object$imatrix[sseg(object$sfpc[object$sto[i]],object$sfpc[object$sto[i]])]
expectratio[i] <- object$ser[object$sto[i]]
if (object$cgen=="auto")
ca[i] <- object$ca[object$sfpc[object$sto[i]]]
}
tskip <- sum(object$nfound)-sum(tf)
if (strn>0)
sim <- simmatrix(object)
else
sim <- NULL
out <- list(means=clist, covs=vlist, stn=strn, stfound=tf, sn=sn,
ser=expectratio,
tskip=tskip, skc=object$skc, tsc=object$tsc, sim=sim,
ca=ca, ca2=object$ca2, calpha=object$calpha,
mer=object$mer, mnc=object$mnc,
method=method, cgen=object$cgen,
pointit=object$pointit)
class(out) <- "summary.mfpc"
out
}
# clusters: List of indicator vectors of representative
# fpcs of significant structures
# in order of stfound.
# Specify x if ind.storage in fixreg was F
fpclusters.mfpc <- function(object, dat=NA, ca=object$ca, p=object$p, ...){
glist <- list()
if (object$method=="fuzzy"){
ca <- object$ca
ca2 <- object$ca2
}
# cat("stn= ",object$stn,"\n")
if (object$stn>0)
for(i in 1:object$stn){
if (identical(object$g,FALSE)){
mc <- object$means[[ object$sfpc[object$sto[i]] ]]
if (p>1){
cm <- solvecov(object$covs[[ object$sfpc[object$sto[i]] ]])
if (object$method=="fuzzy")
glist[[i]] <- wfu(mahalanobis(dat,mc,cm$inv,inverted=TRUE), ca, ca2)
if (object$cgen=="fixed")
glist[[i]] <- (mahalanobis(dat,mc,cm$inv,inverted=TRUE)<=ca)
else
glist[[i]] <- (mahalanobis(dat,mc,cm$inv,inverted=TRUE)<=ca[object$sfpc[object$sto[i]]] )
}
else{
cm <- object$covs[[ object$sfpc[object$sto[i]] ]]
if (object$cgen=="fixed")
glist[[i]] <- as.vector((dat-mc)^2)<=as.double(ca*cm)
else
glist[[i]] <- as.vector((dat-mc)^2)<=as.double(ca[object$sfpc[object$sto[i]]]*cm)
}
} # if g==F
else
glist[[i]] <- object$g[[object$sfpc[object$sto[i]]]]
# print(i)
} # for i
else
warning("No FPC was found often enough!")
glist
}
# Bhattacharyya coordinates of FPC no. no
# bw: black/white
plot.mfpc <- function(x, dat, no, bw=FALSE,
main=c("Representative FPC No. ",no),
xlab=NULL, ylab=NULL,
pch=NULL, col=NULL, ...){
n <- x$n
p <- x$p
ca <- x$ca
sumobj <- summary(x)
if (x$cgen=="auto")
ca <- sumobj$ca[no]
mg <- sumobj$means[[no]]
covg <- sumobj$covs[[no]]
# print("clusters")
if (x$method=="fuzzy"){
g <- fpclusters(x,dat)[[no]]
gv <- as.integer(g>0.5)
gvp <- as.integer(g>0.99)+
3*as.integer(g>0.01 & g<=0.99)
}
else
gvp <- gv <- as.integer(fpclusters(x,dat)[[no]])
if (n==sum(gv)){
cat("Cluster equals whole dataset => No cluster plot is done.\n")
}
else{
if (is.null(pch))
pch <- if (bw) gvp+1 else 1
if (is.null(col))
col <- if (bw) 1 else gvp+1
# print("plot")
if (p==1){
if (is.null(xlab))
xlab <- "Index"
if (is.null(ylab))
ylab <- deparse(substitute(dat))
plot(dat,col=col, pch=pch, main=main, xlab=xlab, ylab=ylab, ...)
abline(c(mg,0))
abline(c(mg-sqrt(ca*covg),0), lty="dotted")
abline(c(mg+sqrt(ca*covg),0), lty="dotted")
}
if (p==2){
xy <- xy.coords(dat,NULL)
if (is.null(xlab))
xlab <- xy$xlab
if (is.null(ylab))
ylab <- xy$ylab
plot(dat,col=col, pch=pch, main=main, xlab=xlab, ylab=ylab, ...)
points(as.vector(mg[1]),as.vector(mg[2]),pch="M",
col=ifelse(bw, col, "blue"), ...)
cge <- eigen(covg, symmetric=TRUE)
poly <- rep(0, times=200)
dim(poly) <- c(100,2)
for(i in 1:100)
poly[i,] <- (mg + c(sin(i*2*pi/100), cos(i*2*pi/100)) %*%
(t(cge$vectors) * sqrt(cge$values))* sqrt(ca))
polygon(poly[,1], poly[,2], border="black")
}
if (p>2){
gc1 <- batcoord(dat,gv)
if (is.null(xlab))
xlab <- "Discriminant projection 1"
if (is.null(ylab))
ylab <- "Discriminant projection 2"
plot(gc1$proj[,1:2],col=col, pch=pch,
main=main, xlab=xlab, ylab=ylab, ...)
points(as.vector((mg %*% gc1$units)[1]),
as.vector((mg %*% gc1$units)[2]), pch="M",
col=ifelse(bw, col, "blue"), ...)
}
}
invisible()
}
print.mfpc <- function(x, ...){
cat("Mahalanobis Fixed Point Cluster object\n")
cat(x$stn," representative stable fixed point clusters\n")
cat(" of totally ",x$nc," found fixed point clusters.\n")
invisible(x)
}
# Summary output of mfpc objects; choose fpcobj <- summary(fpcobj) to
# replace fpcobj by its reduction
print.summary.mfpc <- function(x, maxnc=30, ...){
# print("Beginn Summary")
# print(class(x))
minnc <- min(x$stn,maxnc)
cat(" * Mahalanobis Fixed Point Clusters *\n\n")
cat("Often a clear cluster in the data leads to several similar FPCs.\n")
cat("The summary shows the representative FPCs of groups of similar FPCs.\n\n")
cat("Method ",x$method," was used.\n")
cat("Number of representative FPCs: ", x$stn, "\n\n")
cat("FPCs with less than ",x$mnc," points were skipped.\n")
if (x$mer>0)
cat("FPCs with ratio of times found to number of points less than ",
x$mer," were skipped.\n")
cat(x$tsc+x$tskip, " iteration runs led to ",x$skc+x$tsc," skipped clusters.\n")
if (x$stn>maxnc)
cat("Warning! Only ",maxnc," clusters are displayed. \nSpecify maxnc in call of summary.mfpc if you want enlarged display or \nrun fixmahal with larger ca, calpha, mnc or mer to get fewer clusters. \n")
if (x$method=="fuzzy")
cat(" Weight 1 for r^2<= ",x$ca," weight 0 for r^2> ",
x$ca2," \n")
if (x$cgen=="fixed")
cat(" Constant ca= ", x$ca,
" corresponding to alpha= ",x$calpha,"\n")
cat("\n")
if (x$stn==0)
cat("No FPC was found often enough.\n")
else{
for(i in 1:minnc){
cat(" FPC ",i, "\n")
cat(" Times found (group members): ",x$stfound[i], "\n")
if (x$pointit)
cat(" Ratio to size: ",x$ser[i], "\n")
cat(" Mean:\n")
print(x$means[[i]])
cat(" Covariance matrix:\n")
print(x$covs[[i]])
if (x$cgen=="auto" && x$method!="fuzzy")
cat(" Constant ca= ",x$ca[[i]],"\n")
cat(" Number of points (sum of weights): ",x$sn[[i]], "\n\n")
}
cat("Number of points (rounded weights) in intersection of representative FPCs\n")
sm <- rep(0,times=minnc^2)
dim(sm) <- c(minnc,minnc)
for(i in 1:minnc)
for(j in 1:minnc)
sm[i,j] <- round(x$sim[sseg(i,j)])
print(sm)
}
invisible(x)
}
Try the fpc package in your browser
Any scripts or data that you put into this service are public.
fpc documentation built on Sept. 24, 2024, 9:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions print.summary.mfpc print.mfpc plot.mfpc fpclusters.mfpc summary.mfpc mahalconf wfu cmahal solvecov simmatrix con.comp fpclusters
Documented in cmahal con.comp fpclusters fpclusters.mfpc mahalconf plot.mfpc print.mfpc print.summary.mfpc simmatrix solvecov summary.mfpc wfu
# .packageName <- "fpc"
# .First.lib <- function(lib,pkg){
# require(cluster)
# # if (as.numeric(R.version$major)<=1 & as.numeric(R.version$minor)<9){
# # require(lqs)
# # require(mva)
# # }
# require(MASS)
# }
#
#
# generic functions
#
fpclusters <- function(object, ...)
UseMethod("fpclusters")
#
# general/set distance utilities
#
#
# sseg= position of distance between i and j in smatrix/imatrix
sseg <- function (i,j){
if (i>=j)
out <- (i-1)*i/2+j
else
out <- (j-1)*j/2+i
out
}
# Connectivity components of a graph by depth-first search
# Input: Boolean coincidence matrix, output: Vector of component numbers
con.comp <- function(comat){
nc <- ncol(comat)
ccn <- rep(0, times=nc) # con.comp number for each point
fhist <- rep(FALSE, times=nc) # indicator if point had been under consideration
stn <- 0 # current cc number
pn <- 1 # point no. to which similar objects are looked for
while(pn>0){
stn <- stn+1
repeat{
sm <- 0 # smallest new point no.
ccn[pn] <- stn
fhist[pn] <- TRUE
if(nc>1)
{
for(i in 2:nc)
{
# cat(i, pn, ccn[i], comat[i,pn],"\n")
if((ccn[i]==0) & (comat[i,pn]))
ccn[i] <- stn
if ((sm==0) & (ccn[i]==stn) & (fhist[i]==FALSE))
sm <- i
} # for i
# cat("stn=", stn, "sm=", sm, "\n")
} # if nc>1
if (sm>0)
pn <- sm
else
break
} # repeat
# print("repeat terminated")
pn <- 0
i <- 2
while(i<=nc){
if(ccn[i]==0){
pn <- i
i <- nc
} # if
i <- i+1
} # while i
} # while pn>0 (stn-loop)
ccn
}
# Intersection matrix between significant fpc groups
simmatrix <- function(fpcobj){
sim <- c()
# print(stn)
for(i in 1:fpcobj$stn)
for(j in i:fpcobj$stn){
# cat(i,j,fpcobj$sfpc[fpcobj$sto[i]],fpcobj$sfpc[fpcobj$sto[j]])
sim[sseg(i,j)] <- fpcobj$imatrix[sseg(fpcobj$sfpc[fpcobj$sto[i]],fpcobj$sfpc[fpcobj$sto[j]])]
}
sim
}
#
# fixmahal utilities
#
# Weighted ML-covariance matrix
cov.wml <- function (x, wt = rep(1/nrow(x), nrow(x)),
cor = FALSE, center = TRUE)
{
if (is.data.frame(x))
x <- as.matrix(x)
else if (!is.matrix(x))
stop("x must be a matrix or a data frame")
if (!all(is.finite(x)))
stop("x must contain finite values only")
n <- nrow(x)
if (with.wt <- !missing(wt)) {
if (length(wt) != n)
stop("length of wt must equal the number of rows in x")
if (any(wt < 0) || (s <- sum(wt)) == 0)
stop("weights must be non-negative and not all zero")
wt <- wt/s
}
if (is.logical(center)) {
center <- if (center)
colSums(wt * x)
else 0
}
else {
if (length(center) != ncol(x))
stop("length of center must equal the number of columns in x")
}
x <- sqrt(wt) * sweep(x, 2, center)
cov <- (t(x) %*% x)
y <- list(cov = cov, center = center, n.obs = n)
if (with.wt)
y$wt <- wt
if (cor) {
sdinv <- diag(1/sqrt(diag(cov)), nrow(cov))
y$cor <- sdinv %*% cov %*% sdinv
}
y
}
# inversion of cov-matrices; if singular, Eigenvalues below 1/cmax are set
# to 1/cmax.
solvecov <- function(m, cmax=1e10){
options(show.error.messages = FALSE)
covinv <- try(solve(m))
if(!inherits(covinv,"try-error"))
coll=FALSE
else{
p <- nrow(m)
cove <- eigen(m, symmetric=TRUE)
coll <- TRUE
if (min(cove$values)<1/cmax){
covewi <- diag(p)
for (i in 1:p)
if (cove$values[i]<1/cmax) covewi[i,i] <- cmax
else covewi[i,i] <- 1/cove$values[i]
}
else covewi <- diag(1/cove$values,nrow=length(cove$values))
# cat("covewi ",covewi)
covinv <- cove$vectors %*% covewi %*% t(cove$vectors)
}
options(show.error.messages = TRUE)
out <- list(inv=covinv,coll=coll)
}
# Generation of ca by "recursion for convergence"
cmahal <- function(n,p,nmin,cmin,nc1,c1=cmin,q=1){
cs <- rep(c1,max(n,nc1))
for (i in (nc1-1):nmin)
cs[i] <- cs[i+1]+q*(cs[i+1]-p)/(i+1-q)
cs[1:(nmin-1)] <- rep(cs[nmin],nmin-1)
if (n>nc1)
for (i in (nc1+1):n){
cs[i] <- cs[i-1]-q*(cs[i-1]-p)/i
if (cs[i]<=cmin){
cs[i] <- cmin
i <- n
}
}
cs
}
wfu <- function(md,ca,ca2,a1=1/(ca-ca2),a0=-a1*ca2){
v <- rep(0,length(md))
v[md<=ca] <- 1
mind <- (md>ca) & (md<=ca2)
v[mind] <- a1*md[mind]+a0
v
}
# mahalanofix=vector of mahalanobis distances from all points to
# center of points indexed by gv (default=all)
mahalanofix <- function (x, n=nrow(as.matrix(x)), p=ncol(as.matrix(x)),
gv=rep(1, times=n),
cmax=1e+10, method="ml") {
# print("Mahalanofix")
gv <- as.logical(gv)
ng <- sum(gv)
xg <- x[gv,1:p,drop=FALSE]
if (method=="ml"){
mg <- colMeans(xg)
if (method=="ml")
covg <- (ng-1)*cov(xg)/ng
else
covg <- cov(xg)
}
else{
# if (as.numeric(R.version$major)<=1 & as.numeric(R.version$minor)<9)
# require(lqs)
# require(MASS)
grob <- cov.rob(xg, method=method)
mg <- grob$center
covg <- grob$cov
}
cm <- solvecov(covg,cmax=cmax)
# print("robmahal")
md <- mahalanobis(x,mg,cm$inv,inverted=TRUE)
# print("end robmahal")
out <- list(md=md, mg=mg, covg=covg, covinv=cm$inv, coll=cm$coll)
}
# vector of mahalanobis distances from all points to
# center of points indexed by "fuzzy" gv (default=all)
mahalanofuz <- function (x, n=nrow(as.matrix(x)), p=ncol(as.matrix(x)),
gv=rep(1, times=n),
cmax=1e+10) {
# print("Mahalanofix")
mg <- cov.wml(x,gv)
covg <- mg$cov
mg <- mg$center
cm <- solvecov(covg, cmax=cmax)
# print("Mahalanobis")
md <- mahalanobis(x,mg,cm$inv,inverted=TRUE)
# print("End Mahalanofix")
out <- list(md=md, mg=mg, covg=covg, covinv=cm$inv, coll=cm$coll)
}
# Start configuration for fixmahal, point no. no
# startn points, cov/center of whole data
mahalconf <- function(x, no, startn, covall, plot){
p <- ncol(x)
n <- nrow(x)
gv <- rep(FALSE, n)
# print("mahalconf")
# print(p)
# print(n)
# print(no)
# print(startn)
# print(covall)
om <- order(mahalanobis(x,center=x[no,],cov=solvecov(covall)$inv, inverted=TRUE))
gv[om[1:(p+1)]] <- TRUE
# cat("conf ",gv," ",om,"\n")
if (plot=="start" || plot=="both")
plot(x,col=1+gv+((1:n)==no),pch=1+gv+((1:n)==no))
if (startn>p+1)
for (pn in (p+2):startn){
om <- order(mahalanofix(x,n,p,gv)$md)
j <- 0
while (sum(gv)<pn){
j <- j+1
if (!gv[om[j]])
gv[om[j]] <- TRUE
}
}
# print(sum(gv))
gv
}
#
# fixmahal iteration
#
fpmi <- function (dat, n=nrow(as.matrix(dat)), p=ncol(as.matrix(dat)),
gv, ca, ca2, method="ml", plot,
maxit=5*n, iter=n*1e-6) {
gva <- rep(0, times=n)
change <- TRUE
coll <- FALSE
ir <- 0
while (change) {
ir <- ir+1
if (method=="fuzzy")
change <- sum((gv-gva)^2>iter)
else
change <- !identical(gv,gva)
if (ir>maxit){
change <- FALSE
warning("Maximum number of iterations exceeded.")
}
gva <- gv
# print("mahalanobis")
if (method=="fuzzy")
mfix <- mahalanofuz(x=dat, n=n, p=p, gv=gv)
else
mfix <- mahalanofix(x=dat, n=n, p=p, gv=gv, method=method)
# print(ca)
md <- mfix$md
mg <- mfix$mg
if (mfix$coll) coll=TRUE
covg <- mfix$covg
# print(md)
if (method=="fuzzy"){
gv <- wfu(md, ca, ca2)
# print(gv)
cn <- ca
# plot(x,col=1+gv)
}
else{
ng <- sum(gv)
if (length(ca)==1)
cn <- ca
else
cn <- ca[ng]
# print(ng)
# cat("Target function = ",
# det(covg*(ng-1)/ng)*prod(exp((2-ca)/(1:ng))), "\n")
gv <- md<=cn
}
# print(gv)
if (plot=="iteration" || plot=="both"){
if (method=="fuzzy")
plotgv <- as.integer(gv>0.99)+
3*as.integer(gv>0.01 & gv<=0.99)
else
plotgv <- gv
plot(dat,col=1+plotgv,pch=1+plotgv)
}
} # while change
out <- list(mg=mg, covg=covg, md=md, gv=gv, coll=coll, method=method, ca=cn)
out
}
#
# fixmahal main and output
#
fixmahal <- function (dat, n=nrow(as.matrix(dat)), p=ncol(as.matrix(dat)),
method="fuzzy", cgen="fixed",
ca=NA, ca2=NA,
calpha=ifelse(method=="fuzzy",0.95,0.99),
calpha2=0.995,
pointit=TRUE, subset=n,
nc1=100+20*p,
startn = 18+p, mnc = floor(startn/2),
mer=ifelse(pointit,0.1,0), distcut=0.85,
maxit=5*n, iter=n*1e-5,
init.group=list(),
ind.storage=TRUE, countmode=100,
plot="none"){
# ca2, calpha2 are ignored unless method="fuzzy"
# Initializations, WDC
# print("fixmahal")
if (startn>n) startn <- n
dat <- as.matrix(dat)
if (method=="fuzzy"){
if (is.na(ca)){
ca <- qchisq(calpha,p)
ca2 <- qchisq(calpha2,p)
}
}
else{
if (cgen=="fixed"){
if (is.na(ca))
ca <- qchisq(calpha,p)
else
calpha <- pchisq(ca,p)
}
else{
if (is.na(ca))
ca <- qchisq(calpha,p)
ca <- cmahal(n,p,nmin=startn,cmin=ca,nc1=nc1,q=1)
}
ca2 <- 0
}
dat <- dat[1:n,1:p,drop=FALSE]
tsc <- 0 # number of too small FPCs
ncoll <- 0 # number of iterations leading to collinear data
gv <- rep (1, times=n)
# print("fpmi")
fpc <- fpmi(dat, n, p, gv, ca, ca2, method, plot, maxit, iter) # WDC-iteration
# print("1Ende")
if (cgen=="auto")
cv <- c(fpc$ca)
imatrix <- c(sum(fpc$g)) # matrix of FPC size and intersections
smatrix <- c(sum(fpc$g)) # matrix of FPC size and similarities
nc <- 1 # number of FPCs
clist <- list(fpc$mg) # list of FPC locations
vlist <- list(fpc$covg) # list of FPC covariance matrices
nfound <- c(1) # times found per FPC
expectratio <- c() # nfound/size
if (ind.storage)
glist <- list(fpc$g) # list of FPC indicator vectors
# Standard iterations
if (pointit){
if(subset < n)
itpoints <- sample(1:n,subset)
else{
itpoints <- 1:n
subset <- n
}
if (method=="fuzzy" || method=="classical" || method=="ml")
covall <- cov(dat)
else
covall <- cov.rob(dat, method=method)$cov
for (ni in 1:subset){
i <- itpoints[ni]
if(countmode*round(ni/countmode)==i)
cat("Iteration run no. ", ni, " of ",subset, "\n")
# print(covall)
gv <- mahalconf(dat,i,startn,covall=covall, plot)
# cat("i= ",i," gv= ",(1:366)[gv]," \n")
# print("iteration")
# cat(sum(gv[1:50]),sum(gv[51:100]))
fpc <- fpmi(dat, n, p, gv, ca, ca2, method, plot, maxit, iter)
# print(fpc$mg)
# print("neu?")
neu <- TRUE
j <- 1
while(j<=nc) {
if (sum((clist[[j]]-fpc$mg)^2)<1e-6 && sum((vlist[[j]]-fpc$covg)^2)<1e-5){
neu <- FALSE
cnum <- j
} # if j found
j <- j+1
} # while j
ng <- sum(fpc$g)
if (neu & (ng>=mnc)){
nc <- nc +1
nfound[nc] <- 1
clist[[nc]] <- fpc$mg
vlist[[nc]] <- fpc$covg
if (cgen=="auto")
cv <- c(cv,fpc$ca)
# cat("ind.storage= ",ind.storage, "nc= ", nc, "\n")
if (ind.storage){
glist[[nc]] <- fpc$g
j <- 1
while( j<=(nc-1)) {
imatrix[sseg(nc,j)] <- sum(pmin(fpc$g,glist[[j]]))
# cat(" (storage T) smatrix: ",nc,j,smatrix[sseg(nc,j)],"\n")
j<- j+1
} # for j
} # if ind.storage
else{
for(j in 1:(nc-1)){
if (method=="fuzzy"){
mah <- mahalanobis(dat,center=clist[[j]],
cov=solvecov(vlist[[j]])$inv,inverted=TRUE)
glistj <- wfu(mah, ca, ca2)
}
else{
if (cgen=="fixed"){
glistj <- (mahalanobis(dat,center=clist[[j]],
cov=solvecov(vlist[[j]])$inv,inverted=TRUE) <= ca)
}
else{
glistj <- (mahalanobis(dat,center=clist[[j]],
cov=solvecov(vlist[[j]])$inv,inverted=TRUE) <= cv[j])
}
}
imatrix[sseg(nc,j)] <- sum(pmin(fpc$g,glistj))
# cat(" smatrix: ",nc,j,smatrix[sseg(nc,j)],"\n")
} # for j
} # else (!ind.storage)
imatrix[sseg(nc,nc)] <- sum(fpc$g)
# cat("Neu: Punkt ",i,", Cluster mit ",sum(fpc$g)," Punkten. \n")
} # if neu & ng>=mnc
if (ng<mnc)
tsc <- tsc+1
if (!neu){
nfound[cnum] <- nfound[cnum]+1
}
if (neu & fpc$coll)
ncoll <- ncoll+1
} # for i
} # if pointit
# Iterations with init.group
# print(length(init.group))
grfpc <- FALSE
if(length(init.group)>0){
i <- 1
grfpc <- rep(0, times=length(init.group))
while(i<=length(init.group)){
gv <- init.group[[i]]
# print(sum(gv))
# print(length(gv))
fpc <- fpmi(dat, n, p, gv, ca, ca2, method, plot, maxit, iter)
# print("fpc")
neu <- TRUE
j <- 1
while(j<=nc) {
if (sum((clist[[j]]-fpc$mg)^2)<1e-6 && sum((vlist[[j]]-fpc$covg)^2)<1e-5){
neu <- FALSE
cnum <- j
} # if j found
j <- j+1
} # while j
ng <- sum(fpc$g)
if (neu & (ng>=mnc)){
# print("neu! smatrix")
nc <- nc +1
nfound[nc] <- 1
clist[[nc]] <- fpc$mg
vlist[[nc]] <- fpc$covg
if (cgen=="auto")
cv <- c(cv,fpc$ca)
# cat("ind.storage= ",ind.storage, "nc= ", nc, "\n")
if (ind.storage){
glist[[nc]] <- fpc$g
j <- 1
while( j<=(nc-1)) {
imatrix[sseg(nc,j)] <- sum(pmin(fpc$g,glist[[j]]))
# cat(" (storage T) smatrix: ",nc,j,smatrix[sseg(nc,j)],"\n")
j<- j+1
} # for j
} # if ind.storage
else{
for(j in 1:(nc-1)){
if (method=="fuzzy"){
mah <- mahalanobis(dat,center=clist[[j]],
cov=solvecov(vlist[[j]])$inv,inverted=TRUE)
glistj <- wfu(mah, ca, ca2)
}
else{
if (cgen=="fixed"){
glistj <- (mahalanobis(dat,center=clist[[j]],
cov=solvecov(vlist[[j]])$inv,inverted=TRUE) <= ca)
}
else{
glistj <- (mahalanobis(dat,center=clist[[j]],
cov=solvecov(vlist[[j]])$inv,inverted=TRUE) <= cv[j])
}
}
imatrix[sseg(nc,j)] <- sum(pmin(fpc$g,glistj))
# cat(" smatrix: ",nc,j,smatrix[sseg(nc,j)],"\n")
} # for j
} # else (!ind.storage)
imatrix[sseg(nc,nc)] <- sum(fpc$g)
} # if neu & ng>=mnc
if (ng<mnc)
tsc <- tsc+1
if (!neu)
nfound[cnum] <- nfound[cnum]+1
else
cnum <- nc
if (neu & fpc$coll)
ncoll <- ncoll+1
grfpc[i] <- cnum
i <- i+1
} # while i
} # if init.group
# Find structures
# print(nc)
if(nc>1){
for (i in 1:(nc-1)){
smatrix[sseg(i,i)] <- imatrix[sseg(i,i)]
for (j in (i+1):nc){
smatrix[sseg(i,j)] <- 2*imatrix[sseg(i,j)]/(imatrix[sseg(i,i)]+
imatrix[sseg(j,j)])
# cat("smatrix: ",i,j,smatrix[sseg(i,j)],"\n")
} # for j
}
smatrix[sseg(nc,nc)] <- imatrix[sseg(nc,nc)]
comat <- matrix(nrow=nc,ncol=nc)
for (i in 1:(nc-1))
for(j in (i+1):nc)
comat[i,j] <- comat[j,i] <- (smatrix[sseg(i,j)]>=distcut)
} #if nc>1
else
comat <- matrix(1,1)
# print(comat)
struc <- con.comp(comat)
# print(struc)
stn <- max(struc)
rm(comat)
# print("description of structures")
tf <- rep(0, times=stn) # structure: times found
maxf <- rep(0, times=stn) # structure, rep. fpc: expectratio
sfpc <- rep(0, times=stn) # structure: representative fpc
ser <- rep(0, times=stn) # structure: expectation ratio
for(i in 1:nc){
expectratio[i] <- nfound[i] / imatrix[sseg(i,i)]
tf[struc[i]] <- tf[struc[i]] + nfound[i]
if(expectratio[i]==maxf[struc[i]])
if(imatrix[sseg(i,i)]<imatrix[sseg(sfpc[struc[i]],sfpc[struc[i]])])
sfpc[struc[i]] <- i
if(expectratio[i]>maxf[struc[i]]){
sfpc[struc[i]] <- i
maxf[struc[i]] <- expectratio[i]
}
} # for i
for (i in 1:stn)
ser[i] <- (tf[i]*n)/(imatrix[sseg(sfpc[i],sfpc[i])]*subset)
skc <- stn
stn <- sum(ser>=mer)
skc <- skc-stn
# print(stn)
# Output
cvec <- ca
if (cgen=="auto")
ca <- cv
if (!ind.storage)
glist=FALSE
out <- list(nc=nc, g=glist, means=clist, covs=vlist,
nfound=nfound, er=expectratio,
tsc=tsc, ncoll=ncoll, skc=skc, grto=grfpc,
imatrix=imatrix, smatrix=smatrix,
stn=stn, stfound=tf, ser=ser,
sfpc=sfpc, ssig=sfpc[ser>=mer], sto=order(-ser),
struc=struc,
n=n, p=p, method=method, cgen=cgen,
ca=ca, ca2=ca2, cvec=cvec, calpha=calpha,
pointit=pointit, subset=subset,
mnc=mnc, startn=startn, mer=mer, distcut=distcut)
class(out) <- "mfpc"
out
}
summary.mfpc <- function(object, ...){
# print("Beginn reducemahal")
clist <- list(0)
vlist <- list(0)
expectratio <- c()
tf <- c()
sn <- c()
ca <- object$ca
method <- object$method
strn <- object$stn
if (strn>0)
for(i in 1:strn){
# cat(i, object$coefs[[object$sfpc[object$sto[i]]]], "\n")
clist[[i]] <- object$means[[object$sfpc[object$sto[i]]]]
# cat(i, object$vars[[object$sfpc[object$sto[i]]]], "\n")
vlist[[i]] <- object$covs[[object$sfpc[object$sto[i]]]]
tf[i] <- object$stfound[object$sto[i]]
sn[i] <- object$imatrix[sseg(object$sfpc[object$sto[i]],object$sfpc[object$sto[i]])]
expectratio[i] <- object$ser[object$sto[i]]
if (object$cgen=="auto")
ca[i] <- object$ca[object$sfpc[object$sto[i]]]
}
tskip <- sum(object$nfound)-sum(tf)
if (strn>0)
sim <- simmatrix(object)
else
sim <- NULL
out <- list(means=clist, covs=vlist, stn=strn, stfound=tf, sn=sn,
ser=expectratio,
tskip=tskip, skc=object$skc, tsc=object$tsc, sim=sim,
ca=ca, ca2=object$ca2, calpha=object$calpha,
mer=object$mer, mnc=object$mnc,
method=method, cgen=object$cgen,
pointit=object$pointit)
class(out) <- "summary.mfpc"
out
}
# clusters: List of indicator vectors of representative
# fpcs of significant structures
# in order of stfound.
# Specify x if ind.storage in fixreg was F
fpclusters.mfpc <- function(object, dat=NA, ca=object$ca, p=object$p, ...){
glist <- list()
if (object$method=="fuzzy"){
ca <- object$ca
ca2 <- object$ca2
}
# cat("stn= ",object$stn,"\n")
if (object$stn>0)
for(i in 1:object$stn){
if (identical(object$g,FALSE)){
mc <- object$means[[ object$sfpc[object$sto[i]] ]]
if (p>1){
cm <- solvecov(object$covs[[ object$sfpc[object$sto[i]] ]])
if (object$method=="fuzzy")
glist[[i]] <- wfu(mahalanobis(dat,mc,cm$inv,inverted=TRUE), ca, ca2)
if (object$cgen=="fixed")
glist[[i]] <- (mahalanobis(dat,mc,cm$inv,inverted=TRUE)<=ca)
else
glist[[i]] <- (mahalanobis(dat,mc,cm$inv,inverted=TRUE)<=ca[object$sfpc[object$sto[i]]] )
}
else{
cm <- object$covs[[ object$sfpc[object$sto[i]] ]]
if (object$cgen=="fixed")
glist[[i]] <- as.vector((dat-mc)^2)<=as.double(ca*cm)
else
glist[[i]] <- as.vector((dat-mc)^2)<=as.double(ca[object$sfpc[object$sto[i]]]*cm)
}
} # if g==F
else
glist[[i]] <- object$g[[object$sfpc[object$sto[i]]]]
# print(i)
} # for i
else
warning("No FPC was found often enough!")
glist
}
# Bhattacharyya coordinates of FPC no. no
# bw: black/white
plot.mfpc <- function(x, dat, no, bw=FALSE,
main=c("Representative FPC No. ",no),
xlab=NULL, ylab=NULL,
pch=NULL, col=NULL, ...){
n <- x$n
p <- x$p
ca <- x$ca
sumobj <- summary(x)
if (x$cgen=="auto")
ca <- sumobj$ca[no]
mg <- sumobj$means[[no]]
covg <- sumobj$covs[[no]]
# print("clusters")
if (x$method=="fuzzy"){
g <- fpclusters(x,dat)[[no]]
gv <- as.integer(g>0.5)
gvp <- as.integer(g>0.99)+
3*as.integer(g>0.01 & g<=0.99)
}
else
gvp <- gv <- as.integer(fpclusters(x,dat)[[no]])
if (n==sum(gv)){
cat("Cluster equals whole dataset => No cluster plot is done.\n")
}
else{
if (is.null(pch))
pch <- if (bw) gvp+1 else 1
if (is.null(col))
col <- if (bw) 1 else gvp+1
# print("plot")
if (p==1){
if (is.null(xlab))
xlab <- "Index"
if (is.null(ylab))
ylab <- deparse(substitute(dat))
plot(dat,col=col, pch=pch, main=main, xlab=xlab, ylab=ylab, ...)
abline(c(mg,0))
abline(c(mg-sqrt(ca*covg),0), lty="dotted")
abline(c(mg+sqrt(ca*covg),0), lty="dotted")
}
if (p==2){
xy <- xy.coords(dat,NULL)
if (is.null(xlab))
xlab <- xy$xlab
if (is.null(ylab))
ylab <- xy$ylab
plot(dat,col=col, pch=pch, main=main, xlab=xlab, ylab=ylab, ...)
points(as.vector(mg[1]),as.vector(mg[2]),pch="M",
col=ifelse(bw, col, "blue"), ...)
cge <- eigen(covg, symmetric=TRUE)
poly <- rep(0, times=200)
dim(poly) <- c(100,2)
for(i in 1:100)
poly[i,] <- (mg + c(sin(i*2*pi/100), cos(i*2*pi/100)) %*%
(t(cge$vectors) * sqrt(cge$values))* sqrt(ca))
polygon(poly[,1], poly[,2], border="black")
}
if (p>2){
gc1 <- batcoord(dat,gv)
if (is.null(xlab))
xlab <- "Discriminant projection 1"
if (is.null(ylab))
ylab <- "Discriminant projection 2"
plot(gc1$proj[,1:2],col=col, pch=pch,
main=main, xlab=xlab, ylab=ylab, ...)
points(as.vector((mg %*% gc1$units)[1]),
as.vector((mg %*% gc1$units)[2]), pch="M",
col=ifelse(bw, col, "blue"), ...)
}
}
invisible()
}
print.mfpc <- function(x, ...){
cat("Mahalanobis Fixed Point Cluster object\n")
cat(x$stn," representative stable fixed point clusters\n")
cat(" of totally ",x$nc," found fixed point clusters.\n")
invisible(x)
}
# Summary output of mfpc objects; choose fpcobj <- summary(fpcobj) to
# replace fpcobj by its reduction
print.summary.mfpc <- function(x, maxnc=30, ...){
# print("Beginn Summary")
# print(class(x))
minnc <- min(x$stn,maxnc)
cat(" * Mahalanobis Fixed Point Clusters *\n\n")
cat("Often a clear cluster in the data leads to several similar FPCs.\n")
cat("The summary shows the representative FPCs of groups of similar FPCs.\n\n")
cat("Method ",x$method," was used.\n")
cat("Number of representative FPCs: ", x$stn, "\n\n")
cat("FPCs with less than ",x$mnc," points were skipped.\n")
if (x$mer>0)
cat("FPCs with ratio of times found to number of points less than ",
x$mer," were skipped.\n")
cat(x$tsc+x$tskip, " iteration runs led to ",x$skc+x$tsc," skipped clusters.\n")
if (x$stn>maxnc)
cat("Warning! Only ",maxnc," clusters are displayed. \nSpecify maxnc in call of summary.mfpc if you want enlarged display or \nrun fixmahal with larger ca, calpha, mnc or mer to get fewer clusters. \n")
if (x$method=="fuzzy")
cat(" Weight 1 for r^2<= ",x$ca," weight 0 for r^2> ",
x$ca2," \n")
if (x$cgen=="fixed")
cat(" Constant ca= ", x$ca,
" corresponding to alpha= ",x$calpha,"\n")
cat("\n")
if (x$stn==0)
cat("No FPC was found often enough.\n")
else{
for(i in 1:minnc){
cat(" FPC ",i, "\n")
cat(" Times found (group members): ",x$stfound[i], "\n")
if (x$pointit)
cat(" Ratio to size: ",x$ser[i], "\n")
cat(" Mean:\n")
print(x$means[[i]])
cat(" Covariance matrix:\n")
print(x$covs[[i]])
if (x$cgen=="auto" && x$method!="fuzzy")
cat(" Constant ca= ",x$ca[[i]],"\n")
cat(" Number of points (sum of weights): ",x$sn[[i]], "\n\n")
}
cat("Number of points (rounded weights) in intersection of representative FPCs\n")
sm <- rep(0,times=minnc^2)
dim(sm) <- c(minnc,minnc)
for(i in 1:minnc)
for(j in 1:minnc)
sm[i,j] <- round(x$sim[sseg(i,j)])
print(sm)
}
invisible(x)
}
Try the fpc package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.